IQN: An Incremental Quasi-Newton Method with Local Superlinear Convergence Rate

TL;DR

The paper introduces IQN, a stochastic quasi-Newton method that achieves local superlinear convergence for minimizing sums of smooth, strongly convex functions, bridging the gap between deterministic and stochastic quasi-Newton methods.

Contribution

IQN is the first stochastic quasi-Newton method proven to converge superlinearly locally, using aggregated information, Taylor expansions, and cyclic updates.

Findings

IQN converges superlinearly in a local neighborhood of the solution.

It reduces variance of stochastic approximations through memory and aggregation.

IQN's convergence rate bridges deterministic and stochastic quasi-Newton methods.

Abstract

The problem of minimizing an objective that can be written as the sum of a set of $n$ smooth and strongly convex functions is considered. The Incremental Quasi-Newton (IQN) method proposed here belongs to the family of stochastic and incremental methods that have a cost per iteration independent of $n$. IQN iterations are a stochastic version of BFGS iterations that use memory to reduce the variance of stochastic approximations. The convergence properties of IQN bridge a gap between deterministic and stochastic quasi-Newton methods. Deterministic quasi-Newton methods exploit the possibility of approximating the Newton step using objective gradient differences. They are appealing because they have a smaller computational cost per iteration relative to Newton's method and achieve a superlinear convergence rate under customary regularity assumptions. Stochastic quasi-Newton methods utilize stochastic gradient differences in lieu of actual gradient differences. This makes their computational cost per iteration independent of the number of objective functions $n$. However, existing stochastic quasi-Newton methods have sublinear or linear convergence at best. IQN is the first stochastic quasi-Newton method proven to converge superlinearly in a local neighborhood of the optimal solution. IQN differs from state-of-the-art incremental quasi-Newton methods in three aspects: (i) The use of aggregated information of variables, gradients, and quasi-Newton Hessian approximation matrices to reduce the noise of gradient and Hessian approximations. (ii) The approximation of each individual function by its Taylor's expansion in which the linear and quadratic terms are evaluated with respect to the same iterate. (iii) The use of a cyclic scheme to update the functions in lieu of a random selection routine. We use these fundamental properties of IQN to establish its local superlinear convergence rate.

Figures

Equations

$${\mathbf{x}}^{\ast }:=\underset{\mathbf{x}\in {\mathbb{R}}^p}{\mathrm{argmin}}f(\mathbf{x}):=\underset{\mathbf{x}\in {\mathbb{R}}^p}{\mathrm{argmin}}\frac{1}{n}\sum _{i=1}^n{f}_i(\mathbf{x}).$$

$${\mathbf{x}}^{t+1}={\mathbf{x}}^t-{\eta }^t({\mathbf{B}}^t{)}^{-1}\nabla f({\mathbf{x}}^t),$$

$${\mathbf{s}}^t:={\mathbf{x}}^{t+1}-{\mathbf{x}}^t,{\mathbf{y}}^t:=\nabla f({\mathbf{x}}^{t+1})-\nabla f({\mathbf{x}}^t).$$

$${\mathbf{B}}^{t+1}{\mathbf{s}}^t={\mathbf{y}}^t.$$

$${\mathbf{B}}^{t+1}={\mathbf{B}}^t+\frac{{\mathbf{y}}^t{{\mathbf{y}}^t}^T}{{{\mathbf{y}}^t}^T{\mathbf{s}}^t}-\frac{{\mathbf{B}}^t{\mathbf{s}}^t{{\mathbf{s}}^t}^T{\mathbf{B}}^t}{{{\mathbf{s}}^t}^T{\mathbf{B}}^t{\mathbf{s}}^t}.$$

$$\underset{t\to \infty }{\lim }\frac{\Vert ({\mathbf{B}}^t-{\nabla }^{2}f({\mathbf{x}}^{\ast })){\mathbf{s}}^t\Vert }{\Vert {\mathbf{s}}^t\Vert }=0,$$

$$f_i(\mathbf{x})\approx f_i({\mathbf{z}}_i^t)+\nabla f_i({\mathbf{z}}_i^t{)}^T(\mathbf{x}-{\mathbf{z}}_i^t)+\frac{1}{2}(\mathbf{x}-{\mathbf{z}}_i^t{)}^T{\nabla }^2{f}_i({\mathbf{z}}_i^t)(\mathbf{x}-{\mathbf{z}}_i^t).$$

$$f(\mathbf{x})\approx \frac{1}{n}\sum _{i=1}^n\left[f_i({\mathbf{z}}_i^t)+\nabla f_i({\mathbf{z}}_i^t{)}^T(\mathbf{x}-{\mathbf{z}}_i^t)+\frac{1}{2}(\mathbf{x}-{\mathbf{z}}_i^t{)}^T{\mathbf{B}}_i^t(\mathbf{x}-{\mathbf{z}}_i^t)\right].$$

$${\mathbf{x}}^{t+1}={\left(\frac{1}{n}\sum _{i=1}^n{\mathbf{B}}_i^t\right)}^{-1}\left[\frac{1}{n}\sum _{i=1}^n{\mathbf{B}}_i^t{\mathbf{z}}_i^t-\frac{1}{n}\sum _{i=1}^n\nabla f_i({\mathbf{z}}_i^t)\right].$$

$${\mathbf{x}}^{t+1}={\mathbf{x}}^t-{\left(\frac{1}{n}\sum _{i=1}^n{\mathbf{B}}_i^t\right)}^{-1}\frac{1}{n}\sum _{i=1}^n\nabla f_i({\mathbf{z}}_i^t).$$

$$f(\mathbf{x})\approx \frac{1}{n}\sum _{i=1}^n\left[f_i({\mathbf{z}}_i^t)+\nabla f_i({\mathbf{z}}_i^t{)}^T(\mathbf{x}-{\mathbf{z}}_i^t)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}^t{)}^T{\mathbf{B}}_i^t(\mathbf{x}-{\mathbf{x}}^t)\right].$$

$${\mathbf{z}}_{i_t}^{t+1}={\mathbf{x}}^{t+1},{\mathbf{z}}_i^{t+1}={\mathbf{z}}_i^t\text{ for all }i\ne i_t.$$

$$\nabla f_{i_t}({\mathbf{z}}_i^{t+1})=\nabla f_{i_t}({\mathbf{x}}^{t+1}),\nabla f_i({\mathbf{z}}_i^{t+1})=\nabla f_i({\mathbf{z}}_i^t)\text{ for all }i\ne i_t.$$

$${\mathbf{s}}_i^t:={\mathbf{z}}_i^{t+1}-{\mathbf{z}}_i^t,{\mathbf{y}}_i^t:=\nabla f_i({\mathbf{z}}_i^{t+1})-\nabla f_i({\mathbf{z}}_i^t),$$

$${\mathbf{B}}_i^{t+1}={\mathbf{B}}_i^t+\frac{{\mathbf{y}}_i^t{\mathbf{y}}_i^{tT}}{{\mathbf{y}}_i^{tT}{\mathbf{s}}_i^t}-\frac{{\mathbf{B}}_i^t{\mathbf{s}}_i^t{\mathbf{s}}_i^{tT}{\mathbf{B}}_i^t}{{\mathbf{s}}_i^{tT}{\mathbf{B}}_i^t{\mathbf{s}}_i^t},\text{ for}i=i_t.$$

$${\mathbf{x}}^{t+1}=({\tilde{\mathbf{B}}}^t{)}^{-1}\left({\mathbf{u}}^t-{\mathbf{g}}^t\right),$$

$${\tilde{\mathbf{B}}}^{t+1}$$

$${\mathbf{u}}^{t+1}$$

$${\mathbf{g}}^{t+1}$$

$${\tilde{\mathbf{B}}}^{t+1}={\tilde{\mathbf{B}}}^t+\frac{{\mathbf{y}}_{i_t}^t{\mathbf{y}}_{i_t}^{tT}}{{\mathbf{y}}_i^{tT}{\mathbf{s}}_t^{i_t}}-\frac{{\mathbf{B}}_{i_t}^t{\mathbf{s}}_{i_t}^t{\mathbf{s}}_{i_t}^{tT}{\mathbf{B}}_{i_t}^t}{{\mathbf{s}}_{i_t}^{tT}{\mathbf{B}}_{i_t}^t{\mathbf{s}}_{i_t}^t}.$$

$$({\tilde{\mathbf{B}}}^{t+1}{)}^{-1}$$

$${\mathbf{U}}^t=({\tilde{\mathbf{B}}}^t{)}^{-1}-\frac{({\tilde{\mathbf{B}}}^t{)}^{-1}{\mathbf{y}}_{i_t}^t{\mathbf{y}}_{i_t}^{tT}({\tilde{\mathbf{B}}}^t{)}^{-1}}{{\mathbf{y}}_{i_t}^{tT}{\mathbf{s}}_{i_t}^t+{\mathbf{y}}_{i_t}^{tT}({\tilde{\mathbf{B}}}^t{)}^{-1}{\mathbf{y}}_{i_t}^t}.$$

$$\mu \Vert \mathbf{x}-\hat{\mathbf{x}}{\Vert }^2\le (\nabla f_i(\mathbf{x})-\nabla f_i(\hat{\mathbf{x}}){)}^T(\mathbf{x}-\hat{\mathbf{x}})\le L\Vert \mathbf{x}-\hat{\mathbf{x}}{\Vert }^2.$$

$$\Vert {\nabla }^2{f}_i(\mathbf{x})-{\nabla }^2{f}_i(\hat{\mathbf{x}})\Vert \le \tilde{L}\Vert \mathbf{x}-\hat{\mathbf{x}}\Vert .$$

$$\Vert {\nabla }^2{f}_i(\tilde{\mathbf{x}})(\mathbf{x}-\hat{\mathbf{x}})-(\nabla f_i(\mathbf{x})-\nabla f_i(\hat{\mathbf{x}}))\Vert \le \tilde{L}\Vert \mathbf{x}-\hat{\mathbf{x}}\Vert \max \left\{\Vert \mathbf{x}-\tilde{\mathbf{x}}\Vert ,\Vert \hat{\mathbf{x}}-\tilde{\mathbf{x}}\Vert \right\}.$$

$${\Vert {\mathbf{B}}_i^{t+1}-{\nabla }^2{f}_i({\mathbf{x}}^{\ast })\Vert }_{\mathbf{M}}\le \left[(1-\alpha {{\theta }_i^t}^2{)}^{1/2}+{\alpha }_3{\sigma }_i^t\right]{\Vert {\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast })\Vert }_{\mathbf{M}}+{\alpha }_4{\sigma }_i^t,$$

$${\theta }_i^t=\frac{\Vert \mathbf{M}({\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast })){\mathbf{s}}_i^t\Vert }{\Vert {\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast }){\Vert }_{\mathbf{M}}\Vert {\mathbf{M}}^{-1}{\mathbf{s}}_i^t\Vert }\text{for}{\mathbf{B}}_i^t\ne {\nabla }^2{f}_i({\mathbf{x}}^{\ast }),{\theta }_i^t=0\text{for}{\mathbf{B}}_i^t={\nabla }^2{f}_i({\mathbf{x}}^{\ast }).$$

$$\Vert {\mathbf{B}}_i^{t+1}-{\nabla }^2{f}_i({\mathbf{x}}^{\ast }){\Vert }_{\mathbf{M}}\lesssim (1-\alpha {{\theta }_i^t}^2{)}^{1/2}\Vert {\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast }){\Vert }_{\mathbf{M}}.$$

$$\underset{t\to \infty }{\lim }\frac{\Vert ({\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast })){\mathbf{s}}_i^t\Vert }{\Vert {\mathbf{s}}_i^t\Vert }=0.$$

$$\Vert {\mathbf{x}}^{t+1}-{\mathbf{x}}^{\ast }\Vert \le \frac{\tilde{L}{\Gamma }^t}{n}\sum _{i=1}^n{\Vert {\mathbf{z}}_i^t-{\mathbf{x}}^{\ast }\Vert }^2+\frac{{\Gamma }^t}{n}\sum _{i=1}^n\Vert \left({\mathbf{B}}_i^t-{\nabla }^2{f}_i({\mathbf{x}}^{\ast })\right)\left({\mathbf{z}}_i^t-{\mathbf{x}}^{\ast }\right)\Vert ,$$

$$\Vert {\mathbf{x}}^t-{\mathbf{x}}^{\ast }\Vert \le r^{[\frac{t-1}{n}]+1}\Vert {\mathbf{x}}^0-{\mathbf{x}}^{\ast }\Vert .$$

Full text

{AMS}

90C06, 90C25, 90C30, 90C52

1 Introduction

We study large scale optimization problems with objective functions expressed as the sum of a set of components which arise often in application domains such as machine learning [3, 2, 30, 8], control [5, 7, 15], and wireless communications [27, 23, 24]. Formally, we consider a variable ${\mathbf{x}}\in{\mathbb{R}}^{p}$ and a function $f$ which is defined as the average of $n$ smooth and strongly convex functions labelled $f_{i}:\mathbb{R}^{p}\rightarrow\mathbb{R}$ for $i=1,\ldots,n$. We refer to individual functions $f_{i}$ as sample functions and to the total number of functions $n$ as the sample size. Our goal is to find the optimal argument ${\mathbf{x}}^{*}$ that solves the strongly convex program

[TABLE]

We restrict attention to cases where the component functions $f_{i}$ are strongly convex and their gradients are Lipschitz continuous. We further focus in problems where $n$ is large enough so as to warrant application of stochastic or iterative methods. Our goal is to propose an iterative quasi-Newton method to solve (1) which is shown to exhibit a local superlinear convergence rate. This is achieved while performing local iterations with a cost of order $\mathcal{O}(p^{2})$ independent of the number of samples $n$.

Setting temporarily aside the complications related to the number of component functions, the minimization of $f$ in (1) can be carried out using iterative descent algorithms. A simple solution is to use gradient descent (GD) which iteratively descends along gradient directions $\nabla f({\mathbf{x}})=(1/n)\sum_{i=1}^{n}\nabla f_{i}({\mathbf{x}})$. GD incurs a per iteration computational cost of order $\mathcal{O}(np)$ and is known to converge at a linear rate towards ${\mathbf{x}}^{*}$ under the hypotheses we have placed on $f$. Whether the linear convergence rate of GD is acceptable depends on the desired accuracy and on the condition number of $f$ which, when large, can make the convergence constant close to one. As one or both of these properties often limit the applicability of GD, classical alternatives to improve convergence rates have been developed. Newton’s method adapts to the curvature of the objective by computing Hessian inverses and converges at a quadratic rate in a local neighborhood of the optimal argument irrespective of the problem’s condition number. To achieve this quadratic convergence rate, we must evaluate and invert Hessians resulting in a per iteration cost of order $\mathcal{O}(np^{2}+p^{3})$. Quasi-Newton methods build on the idea of approximating the Newton step using first-order information of the objective function and exhibit local superlinear convergence [4, 22, 12]. An important feature of quasi-Newton methods is that they have a per iteration cost of order $\mathcal{O}(np+p^{2})$, where the term $\mathcal{O}(np)$ corresponds to the cost of gradient computation and the cost $\mathcal{O}(p^{2})$ indicates the computational complexity of updating the approximate Hessian inverse matrix.

The combination of a local superlinear convergence rate and the smaller computational cost per iteration relative to Newton – a reduction by a factor of $p$ operations per iteration – make quasi-Newton methods an appealing choice. In the context of optimization problems having the form in (1), quasi-Newton methods also have the advantage that curvature is estimated using gradient evaluations. To see why this is meaningful we must recall that the customary approach to avoid the $\mathcal{O}(np)$ computational cost of GD iterations is to replace gradients $\nabla f({\mathbf{x}})$ by their stochastic approximations $\nabla f_{i}({\mathbf{x}})$, which can be evaluated with a cost of order $\mathcal{O}(p)$. One can then think of using stochastic versions of these gradients to develop stochastic quasi-Newton methods with per iterations cost of order $\mathcal{O}(p+p^{2})$. This idea was demonstrated to be feasible in [29] which introduces a stochastic (online) version of the BFGS quasi-Newton method as well as a stochastic version of its limited memory variant. Although [29] provides numerical experiments illustrating significant improvements in convergence times relative to stochastic (S) GD, theoretical guarantees are not established.

The issue of proving convergence of stochastic quasi-Newton methods is tackled in [18] and [19]. In [18] the authors show that stochastic BFGS may not be convergent because the Hessian approximation matrices can become close to singular. A regularized stochastic BFGS (RES) method is proposed by changing the proximity condition of BFGS to ensure that the eigenvalues of the Hessian inverse approximation are uniformly bounded. Enforcing this property yields a provably convergent algorithm. In [19] the authors show that the limited memory version of stochastic (online) BFGS proposed in [29] is almost surely convergent and has a sublinear convergence rate in expectation. This is achieved without using regularizations. An alternative provably convergent stochastic quasi-Newton method is proposed in [6]. This method differs from those in [29, 18, 19] in that it collects (stochastic) second order information to estimate the objective’s curvature. This is in contrast to estimating curvature using the difference of two consecutive stochastic gradients.

Although the methods in [29, 18, 19, 6] are successful in expanding the application of quasi-Newton methods to stochastic settings, their convergence rate is sublinear. This is not better than the convergence rate of SGD and, as is also the case in SGD, is a consequence of the stochastic approximation noise which necessitates the use of diminishing stepsizes. The stochastic quasi-Newton methods in [16, 20] resolve this issue by using the variance reduction technique proposed in [13]. The fundamental idea of the work in [13] is to reduce the noise of the stochastic gradient approximation by computing the exact gradient in an outer loop to use it in an inner loop for gradient approximation. The methods in [16, 20], which incorporate the variance reduction scheme presented in [13] into the update of quasi-Newton methods, are successful in achieving a linear convergence rate.

At this point, we must remark on an interesting mismatch. The convergence rate of SGD is sublinear, and the convergence rate of deterministic GD is linear. The use of variance reduction techniques in SGD recovers the linear convergence rate of GD, [13]. On the other hand, the convergence rate of stochastic quasi-Newton methods is sublinear, and the convergence rate of deterministic quasi-Newton methods is superlinear. The use of variance reduction in stochastic quasi-Newton methods achieves linear convergence but does not recover a superlinear rate. Hence, a fundamental question remains unanswered: Is it possible to design an incremental quasi-Newton method that recovers the superlinear convergence rate of deterministic quasi-Newton algorithms? In this paper, we show that the answer to this open problem is positive by proposing an incremental quasi-Newton method (IQN) with a local superlinear convergence rate. This is the first quasi-Newton method to achieve superlinear convergence while having a per iteration cost independent of the number of functions $n$ – the cost per iteration is of order $\mathcal{O}(p^{2})$.

There are three major differences between the IQN method and state-of-the-art incremental (stochastic) quasi-Newton methods that lead to the former’s superlinear convergence rate. First, the proposed IQN method uses the aggregated information of variables, gradients, and Hessian approximation matrices to reduce the noise of approximation for both gradients and Hessian approximation matrices. This is different to the variance-reduced stochastic quasi-Newton methods in [16, 20] that attempt to reduce only the noise of gradient approximations. Second, in IQN the index of the updated function is chosen in a cyclic fashion, rather than the random selection scheme used in the incremental methods in [29, 18, 19, 6]. The cyclic routine in IQN allows to bound the error at each iteration as a function of the errors of the last $n$ iterates, something that is not possible when using a random scheme. To explain the third and most important difference we point out that the form of quasi-Newton updates is the solution of a local second order Taylor approximation of the objective. It is possible to understand stochastic quasi-Newton methods as an analogous approximation of individual sample functions. However, it turns out that the state-of-the-art stochastic quasi-methods evaluate the linear and quadratic terms of the Taylor’s expansion at different points yielding and inconsistent approximation (Remark 7). The IQN method utilizes a consistent Taylor series which yields a more involved update which we nonetheless show can be implemented with the same computational cost. These three properties together lead to an incremental quasi-Newton method with a local superlinear convergence rate.

We start the paper by recapping the BFGS quasi-Newton method and the Dennis-Moré condition which is sufficient and necessary to prove superlinear convergence rate of the BFGS method (Section 2). Then, we present the proposed Incremental Quasi-Newton method (IQN) as an incremental aggregated version of the traditional BFGS method (Section 3). We first explain the difference between the Taylor’s expansion used in IQN and state-of-the-art incremental (stochastic) quasi-Newton methods. Further, we explain the mechanism for aggregation of the functions information and the scheme for updating the stored information. Moreover, we present an efficient implementation of the proposed IQN method with a computational complexity of the order $\mathcal{O}(p^{2})$ (Section 3.1). The convergence analysis of the IQN method is then presented (Section 4). We use the classic analysis of quasi-Newton methods to show that in a local neighborhood of the optimal solution the sequence of variables converges to the optimal argument ${\mathbf{x}}^{*}$ linearly after each pass over the set of functions (Lemma 3). We use this result to show that for each component function $f_{i}$ the Dennis-Moré condition holds (Proposition 4.7). However, this condition is not sufficient to prove superlinear convergence of the sequence of errors $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$, since it does not guarantee the Dennis-Moré condition for the global objective $f$. To overcome this issue we introduce a novel convergence analysis approach which exploits the local linear convergence of IQN to present a more general version of the Dennis-Moré condition for each component function $f_{i}$ (Lemma 4.9). We exploit this result to establish superlinear convergence of the iterates generated by IQN (Theorem 4.11). In Section 6, we present numerical simulation results, comparing the performance of IQN to that of first-order incremental and stochastic methods. We test the performance on a set of large-scale regression problems and observe strong numerical gain in total computation time relative to existing methods.

1.1 Notation

Vectors are written as lowercase ${\mathbf{x}}\in{\mathbb{R}}^{p}$ and matrices as uppercase ${\mathbf{A}}\in{\mathbb{R}}^{p\times p}$. We use $\|{\mathbf{x}}\|$ and $\|{\mathbf{A}}\|$ to denote the Euclidean norm of vector ${\mathbf{x}}$ and matrix ${\mathbf{A}}$, respectively. Given a positive definite matrix ${\mathbf{M}}$, the weighted matrix norm $\|{\mathbf{A}}\|_{\mathbf{M}}$ is defined as $\|{\mathbf{A}}\|_{\mathbf{M}}:=\|{\mathbf{M}}{\mathbf{A}}{\mathbf{M}}\|_{\mathbf{F}}$, where $\|.\|_{\mathbf{F}}$ is the Frobenius norm. Given a function $f$ its gradient and Hessian at point ${\mathbf{x}}$ are denoted as $\nabla f({\mathbf{x}})$ and $\nabla^{2}f({\mathbf{x}})$, respectively.

2 BFGS Quasi-Newton Method

Consider the problem in (1) for relatively large $n$. In a conventional optimization setting, this can be solved using a quasi-Newton method that iteratively updates a variable ${\mathbf{x}}^{t}$ for $t=0,1,\ldots$ based on the general recursive expression

[TABLE]

where $\eta^{t}$ is a scalar stepsize and ${{\mathbf{B}}^{t}}$ is a positive definite matrix that approximates the exact Hessian of the objective function $\nabla^{2}f({\mathbf{x}}^{t})$. The stepsize $\eta^{t}$ is evaluated based on a line search routine for the global convergence of quasi-Newton methods. Our focus in this paper, however, is on the local convergence of quasi-Newton methods, which requires the unit stepsize $\eta^{t}=1$. Therefore, throughout the paper we assume that the variable ${\mathbf{x}}^{t}$ is close to the optimal solution ${\mathbf{x}}^{*}$ – we will formalize the notion of being close to the optimal solution – and the stepsize is $\eta^{t}=1$.

The goal of quasi-Newton methods is to compute the Hessian approximation matrix ${{\mathbf{B}}^{t}}$ and its inverse ${({\mathbf{B}}^{t})}^{-1}$ by using only the first-order information, i.e., gradients, of the objective. Their use is widespread due to the many applications in which the Hessian information required in Newton’s method is either unavailable or computationally intensive. There are various approaches to approximate the Hessian, but the common feature among quasi-Newton methods is that the Hessian approximation must satisfy the secant condition. To be more precise, consider ${\mathbf{s}}^{t}$ and ${\mathbf{y}}^{t}$ as the variable and gradient variations, explicitly defined as

[TABLE]

Then, given the variable variation ${\mathbf{s}}^{t}$ and gradient variation ${\mathbf{y}}^{t}$, the Hessian approximation matrix in all quasi-Newton methods must satisfy the secant condition

[TABLE]

This condition is fundamental in quasi-Newton methods because the exact Hessian $\nabla^{2}f({\mathbf{x}}^{t})$ satisfies this equality when the iterates ${\mathbf{x}}^{t+1}$ and ${\mathbf{x}}^{t}$ are close to each other. If we consider the matrix ${\mathbf{B}}^{t+1}$ as the unknown matrix, the system of equations in (4) does not have a unique solution. Different quasi-Newton methods enforce different conditions on the matrix ${\mathbf{B}}^{t+1}$ to come up with a unique update. This extra condition is typically a proximity condition that ensures that ${\mathbf{B}}^{t+1}$ is close to the previous Hessian approximation matrix ${\mathbf{B}}^{t}$ [4, 22, 12]. In particular, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method defines the update of Hessian approximation matrix as

[TABLE]

The BFGS method is popular not only for its strong numerical performance relative to the gradient descent method, but also because it is shown to exhibit a superlinear convergence rate [4], thereby providing a theoretical guarantee of superior performance. In fact, it can be shown that, the BFGS update satisfies the condition

[TABLE]

known as the Dennis-Moré condition, which is both necessary and sufficient for superlinear convergence [12]. This result solidifies quasi-Newton methods as a strong alternative to first order methods when exact second-order information is unavailable. However, implementation of the BFGS method is not feasible when the number of functions $n$ is large, due to its high computational complexity on the order $\mathcal{O}(np+p^{2})$. In the following section, we propose a novel incremental BFGS method that has the computational complexity of $\mathcal{O}(p^{2})$ per iteration and converges at a superlinear rate.

3 IQN: Incremental aggregated BFGS

We propose an incremental aggregated BFGS algorithm, which we call the Incremental Quasi-Newton (IQN) method. The IQN method is incremental in that, at each iteration, only the information associated with a single function $f_{i}$ is updated. The particular function is chosen by cyclicly iterating through the $n$ functions. The IQN method is aggregated in that the aggregate of the most recently observed information of all functions $f_{1},\dots,f_{n}$ is used to compute the updated variable ${\mathbf{x}}^{t+1}$.

In the proposed method, we consider ${\mathbf{z}}_{1}^{t},\ldots,{\mathbf{z}}_{n}^{t}$ as the copies of the variable ${\mathbf{x}}$ at time $t$ associated with the functions $f_{1},\dots,f_{n}$, respectively. Likewise, define $\nabla f_{i}({\mathbf{z}}_{i}^{t})$ as the gradient corresponding to the $i$-th function. Further, consider ${\mathbf{B}}_{i}^{t}$ as a positive definite matrix which approximates the $i$-th component Hessian $\nabla^{2}f_{i}({\mathbf{x}}^{t})$. We refer to ${\mathbf{z}}_{i}^{t}$, $\nabla f_{i}({\mathbf{z}}_{i}^{t})$, and ${\mathbf{B}}_{i}^{t}$ as the information corresponding to the $i$-th function $f_{i}$ at step $t$. Note that the functions’ information is stored in a shared memory as shown in Fig. 1. To introduce the IQN method, we first explain the mechanism for computing the updated variable ${\mathbf{x}}^{t+1}$ using the stored information $\{{\mathbf{z}}_{i}^{t},\nabla f_{i}({\mathbf{z}}_{i}^{t}),{\mathbf{B}}_{i}^{t}\}_{i=1}^{n}$. Then, we elaborate on the scheme for updating the information of the functions.

To derive the full variable update, consider the second order approximation of the objective function $f_{i}({\mathbf{x}})$ centered around its current iterate ${\mathbf{z}}_{i}^{t}$,

[TABLE]

As in traditional quasi-Newton methods, we replace the $i$-th Hessian $\nabla^{2}f_{i}({\mathbf{z}}_{i}^{t})$ by ${\mathbf{B}}_{i}^{t}$. Using the approximation matrices in place of Hessians, the complete (aggregate) function $f({\mathbf{x}})$ can be approximated with

[TABLE]

Note that the right hand side of (8) is a quadratic approximation of the function $f$ based on the available information at step $t$. Hence, the updated iterate ${\mathbf{x}}^{t+1}$ can be defined as the minimizer of the quadratic program in (8), explicitly given by

[TABLE]

First note that the update in (9) shows that the updated variable ${\mathbf{x}}^{t+1}$ is a function of the stored information of all functions $f_{1},\dots,f_{n}$. Furthermore, we use the aggregated information of variables, gradients, and the quasi-Newton Hessian approximations to evaluate the updated variable. This is done to vanish the noise in approximating both gradients and Hessians as the sequence approaches the optimal argument.

Remark 1**.**

Given the BFGS Hessian approximation matrices $\{{\mathbf{B}}_{i}^{t}\}_{i=1}^{n}$ and gradients $\{\nabla f_{i}({\mathbf{z}}_{i}^{t})\}_{i=1}^{n}$, one may consider an update more akin to traditional descent-based methods, i.e.

[TABLE]

To evaluate the advantage of the proposed update for IQN in (9) relative to the update in (10), we proceed to study the Taylor’s expansion that leads to the update in (10). It can be shown that the update in (10) is the outcome of the following approximation

[TABLE]

*Observe that the linear term in (11) is centered at ${\mathbf{z}}_{i}^{t}$, while the quadratic term is approximated near the iterate ${\mathbf{x}}^{t}$. This inconsistency in the Taylor’s expansion of each function $f_{i}$ leads to an inaccurate second-order approximation, and subsequently a slower incremental quasi-Newton method. *

Thus far we have discussed the procedure to compute the updated variable ${\mathbf{x}}^{t+1}$ given the local iterates, gradients, and Hessian approximations at time $t$. Now, it remains to show how we update the local information of functions $f_{1},\dots,f_{n}$ using the variable ${\mathbf{x}}^{t+1}$. In each iteration of the IQN method, we update the local information of only a single function, chosen in a cyclic manner. Defining $i_{t}$ to be the index of the function selected at time $t$, we update the local variables ${\mathbf{z}}_{i_{t}}^{t+1}$, $\nabla f_{i_{t}}({\mathbf{z}}_{i}^{t+1})$, and ${\mathbf{B}}_{i}^{t+1}$ using the updated variable ${\mathbf{x}}^{t+1}$ while all other local variables remain unchanged. In particular, the variables ${\mathbf{z}}_{i}$ are updated as

[TABLE]

Observe in the update in (12) that the variable associated with the function $f_{i_{t}}$ is set to be the updated variable ${\mathbf{x}}^{t+1}$ while the other iterates are simply kept as their previous value. Likewise, we update the table of gradients accordingly with the gradient of $f_{i_{t}}$ evaluated at the new variable ${\mathbf{x}}^{t+1}$. The rest of gradients stored in the memory will stay unchanged, i.e.,

[TABLE]

To update the curvature information, it would be ideal to compute the Hessian $\nabla^{2}f_{i_{t}}({\mathbf{x}}^{t+1})$ and update the curvature information following the schemes for variables in (12) and gradients in (13). However, our focus is on the applications that the computation of the Hessian is either impossible or computationally expensive. Hence, to the update curvature approximation matrix ${\mathbf{B}}_{i_{t}}^{t}$ corresponding to the function $f_{i_{t}}$, we use the steps of BFGS in (5). To do so, we define variable and gradient variations associated with each individual function $f_{i}$ as

[TABLE]

respectively. The Hessian approximation ${\mathbf{B}}_{i_{t}}^{t}$ corresponding to the function $f_{i_{t}}$ can be computed using the update of BFGS as

[TABLE]

Again, the Hessian approximation matrices for all other functions remain unchanged, i.e., ${\mathbf{B}}_{i}^{t+1}={\mathbf{B}}_{i}^{t}$ for $i\neq i_{t}$. The system of updates in (12)-(15) explains the mechanism of updating the information of the function $f_{i_{t}}$ at step $t$. Notice that to update the Hessian approximation matrix for the $i_{t}$-th function there is no need to store the variations in (14), since the old variables ${\mathbf{z}}_{i}^{t}$ and $\nabla f_{i}({\mathbf{z}}_{i}^{t})$ are available in memory and the updated versions ${\mathbf{z}}_{i}^{t+1}={\mathbf{x}}^{t+1}$ and $\nabla f_{i}({\mathbf{z}}_{i}^{t+1})=\nabla f_{i}({\mathbf{x}}^{t+1})$ are evaluated at step $t$; see Fig. 1 for more details.

Because of the cyclic update scheme, the set of iterates $\{{\mathbf{z}}_{1}^{t},{\mathbf{z}}_{2}^{t},\ldots,{\mathbf{z}}_{n}^{t}\}$ is equal to the set $\{{\mathbf{x}}^{t},{\mathbf{x}}^{t-1},\ldots,{\mathbf{x}}^{t-n+1}\}$, and, therefore, the set of variables used in the update of IQN is the set of the last $n$ iterates. The update of IQN in (9) incorporates the information of all the functions $f_{1},\dots,f_{n}$ to compute the updated variable ${\mathbf{x}}^{t+1}$; however, it uses delayed variables, gradients, and Hessian approximations rather than the the updated variable ${\mathbf{x}}^{t+1}$ for all functions as in classic quasi-Newton methods. The use of delay allows IQN to update the information of a single function at each iteration, thus reducing the computational complexity relative to classic quasi-Newton methods.

Although the update in (9) is helpful in understanding the rationale behind the IQN method, it cannot be implemented at a low computation cost, since it requires computation of the sums $\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t}$, $\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t}{\mathbf{z}}_{i}^{t}$, and $\sum_{i=1}^{n}\nabla f_{i}({\mathbf{z}}_{i}^{t})$ as well as computing the inversion $(\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t})^{-1}$. In the following section, we introduce an efficient implementation of the IQN method that has the computational complexity of $\mathcal{O}(p^{2})$.

3.1 Efficient implementation of IQN

To see that the updating scheme in (9) requires evaluation of only a single gradient and Hessian approximation matrix per iteration, consider writing the update as

[TABLE]

where we define ${\tilde{\mathbf{B}}}^{t}:=\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t}$ as the aggregate Hessian approximation, ${\mathbf{u}}^{t}:=\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t}{\mathbf{z}}_{i}^{t}$ as the aggregate Hessian-variable product, and ${\mathbf{g}}^{t}:=\sum_{i=1}^{n}\nabla f_{i}({\mathbf{z}}_{i}^{t})$ as the aggregate gradient. Then, given that at step $t$ only a single index $i_{t}$ is updated, we can evaluate these variables for step $t+1$ as

[TABLE]

Thus, only ${\mathbf{B}}_{i_{t}}^{t+1}$ and $\nabla f_{i_{t}}({\mathbf{z}}_{i_{t}}^{t+1})$ are required to be computed at step $t$.

Although the updates in (17)-(19) have low computational complexity, the update in (16) requires computing $({\tilde{\mathbf{B}}}^{t})^{-1}$ which has a computational complexity of ${\mathcal{O}}(p^{3})$. This inversion can be avoided by simplifying the update in (17) as

[TABLE]

To derive the expression in (20) we have substituted the difference ${\mathbf{B}}_{i_{t}}^{t+1}-{\mathbf{B}}_{i_{t}}^{t}$ by its rank two expression in (15). Given the matrix $({\tilde{\mathbf{B}}}^{t})^{-1}$, by applying the Sherman-Morrison formula twice to the update in (20) we can compute $({\tilde{\mathbf{B}}}^{t+1})^{-1}$ as

[TABLE]

where the matrix ${\mathbf{U}}^{t}$ is evaluated as

[TABLE]

The computational complexity of the updates in (21) and (22) is of the order ${\mathcal{O}}(p^{2})$ rather than the ${\mathcal{O}}(p^{3})$ cost of computing the inverse directly. Therefore, the overall cost of IQN is of the order ${\mathcal{O}}(p^{2})$ which is substantially lower than ${\mathcal{O}}(np^{2})$ of deterministic quasi-Newton methods.

The complete IQN algorithm is outlined in Algorithm 1. Beginning with initial variable ${\mathbf{x}}^{0}$ and gradient and Hessian estimates $\nabla f_{i}({\mathbf{x}}^{0})$ and ${\mathbf{B}}_{i}^{0}$ for all $i$, each variable copy ${\mathbf{z}}^{0}_{i}$ is set to ${\mathbf{x}}^{0}$ in Step 1 and initial values are set for ${\mathbf{u}}^{0}$, ${\mathbf{g}}^{0}$ and $({\tilde{\mathbf{B}}}^{0})^{-1}$ in Step 2. For all $t$, in Step 4 the index $i_{t}$ of the next function to update is selected cyclically. The variable ${\mathbf{x}}^{t+1}$ is computed according to the update in (16) in Step 5. In Step 6, the variable ${\mathbf{s}}_{i_{t}}^{t+1}$ and gradient ${\mathbf{y}}_{i_{t}}^{t+1}$ variations are evaluated as in (14) to compute the BFGS matrix ${\mathbf{B}}_{i_{t}}^{t+1}$ from the update in (15). This information, as well as the updated variable and its gradient, are used in Step 7 to update ${\mathbf{u}}^{t+1}$ and ${\mathbf{g}}^{t+1}$ as in (18) and (19), respectively. The inverse matrix $({\tilde{\mathbf{B}}}^{t+1})^{-1}$ is also computed by following the expressions in (21) and (22). Finally in Step 8, we update the variable, gradient, and Hessian approximation tables based on the policies in (12), (13), and (15), respectively.

4 Convergence Analysis

In this section, we study the convergence rate of the proposed IQN method. We first establish its local linear convergence rate, then demonstrate limit properties of the Hessian approximations, and finally show that in a region local to the optimal point the sequence of residuals converges at a superlinear rate. To prove these results we make two main assumptions, both of which are standard in the analysis of quasi-Newton methods.

Assumption 1**.**

There exist positive constants $0<\mu\leq L$ such that, for all $i$ and ${\mathbf{x}},{\hat{\mathbf{x}}}\in{\mathbb{R}}^{p}$, we can write

[TABLE]

Assumption 2**.**

There exists a positive constant $0<\tilde{L}$ such that, for all $i$ and ${\mathbf{x}},{\hat{\mathbf{x}}}\in{\mathbb{R}}^{p}$, we can write

[TABLE]

The lower bound in (23) implies that the functions $f_{i}$ are strongly convex with constant $\mu$, and the upper bound shows that the gradients $\nabla f_{i}$ are Lipschitz continuous with parameter $L$.

The condition in Assumption 2, states that the Hessians $\nabla^{2}f_{i}$ are Lipschitz continuous with constant $\tilde{L}$. This assumption is commonly made in the analyses of Newton’s method [21] and quasi-Newton algorithms [4, 22, 12]. According to Lemma 3.1 in [4], Lipschitz continuity of the Hessians with constant $\tilde{L}$ implies that for $i=1,\dots,n$ and arbitrary vectors ${\mathbf{x}},{\tilde{\mathbf{x}}},{\hat{\mathbf{x}}}\in{\mathbb{R}}^{p}$ we can write

[TABLE]

We use the inequality in (25) in the process of proving the convergence of IQN.

The goal of BFGS quasi-Newton methods is to approximate the objective function Hessian using the first-order information. Likewise, in the incremental BFGS method, we aim to show that the Hessian approximation matrices for all functions $f_{1},\dots,f_{n}$ are close to the exact Hessian. In the following lemma, we study the difference between the $i$-th optimal Hessian $\nabla^{2}f_{i}({\mathbf{x}}^{*})$ and its approximation ${\mathbf{B}}_{i}^{t}$ over time.

Lemma 4.1**.**

Consider the proposed IQN method in (9). Further, let $i$ be the index of the updated function at step $t$, i.e., $i=i_{t}$. Define the residual sequence for function $f_{i}$ as $\sigma_{i}^{t}:=\max\{\|{\mathbf{z}}_{i}^{t+1}-{\mathbf{x}}^{*}\|,\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|\}$ and set ${\mathbf{M}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$. If Assumptions 1 and 2 hold and the condition $\sigma_{i}^{t}<m/(3\tilde{L})$ is satisfied then

[TABLE]

where $\alpha,\alpha_{3}$, and $\alpha_{4}$ are some positive bounded constants and

[TABLE]

Proof 4.2**.**

*See Appendix A. *

The result in (26) establishes an upper bound for the weighted norm $\|{\mathbf{B}}_{i}^{t+1}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ with respect to its previous value $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ and the sequence $\sigma_{i}^{t}:=\max\{\|{\mathbf{z}}_{i}^{t+1}-{\mathbf{x}}^{*}\|,\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|\}$, when the variables are in a neighborhood of the optimal solution such that $\sigma_{i}^{t}<m/(3\tilde{L})$. Indeed, the result in (26) holds only for the index $i=i_{t}$ and for the rest of indices we have $\|{\mathbf{B}}_{i}^{t+1}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}=\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ simply by definition of the cyclic update. Note that if the residual sequence $\sigma_{i}^{t}$ associated with $f_{i}$ approaches zero, we can simplify (26) as

[TABLE]

The equation in (28) implies that if ${\theta_{i}^{t}}$ is always strictly larger than zero, the sequence $\|{\mathbf{B}}_{i}^{t+1}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ approaches zero. If not, then the sequence ${\theta_{i}^{t}}$ converges to zero which implies the Dennis-Moré condition from (6), i.e.

[TABLE]

Therefore, under both conditions the result in (29) holds. This is true since the limit $\lim_{t\to\infty}\|{\mathbf{B}}_{i}^{t+1}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}=0$ yields the result in (29).

Based on this intuition, we proceed to show that the sequence $\sigma_{i}^{t}$ converges to zero for all $i=1,\dots,n$. To do so, we show that the sequence $\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ is linearly convergent for all $i=1,\dots,n$. To achieve this goal we first prove an upper bound for the error $\|{\mathbf{x}}^{t+1}-{\mathbf{x}}^{*}\|$ of IQN in the following lemma.

Lemma 4.3**.**

*Consider the proposed IQN method in (9). If the conditions in Assumptions 1 and 2 hold, then the sequence of iterates generated by IQN satisfies *

[TABLE]

*where $\Gamma^{t}:=\|((1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t})^{-1}\|$. *

Proof 4.4**.**

*See Appendix B. *

Lemma 4.3 shows that the residual $\|{\mathbf{x}}^{t+1}-{\mathbf{x}}^{*}\|$ is bounded above by a sum of quadratic and linear terms of the last $n$ residuals. This can eventually lead to a superlinear convergence rate by establishing the linear term converges to zero at a fast rate, leaving us with an upper bound of quadratic terms only. First, however, we establish a local linear convergence rate in the proceeding theorem to show that the sequence $\sigma_{i}^{t}$ converges to zero.

Lemma 4.5**.**

Consider the proposed IQN method in (9). If Assumptions 1 and 2 hold, then, for any $r\in(0,1)$ there are positive constants $\epsilon(r)$ and $\delta(r)$ such that if $\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|<\epsilon(r)$ and $\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}<\delta(r)$ for ${\mathbf{M}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$ and $i=1,2,\dots,n$, the sequence of iterates generated by IQN satisfies

[TABLE]

*Moreover, the sequences of norms $\{\|{\mathbf{B}}_{i}^{t}\|\}$ and $\{\|({{\mathbf{B}}_{i}^{t}})^{-1}\|\}$ are uniformly bounded. *

Proof 4.6**.**

*See Appendix C. *

The result in Lemma 4.5 shows that the sequence of iterates generated by IQN has a local linear convergence rate after each pass over all functions. Consequently, we obtain that the $i$-th residual sequence $\sigma_{i}^{t}$ is linearly convergent for all $i$. Note that Lemma 4.5 can be considered as an extension of Theorem 3.2 in [4] for incremental settings. Following the arguments in (28) and (29), we use the summability of the sequence $\sigma_{i}^{t}$ along with the result in Lemma 4.1 to prove Dennis-Moré condition for all functions $f_{i}$.

Proposition 4.7**.**

Consider the proposed IQN method in (9). Assume that the hypotheses in Lemmata 4.1 and 4.5 are satisfied. Then, for all $i=1,\dots,n$ it holds,

[TABLE]

Proof 4.8**.**

*See Appendix D. *

The statement in Proposition 4.7 indicates that for each function $f_{i}$ the Dennis-Moré condition holds. In the tradition quasi-Newton methods the Dennis-Moré condition is sufficient to show that the method is superlinearly convergent. However, the same argument does not hold for the proposed IQN method, since we can’t recover the Dennis-Moré condition for the global objective function $f$ from the result in Proposition 4.7. In other words, the result in (32) does not imply the limit in (6) required in the superlinear convergence analysis of quasi-Newton methods. Therefore, here we pursue a different approach and seek to prove that the linear terms $({\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*}))({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*})$ in (30) converge to zero at a superlinear rate, i.e., for all $i$ we can write $lim_{t\to\infty}\|({\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*}))({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*})\|/\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|=0$. If we establish this result, it follows from the result in Lemma 4.3 that the sequence of residuals $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ converges to zero superlinearly.

We continue the analysis of the proposed IQN method by establishing a generalized limit property that follows from the Dennis-Moré criterion in (6). In the following lemma, we leverage the local linear convergence of the iterates ${\mathbf{x}}^{t}$ to show that that the vector ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ lies in the null space of ${\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})$ as $t$ approaches infinity.

Lemma 4.9**.**

Consider the proposed IQN method in (9). Assume that the hypotheses in Lemmata 4.1 and 4.5 are satisfied. As $t$ approaches infinity, the following holds for all $i$,

[TABLE]

Proof 4.10**.**

*See Appendix E. *

The result in Lemma 4.9 can thus be used in conjunction with Lemma 4.3 to show that the residual $\|{\mathbf{x}}^{t+1}-{\mathbf{x}}^{*}\|$ is bounded by a sum of quadratic terms of previous residuals and a term that converges to zero superlinearly. This result leads us to the following result, namely the local superlinear convergence of the sequence of residuals with respect to the average sequence, stated in the following theorem.

Theorem 4.11**.**

Consider the proposed IQN method in (9). Suppose that the conditions in the hypotheses of Lemmata 4.1 and 4.5 are valid. Then, the sequence of residuals $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ satisfies

[TABLE]

Proof 4.12**.**

*See Appendix F. *

The result in (34) shows a mean-superlinear convergence rate for the sequence of iterates generated by IQN. To be more precise, it shows that the ratio that captures the error at step $t$ divided by the average of last $n$ errors converges to zero. This is not equivalent to the classic Q-superlinear convergence for full-batch quasi-Newton methods, i.e., $\lim_{t\to\infty}\|{\mathbf{x}}^{t+1}-{\mathbf{x}}^{*}\|/\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|=0$. Although Q-superlinear convergence of the residuals $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ is not provable, we can show that there exists a subsequence of the sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ that converges to zero superlinearly. In addition, there exists a superlinearly convergent sequence that is an upper bound for the original sequence of errors $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$. We formalize these results in the following theorem.

Theorem 4.13**.**

Consider the proposed IQN method in (9). Suppose that the conditions in the hypotheses of Lemmata 4.1 and 4.5 are valid. Then, there exists a subsequence of $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ that converges to zero superlinearly. Moreover, there exists a sequence $\zeta^{t}$ such that $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq\zeta^{t}$ for all $t\geq 0$, and the sequence $\zeta^{t}$ converges to zero at a superlinear rate, i.e.,

[TABLE]

Proof 4.14**.**

*See Appendix G. *

The first result in Theorem 4.13 states that although the whole sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ is not necessarily superlinearly convergent, there exists a subsequence of the sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ that converges at a superlinear rate. The second claim in Theorem 4.13 establishes R-superlinear convergence rate of the whole sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$. In other words, it guarantees that $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ is upper bounded by a superlinearly convergent sequence.

5 Related Works

Various methods have been studied in the literature to improve the performance of traditional full-batch optimization algorithms. The most famous method for reducing the computational complexity of gradient descent (GD) is stochastic gradient descent (SGD), which uses the gradient of a single randomly chosen function to approximate the full-gradient [2]. Incremental gradient descent method (IGD) is similar to SGD except the function is chosen in a cyclic routine [1]. Both SGD and IGD suffer from slow sublinear convergence rate because of the noise of gradient approximation. The incremental aggregated methods, which use memory to aggregate the gradients of all $n$ functions, are successful in reducing the noise of gradient approximation to achieve linear convergence rate [26, 28, 9, 13]. The work in [26] suggests a random selection of functions which leads to stochastic average gradient method (SAG), while the works in [1, 11, 17] use a cyclic scheme.

Moving beyond first order information, there have been stochastic quasi-Newton methods to approximate Hessian information [29, 18, 19, 20, 10]. All of these stochastic quasi-Newton methods reduce computational cost of quasi-Newton methods by updating only a randomly chosen single or small subset of gradients at each iteration. However, they are not able to recover the superlinear convergence rate of quasi-Newton methods [4, 22, 12]. The incremental Newton method (NIM) in [25] is the only incremental method shown to have a superlinear convergence rate; however, the Hessian function is not always available or computationally feasible. Moreover, the implementation of NIM requires computation of the incremental aggregated Hessian inverse which has the computational complexity of the order $\mathcal{O}(p^{3})$.

6 Numerical Results

We proceed by simulating the performance of IQN on a variety of machine learning problems on both artificial and real datasets. We compare the performance of IQN against a collection of well known first order stochastic and incremental algorithms—namely SAG, SAGA, and IAG. To begin, we look at a simple quadratic program, also equivalent to the solution of linear least squares estimation problem. Consider the objective function to be minimized,

[TABLE]

We generate $\mathbf{A}_{i}\in{\mathbb{R}}^{p\times p}$ as a random positive definite matrix and $\mathbf{b}_{i}\in{\mathbb{R}}^{p}$ as a random vector for all $i$. In particular we set the matrices $\mathbf{A}_{i}:=\text{diag}\{\mathbf{a}_{i}\}$ and generate random vectors ${\mathbf{a}}_{i}$ with the first $p/2$ elements chosen from $[1,10^{\xi/2}]$ and last $p/2$ elements chosen from $[10^{-\xi/2},1]$. The parameter $\xi$ is used to manually set the condition number for the quadratic program in (36), ranging from $\xi=1$ (i.e. small condition number $10^{2}$) and $\xi=2$ (i.e. large condition number $10^{4}$). The vectors $\mathbf{b}_{i}$ are chosen uniformly and randomly from the box $[0,10^{3}]^{p}$. The variable dimension is set to be $p=10$ and number of functions $n=1000$. Given that we focus on local convergence, we use a constant step size of $\eta=1$ for the proposed IQN method while choosing the largest step size allowable by the other methods to converge.

In Figure 2 we present a simulation of the convergence path of the normalized error $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|/\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$ for the quadratic program. In the the left image, we show a sample simulation path for all methods on the quadratic problem with a small condition number. Step sizes of $\eta=5\times 10^{-5}$, $\eta=10^{-4}$ and $\eta=10^{-6}$ were used for SAG, SAGA, and IAG, respectively. These step sizes are tuned to compare the best performance of these methods with IQN. The proposed method reaches a error of $10^{-10}$ after 10 passes through the data. Alternatively, SAGA achieves the same error of $10^{-5}$ after 30 passes, while SAG and IAG do not reach $10^{-5}$ after 40 passes.

In the right image of Figure 2, we repeat the same simulation but with larger condition number. In this case, SAG uses stepsize $\eta=2\times 10^{-4}$ while others remain the same. Observe that while the performance of IQN does not degrade with larger condition number, the first order methods all suffer large degradation. SAG, SAGA, and IAG reach after 40 passes a normalized error of $6.5\times 10^{-3}$, $5.5\times 10^{-2}$, and $9.6\times 10^{-1}$, respectively. It can be seen that IQN significantly outperforms the first order method for both condition number sizes, with the outperformance increasing for larger condition number. This is an expected result, as first order methods often do not perform well for ill conditioned problems.

6.1 Logistic regression

We proceed to numerically evaluate the performance of IQN relative to existing methods on the classification of handwritten digits in the MNIST database [14]. In particular, we solve the binary logistic regression problem. A logistic regression takes as inputs $n$ training feature vectors ${\mathbf{u}}_{i}\in{\mathbb{R}}^{p}$ with associated labels $v_{i}\in\{-1,1\}$ and outputs a linear classifier ${\mathbf{x}}$ to predict the label of unknown feature vectors. For the digit classification problem, each feature vector ${\mathbf{u}}_{i}$ represents a vectorized image and label $v_{i}$ its label as one of two digits. We evaluate for any training sample $i$ the probability of a label $v_{i}=1$ given image ${\mathbf{u}}_{i}$ as $P(v=1|{\mathbf{u}})=1/(1+\exp(-{\mathbf{u}}^{T}{\mathbf{x}}))$. The classifier ${\mathbf{x}}$ is chosen to be the vector which maximizes the log likelihood across all $n$ samples. Given $n$ images ${\mathbf{u}}_{i}$ with associated labels $v_{i}$, the optimization problem for logistic regression is written as

[TABLE]

where the first term is a regularization term parametrized by $\lambda\geq 0$.

For our simulations we select from the MNIST dataset $n=1000$ images with dimension $p=784$ labelled as one of the digits “0” or “8’ and fix the regularization parameter as $\lambda=1/n$ and stepsize $\eta=0.01$ for all first order methods. In Figure 3 we present the convergence path of IQN relative to existing methods in terms of the norm of the gradient. As in the case of the quadratic program, the IQN performs all gradient-based methods. IQN reaches a gradient magnitude of $4.8\times 10^{-8}$ after 60 passes through the data while the SAGA reaches only a magnitude of $7.4\times 10^{-5}$ (all other methods perform even worse). Further note that while the first order methods begin to level out after 60 passes, the IQN method continues to descend. These results demonstrate the effectiveness of IQN on a practical machine learning problem with real world data.

Appendix A Proof of Lemma 4.1

To prove the claim in Lemma 4.1, we first prove the the following lemma which is based on the result in [4, Lemma 5.2].

Lemma A.1**.**

Consider the proposed IQN method in (9). Let ${\mathbf{M}}$ be a nonsingular symmetric matrix such that

[TABLE]

for some $\beta\in[0,1/3]$ and vectors ${\mathbf{s}}_{i}^{t}$ and ${\mathbf{y}}_{i}^{t}$ in ${\mathbb{R}}^{p}$ with ${\mathbf{s}}_{i}^{t}\neq{\mathbf{0}}$. Consider $i$ as the index of the updated function at step $t$, i.e., $i=i_{t}$, and let ${\mathbf{B}}_{i}^{t}$ be symmetric and computed according to the update in (15). Then, there exist positive constants $\alpha$, $\alpha_{1}$, and $\alpha_{2}$ such that, for any symmetric ${\mathbf{A}}\in{\mathbb{R}}^{p\times p}$ we have,

[TABLE]

where $\alpha=(1-2\beta)/(1-\beta^{2})\in[3/8,1]$, $\alpha_{1}=2.5(1-\beta)^{-1}$, $\alpha_{2}=2(1+2\sqrt{p})\|{\mathbf{M}}\|_{\mathbf{F}}$, and

[TABLE]

Proof A.2**.**

*First note that the Hessian approximation ${\mathbf{B}}_{i}^{t+n}$ is equal to ${\mathbf{B}}_{i}^{t+1}$ if the function $f_{i}$ is updated at step $t$. Considering this observation and the result of Lemma 5.2. in [4] the claim in (39) follows. *

The result in Lemma A.1 provides an upper bound for the difference between the Hessian approximation matrix ${\mathbf{B}}_{i}^{t+n}$ and any positive definite matrix ${\mathbf{A}}$ with respect to the difference between the previous Hessian approximation ${\mathbf{B}}_{i}^{t}$ and the matrix ${\mathbf{A}}$. The interesting choice for the arbitrary matrix ${\mathbf{A}}$ is the Hessian of the $i$-th function at the optimal argument, i.e., ${\mathbf{A}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})$, which allows us to capture the difference between the sequence of Hessian approximation matrices for function $f_{i}$ and the Hessian $\nabla^{2}f_{i}({\mathbf{x}}^{*})$ at the optimal argument. We proceed to use the result in Lemma A.1 for ${\mathbf{M}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$ and ${\mathbf{A}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})$ to prove the claim in (26). To do so, we first need to show that the condition in (38) is satisfied. Note that according to the condition in Assumptions 1 and 2 we can write

[TABLE]

This observation implies that the left hand side of the condition in (38) for ${\mathbf{M}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$ is bounded above by

[TABLE]

Thus, the condition in (38) is satisfied since $\tilde{L}\sigma_{i}^{t}/m<1/3.$ Replacing the upper bounds in (41) and (42) into the expression in (39) implies the claim in (26) with

[TABLE]

and the proof is complete.

Appendix B Proof of Lemma 4.3

Start by subtracting ${\mathbf{x}}^{*}$ from both sides of (9) to obtain

[TABLE]

As the gradient of $f$ at the optimal point is the vector zero, i.e., $(1/n)\sum_{i=1}^{n}\nabla f_{i}({\mathbf{x}}^{*})=\mathbf{0}$, we can subtract $(1/n)\sum_{i=1}^{n}\nabla f_{i}({\mathbf{x}}^{*})$ from the right hand side of (44) and rearrange terms to obtain

[TABLE]

The expression in (45) relates the residual at time $t+1$ to the previous $n$ residuals and the Hessian approximations ${\mathbf{B}}_{i}^{t}$. To analyze this further, we can replace the Hessian approximations ${\mathbf{B}}_{i}^{t}$ with the actual Hessians $\nabla^{2}f_{i}({\mathbf{x}}^{*})$ and the approximation difference $\nabla^{2}f_{i}({\mathbf{x}}^{*})-{\mathbf{B}}_{i}^{t}$. To do so, we add and subtract $({1}/{n})\sum_{i=1}^{n}\nabla^{2}f_{i}({\mathbf{x}}^{*})\left({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\right)$ to the right hand side of (45) and rearrange terms to obtain

[TABLE]

We proceed to take the norms of both sides and use the triangle inequality to obtain an upper bound on the norm of the residual $\|{\mathbf{x}}^{t+1}-{\mathbf{x}}^{*}\|$,

[TABLE]

To obtain the quadratic term in (30) from the first term in (B), we use the Lipschitz continuity of the Hessians $\nabla^{2}f_{i}$ which leads to the inequality

[TABLE]

Replacing the expression $\|\nabla^{2}f_{i}({\mathbf{x}}^{*})\left({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\right)-\left(\nabla f_{i}({\mathbf{z}}_{i}^{t})-\nabla f_{i}({\mathbf{x}}^{*})\right)\|$ in (B) by the upper bound in (48), the claim in (30) follows.

Appendix C Proof of Lemma 4.5

In this proof we use some steps in the proof of [4, Theorem 3.2]. To start we use the fact that in a finite-dimensional vector space there always exists a constant $\eta>0$ such that $\|{\mathbf{A}}\|\leq\eta\|{\mathbf{A}}\|_{\mathbf{M}}$. Consider $\gamma=1/m$ is an upper bound for the norm $\|\nabla^{2}f({\mathbf{x}}^{*})^{-1}\|$. Assume that $\epsilon(r)=\epsilon$ and $\delta(r)=\delta$ are chosen such that

[TABLE]

Based on the assumption that $\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}\leq\delta$ we can derive the upper bound $\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq\eta\delta$. This observation along with the inequality $\|\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq L$ implies that $\|{\mathbf{B}}_{i}^{0}\|\leq\eta\delta+L$. Therefore, we obtain $\|(1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{0}\|\leq\eta\delta+L$. The second inequality in (49) implies that $2\gamma(1+r)\eta\delta\leq r$. Based on this observation and the inequalities $\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq\eta\delta<2\eta\delta$ and $\gamma\geq\|\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1}\|$, we obtain from Banach Lemma that $\|({\mathbf{B}}_{i}^{0})^{-1}\|\leq(1+r)\gamma$. Following the same argument for the matrix $((1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{0})^{-1}$ with the inequalities $\|(1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{0}-(1/n)\sum_{i=1}^{n}\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq(1/n)\sum_{i=1}^{n}\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq\eta\delta$ and $\|\nabla^{2}f({\mathbf{x}}^{*})^{-1}\|\leq\gamma$ we obtain that

[TABLE]

This upper bound in conjunction with the result in (30) yields

[TABLE]

Considering the assumptions that $\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|\leq\epsilon$ and $\|{\mathbf{B}}_{i}^{0}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq\eta\delta<2\eta\delta$ we can write

[TABLE]

where the second inequality follows from the second condition in (49). Without loss of generality, assume that $i_{0}=1$. Then, based on the result in (26) we obtain

[TABLE]

We proceed to the next iteration which leads to the inequality

[TABLE]

And since the updated index is $i_{1}=2$ we obtain

[TABLE]

With the same argument we can show that all $\left\|{\mathbf{B}}_{t}^{t}-\nabla^{2}f_{t}({\mathbf{x}}^{*})\right\|_{\mathbf{M}}\leq 2\delta$ and $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq\epsilon$, for all iterates $t=1,\dots,n$. Moreover, we have $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq r\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$ for $t=1,\dots,n$.

Now we use the results for iterates $t=1,\dots,n$ as the base of our induction argument. To be more precise, let’s assume that for iterates $t=jn+1,jn+2,\dots,jn+n$ we know that the residuals are bounded above by $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq r^{j+1}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$ and the Hessian approximation matrices ${\mathbf{B}}_{i}^{t}$ satisfy the inequalities $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq 2\eta\delta$. Our goal is to show that for iterates $t=(j+1)n+1,(j+1)n+2,\dots,(j+1)n+n$ the inequalities $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq r^{j+2}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$ and $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq 2\eta\delta$ hold.

Based on the inequalities $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq 2\eta\delta$ and $\|\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1}\|\leq\gamma$ we can show that for all $t=jn+1,jn+2,\dots,jn+n$ we have

[TABLE]

Using (56) and the inequality in (26) for the iterate $t=(j+1)n+1$, we obtain

[TABLE]

Since the variables are updated in a cyclic fashion the set of variables $\{{\mathbf{z}}_{i}^{(j+1)n}\}_{i=1}^{i=n}$ is equal to the set $\{{\mathbf{x}}^{(j+1)n-i}\}_{i=0}^{i=n-1}$. By considering this relation and replacing the norms $\|[{\mathbf{B}}_{i}^{(j+1)n}-\nabla^{2}f_{i}({\mathbf{x}}^{*})]({\mathbf{z}}_{i}^{(j+1)n}-{\mathbf{x}}^{*})\|$ by their upper bounds $2\eta\delta\|{\mathbf{z}}_{i}^{(j+1)n}-{\mathbf{x}}^{*}\|$ we can simplify the right hand side of (C) as

[TABLE]

Since $\|{\mathbf{x}}^{jn+i}-{\mathbf{x}}^{*}\|\leq\epsilon$ for all $j=1,\dots,n$, we obtain

[TABLE]

According to the second inequality in (49) and the assumption that for iterates $t=jn+1,jn+2,\dots,jn+n$ we know that $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq r^{j+1}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$, we can replace the right hand side of (59) by the following upper bound

[TABLE]

Now we show that the updated Hessian approximation ${\mathbf{B}}_{i_{t}}^{(j+1)n+1}$ for $t=(j+1)n+1$ satisfies the inequality $\|{\mathbf{B}}_{i_{t}}^{(j+1)n+1}-\nabla^{2}f_{i_{t}}({\mathbf{x}}^{*})\|_{\mathbf{M}}\leq 2\delta$. According to the result in (26), we can write

[TABLE]

Now observe that $\sigma_{i_{t}}^{jn+1}=\max\{\|{\mathbf{x}}^{(j+1)n+1}-{\mathbf{x}}^{*}\|,\|{\mathbf{x}}^{jn+1}-{\mathbf{x}}^{*}\|\}$ is bounded above by $r^{j+1}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$. Applying this substitution into (C) and considering the conditions $\|{\mathbf{B}}_{i_{t}}^{jn+1}-\nabla^{2}f_{i_{t}}({\mathbf{x}}^{*})\|_{\mathbf{M}}\leq 2\delta$ and $\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|\leq\epsilon$ lead to the inequality

[TABLE]

By writing the expression in (62) for previous iterations and using a recursive logic we obtain that

[TABLE]

Based on the first inequality in (49), the right hand side of (63) is bounded above by $\delta$. Moreover, the norm $\|{\mathbf{B}}_{i_{t}}^{0}-\nabla^{2}f_{i_{t}}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ is also upper bounded by $\delta$. These two bounds imply that

[TABLE]

and consequently $\|{\mathbf{B}}_{i_{t}}^{(j+1)n+1}-\nabla^{2}f_{i_{t}}({\mathbf{x}}^{*})\|\leq 2\eta\delta$. By following the steps from (C) to (64), we can show for all iterates $t=(j+1)n+1,(j+1)n+2,\dots,(j+1)n+n$ the inequalities $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq r^{j+2}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|$ and $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq 2\eta\delta$ hold. The induction proof is complete and (31) holds. Moreover, the inequality $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|\leq 2\eta\delta$ holds for all $i$ and steps $t$. Hence, the norms $\|{\mathbf{B}}_{i}^{t}\|$ and $\|({{\mathbf{B}}_{i}^{t}})^{-1}\|$, and consequently $\|(1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t}\|$ and $\|((1/n)\sum_{i=1}^{n}{\mathbf{B}}_{i}^{t})^{-1}\|$ are uniformly bounded.

Appendix D Proof of Proposition 4.7

According to the result in Lemma 4.5, we can show that the sequence of errors $\sigma_{i}^{t}=\max\{\|{\mathbf{z}}_{i}^{t+1}-{\mathbf{x}}^{*}\|,\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|\}$ is summable for all $i$. To do so, consider the sum of the sequence $\sigma_{i}^{t}$ which is upper bounded by

[TABLE]

Note that the last time that the index $i$ is chosen before time $t$ should be in the set $\{t-1,\dots,t-n\}$. This observation in association with the result in (31) implies that

[TABLE]

Simplifying the sum in the right hand side of (66) yields

[TABLE]

Thus, the sequence $\sigma_{i}^{t}$ is summable for all $i=1,\dots,n$. To complete the proof we use the following result from Lemma 3.3 in [12].

Lemma D.1**.**

Let $\{\phi^{t}\}$ and $\{\delta^{t}\}$ be sequences of nonnegative numbers such that

[TABLE]

*Then, the sequence $\{\phi^{t}\}$ converges. *

Considering the results in Lemmata 4.1 and D.1, and the fact that $\sigma_{i}^{t}$ is summable as shown in(67), we obtain that the sequence $\left\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\right\|_{\mathbf{M}}$ for ${\mathbf{M}}:=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$ is convergent and the following limit exists

[TABLE]

where $l$ is a nonnegative constant. Moreover, following the proof of Theorem 3.4 in [12] we can show that

[TABLE]

and, therefore, summing both sides implies,

[TABLE]

Replacing $\theta_{i}^{t}$ in (71) by its definition in (27) results in

[TABLE]

Since the norm $\|{\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})\|_{\mathbf{M}}$ is upper bounded and the eigenvalues of the matrix ${\mathbf{M}}=\nabla^{2}f_{i}({\mathbf{x}}^{*})^{-1/2}$ are uniformly lower and upper bounded, we conclude from the result in (72) that

[TABLE]

which yields the claim in (32).

Appendix E Proof of Lemma 4.9

Consider the sets of variable variations ${\mathcal{S}}_{1}=\{{\mathbf{s}}_{i}^{t+n\tau}\}_{\tau=0}^{\tau=T}$ and ${\mathcal{S}}_{2}=\{{\mathbf{s}}_{i}^{t+n\tau}\}_{\tau=0}^{\tau=\infty}$. It is trivial to show that ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ is in the span of the set ${\mathcal{S}}_{2}$, since the sequences of variables ${\mathbf{x}}^{t}$ and ${\mathbf{z}}_{i}^{t}$ converge to ${\mathbf{x}}^{*}$ and we can write ${\mathbf{x}}^{*}-{\mathbf{z}}_{i}^{t}=\sum_{\tau=0}^{\infty}{\mathbf{s}}_{i}^{t+n\tau}$. We proceed to show that the vector ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ is also in the span of the set ${\mathcal{S}}_{1}$ when $T$ is sufficiently large. To do so, we use a contradiction argument. Let’s assume that the vector ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ does not lie in the span of the set ${\mathcal{S}}_{1}$, and, therefore, it can be decomposed as the sum of two non-zero vectors given by

[TABLE]

where ${\mathbf{v}}^{t}_{\parallel}$ lies in the span of ${\mathcal{S}}_{1}$ and ${\mathbf{v}}^{t}_{\perp}$ is orthogonal to the span of ${\mathcal{S}}_{1}$. Since we assume that ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ does not lie in the span of ${\mathcal{S}}_{1}$, we obtain that ${\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}$ also does not lie in this span, since ${\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}$ can be written as the sum ${\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}={\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}+\sum_{\tau=0}^{T}{\mathbf{s}}_{i}^{t+n\tau}$. These observations imply that we can also decompose the vector ${\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}$ as

[TABLE]

where ${\mathbf{v}}^{t+nT}_{\parallel}$ lies in the span of ${\mathcal{S}}_{1}$ and ${\mathbf{v}}^{t+nT}_{\perp}$ is orthogonal to the span of ${\mathcal{S}}_{1}$. Moreover, we obtain that ${\mathbf{v}}^{t+nT}_{\perp}$ is equal to ${\mathbf{v}}^{t}_{\perp}$, i.e.,

[TABLE]

This is true since ${\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}$ can be written as the sum of ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ and a group of vectors that lie in the span of ${\mathcal{S}}_{1}$. We assume that the norm $\|{\mathbf{v}}^{t+nT}_{\perp}\|=\|{\mathbf{v}}^{t}_{\perp}\|=\epsilon$ where $\epsilon>0$ is a strictly positive constant. According to the linear convergence of the sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ in Lemma 4.5 we know that

[TABLE]

If we pick large enough $T$ such that $r^{T}\|{\mathbf{x}}^{0}-{\mathbf{x}}^{*}\|<\epsilon$, then we obtain $\|{\mathbf{z}}_{i}^{t+nT}-{\mathbf{x}}^{*}\|<\epsilon$ which contradicts the assumption $\|{\mathbf{v}}^{t}_{\perp}\|=\epsilon$. Thus, we obtain that the vector ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ is also in the span of set ${\mathcal{S}}_{1}$.

Since the vector ${\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}$ is in the span of ${\mathcal{S}}_{1}$, we can write the normalized vector $({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*})/\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ as a linear combination of the set of normalized vectors $\{{\mathbf{s}}_{i}^{t+n\tau}/\|{\mathbf{s}}_{i}^{t+n\tau}\|\}_{\tau=0}^{\tau=T}$. This property allows to write

[TABLE]

where $a_{\tau}$ is coefficient of the vector ${\mathbf{s}}_{i}^{t+n\tau}$ when we write $({\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*})/\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ as the linear combination of the normalized vectors $\{{\mathbf{s}}_{i}^{t+n\tau}/\|{\mathbf{s}}_{i}^{t+n\tau}\|\}_{\tau=0}^{\tau=T}$. Now since the index of the difference ${\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})$ does not match with the descent directions ${\mathbf{s}}_{i}^{t}+n\tau$. We add and subtract the term ${\mathbf{B}}_{i}^{t+n\tau}$ to the expression ${\mathbf{B}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*})$ and use the triangle inequality to write

[TABLE]

We first simplify the first limit in the right hand side of (E). Using the Cauchy-Schwarz inequality and the result in Proposition 4.7 we can write

[TABLE]

Based on the results in (E) and (E), to prove the claim in (33) it remains to show

[TABLE]

To reach this goal, we first study the limit of the difference between two consecutive update Hessian approximation matrices $\lim_{t\rightarrow\infty}\|{\mathbf{B}}_{i}^{t}-{\mathbf{B}}_{i}^{t+n}\|$. Note that if we set ${\mathbf{A}}={\mathbf{B}}_{i}^{t}$ in (39), we obtain that

[TABLE]

where ${\mathbf{M}}=(\nabla^{2}f_{i}({\mathbf{x}}^{*}))^{-1/2}$. By adding and subtracting the term $\nabla^{2}f_{i}({\mathbf{x}}^{*}){\mathbf{s}}_{i}^{t}$ and using the result in (32), we can show that the difference $\|{\mathbf{B}}_{i}^{t+n}-{\mathbf{B}}_{i}^{t}\|_{\mathbf{M}}$ approaches zero asymptotically. In particular,

[TABLE]

Since $\|{\mathbf{y}}_{i}^{t}-\nabla^{2}f_{i}({\mathbf{x}}^{*}){\mathbf{s}}_{i}^{t}\|$ is bounded above by $\tilde{L}\|{\mathbf{s}}_{i}^{t}\|\max\{\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|,\|{\mathbf{z}}_{i}^{t+1}-{\mathbf{x}}^{*}\|\}$ and the eigenvalues of the matrix ${\mathbf{M}}$ are uniformly bounded we obtain that the first limit in the right hand side of (E) converges to zero. Further, the result in (32) shows that the second limit in the right hand side of (E) also converges to zero. Therefore,

[TABLE]

Following the same argument we can show that for any two consecutive Hessian approximation matrices the difference approaches zero asymptotically. Thus, we obtain

[TABLE]

Observing the result in (E) we can show that

[TABLE]

Therefore, the result in (81) holds. The claim in (33) follows by combining the results in (E), (E), and (81).

Appendix F Proof of Theorem 4.11

The result in Lemma 4.3 implies

[TABLE]

Divide both sides of (87) by $(1/n)\sum_{i=1}^{n}\left\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\right\|$ to obtain

[TABLE]

Since the error $\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ is a lower bound for the sum of errors $\sum_{i=1}^{n}\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$, we can replace $\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ for $\sum_{i=1}^{n}\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\|$ into (88) which implies

[TABLE]

Since $\Gamma^{t}$ is bounded above, computing the limit of both sides in (F) yields

[TABLE]

The result in (90) in association with the simplification for the sum $\sum_{i=1}^{n}\left\|{\mathbf{z}}_{i}^{t}-{\mathbf{x}}^{*}\right\|=\sum_{i=0}^{n-1}\left\|{\mathbf{x}}^{t-i}-{\mathbf{x}}^{*}\right\|$ leads to the claim in (34).

Appendix G Proof of Theorem 4.13

Consider the definition of the sequence ${\tilde{\mathbf{x}}}^{t}=\operatornamewithlimits{argmax}_{u\in\{tn,\dots,tn+n-1\}}\{\|{\mathbf{x}}^{u}-{\mathbf{x}}^{*}\|\}$ which is a subsequence of the sequence $\{{\mathbf{x}}^{t}\}_{t=0}^{\infty}$. Our goal is to show this subsequence converges superlinearly to ${\mathbf{x}}^{*}$, i.e., $\lim_{t\to\infty}\frac{\|{\tilde{\mathbf{x}}}^{t+1}-{\mathbf{x}}^{*}\|}{\|{\tilde{\mathbf{x}}}^{t}-{\mathbf{x}}^{*}\|}=0.$ To do so, first note that the result in Theorem 4.11 implies that

[TABLE]

which follows from the inequality $\max\{\|{\mathbf{x}}^{t-1}-{\mathbf{x}}^{*}\|,\dots,\|{\mathbf{x}}^{t-n}-{\mathbf{x}}^{*}\|\}\geq(1/n)(\|{\mathbf{x}}^{t-1}-{\mathbf{x}}^{*}\|+\dots+\|{\mathbf{x}}^{t-n}-{\mathbf{x}}^{*}\|)$. Based on the limit in (91), there exists a large enough $t_{0}$ such that for all $t\geq t_{0}$ the following inequality holds,

[TABLE]

Combining the inequality in (92) with the inequalities $\|{\mathbf{x}}^{t-i}-{\mathbf{x}}^{*}\|\leq\max\{\|{\mathbf{x}}^{t-1}-{\mathbf{x}}^{*}\|,\dots,\|{\mathbf{x}}^{t-n}-{\mathbf{x}}^{*}\|\}$ for $i=1,\dots,n-1$ yields

[TABLE]

and consequently we can generalize this result to obtain

[TABLE]

for any positive integer $\tau$ such that $t-\tau\geq t_{0}$.

We use the result in (94) to build a superlinearly convergent subsequence of the residuals sequence $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$. If we define ${\mathbf{x}}^{tn+u_{t}^{*}}$ as the iterate that has the largest error among the iterates in the $t+1$-th pass, i.e.,

[TABLE]

then it follows that ${\tilde{\mathbf{x}}}^{t}={\mathbf{x}}^{tn+u_{t}^{*}}$, where $u_{t}^{*}\in\{0,1,\dots,n-1\}$. Moreover, we obtain

[TABLE]

The equality follows from the definition of the iterate ${\tilde{\mathbf{x}}}^{t}$ and the definition in (95), and the inequality holds because of the result in (94). Considering the result in (91), computing the limit of both sides leads to the conclusion that the sequence $\|{\tilde{\mathbf{x}}}^{t}-{\mathbf{x}}^{*}\|$ is superlinearly convergent. In other words, we obtain that the subsequence $\{\|{\mathbf{x}}^{tn+u_{t}^{*}}-{\mathbf{x}}^{*}\|\}_{t=0}^{t=\infty}$ superlinearly converges to zero.

Let’s define the sequence $q^{t}$ such that $q^{kn}=\dots=q^{kn+n-1}=\|{\tilde{\mathbf{x}}}^{k}-{\mathbf{x}}^{*}\|$ for $k=0,1,2,\dots$, which means that the value of the sequence $q^{t}$ is fixed for each pass and is equal to the max error of the corresponding pass. Therefore, it is trivial to show that $q^{t}$ is always larger than or equal to $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$, i.e., $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|\leq q^{t}$ for all $t\geq 0$. Now define the sequence $\zeta^{t}$ such that $\zeta^{t}=q^{t}$ for $t=0,\dots,n-1$, and for $t\geq n$

[TABLE]

According to this definition we can verify that $\zeta^{t}$ is an upper bound for the sequence $q^{t}$ and, consequently, an upper bound for the sequence of errors $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$. Based on the definition of the sequence $\zeta^{t}$ in (97), the ratio $\zeta^{t+1}/\zeta^{t}$ is given by $(q^{\left\lfloor\frac{t+1}{n}\right\rfloor n+n-1}/q^{\left\lfloor\frac{t+1}{n}\right\rfloor n-1})^{1/n}$. This simplification in association with the definitions of the sequences $\|{\mathbf{x}}^{t}-{\mathbf{x}}^{*}\|$ and $\|{\tilde{\mathbf{x}}}^{t}-{\mathbf{x}}^{*}\|$ implies that

[TABLE]

which leads to the claim in (35).