Quasi-Newton Methods for Machine Learning: Forget the Past, Just Sample

TL;DR

This paper introduces sampled quasi-Newton methods that use random sampling around the current point to build more reliable Hessian approximations, improving efficiency and convergence in machine learning optimization tasks.

Contribution

The paper proposes novel sampled LBFGS and sampled LSR1 algorithms that utilize local sampling for Hessian approximation, enhancing convergence and parallelization capabilities.

Findings

Outperform classical quasi-Newton methods in classification and neural network training.

Require fewer data points (epochs) for convergence.

Effective in parallel and distributed computing environments.

Abstract

We present two sampled quasi-Newton methods (sampled LBFGS and sampled LSR1) for solving empirical risk minimization problems that arise in machine learning. Contrary to the classical variants of these methods that sequentially build Hessian or inverse Hessian approximations as the optimization progresses, our proposed methods sample points randomly around the current iterate at every iteration to produce these approximations. As a result, the approximations constructed make use of more reliable (recent and local) information, and do not depend on past iterate information that could be significantly stale. Our proposed algorithms are efficient in terms of accessed data points (epochs) and have enough concurrency to take advantage of parallel/distributed computing environments. We provide convergence guarantees for our proposed methods. Numerical tests on a toy classification problem as well as on popular benchmarking binary classification and neural network training tasks reveal that the methods outperform their classical variants.

Figures

Tables

Table 1. Table 1: Deep Neural Networks used in the experiments.

model	$d$	input	# classes
LeNet	3.2M	$28\times 28\times 3$	10
\hdashlinealexnet v2	50.3M	$224\times 224\times 3$	1,000
\hdashlinevgg a	132.8M	$224\times 224\times 3$	1,000

Table 2. Table 2: Summary of Computational Cost and Storage (per iteration) for different Quasi-Newton methods.

method	computational cost	storage
BFGS	$nd+d^2+{\kappa }_{ls}nd$	$d^2$
\hdashlineLBFGS	$nd+4md+{\kappa }_{ls}nd$	$2md$
\hdashlineS-LBFGS	$nd+mnd+4md+{\kappa }_{ls}nd$	$2md$
SR1	$nd+d^2+nd+{\kappa }_{tr}{d}^2$	$d^2$
\hdashlineLSR1	$nd+nd+{\kappa }_{tr}md$	$2md$
\hdashlineS-LSR1	$nd+mnd+nd+{\kappa }_{tr}md$	$2md$

Table 3. Table 3: Problem details for empirical evaluation of probabilities.

problem	figure
$\lambda =(\underset{d/2\text{ -1’s}}{\underbrace{-1,\dots ,-1}},\underset{d/2\text{ 1’s}}{\underbrace{1,\dots ,1}})$	Figures 7a & 8a
\hdashline$\lambda ={10}^{-4}(\underset{d/2\text{ -1’s}}{\underbrace{-1,\dots ,-1}},\underset{d/2\text{ 1’s}}{\underbrace{1,\dots ,1}})$	Figures 7b & 8b
\hdashlineToy Problem: (small network)	Figures 7c & 8c
\hdashlineToy Problem: (medium network)	Figures 7d & 8d

Table 4. Table 4: Toy Classification Problem: Neural Network Details

network	structure	$d$
small	2-2-2-2-2-2	$36$
\hdashlinemedium	2-4-8-8-4-2	$176$
\hdashlinelarge	2-10-20-20-10-2	$908$

Table 5. Table 5: Details for MNIST and CIFAR10 Problems.

problem	network	structure	$d$
MNIST	Net1	$784-C_{5,3}-C_{5,5}-10-10$	$1378$
\hdashlineMNIST	Net2	$784-C_{5,6}-C_{5,16}-C_{4,120}-82-10$	$44164$
CIFAR10	Net3	$1024,3-C_{5,3}-C_{5,4}-10-10$	$1652$
\hdashlineCIFAR10	Net4	$1024,3-C_{5,3}-C_{5,5}-16-32-10$	$2312$

Table 6. Table 6: Summary of datasets.

dataset	# of samples	# of features	# of classes
ijcnn1	35,000	22	2
\hdashlinercv1	20,242	47,326	2
\hdashlinegisette	6,000	5,000	2
\hdashlinew8a	49,749	300	2
\hdashline

Equations

$$\underset{w\in {\mathbb{R}}^d}{\min }F(w):=\frac{1}{n}\sum _{i=1}^{n}f(w;x^i,y^i)=\frac{1}{n}\sum _{i=1}^n{f}_i(w)$$

$$w_{k+1}=w_k-{\alpha }_k{H}_k\nabla F(w_k),$$

$$H_{k+1}=V_k^T{H}_k{V}_k+{\rho }_k{s}_k{s}_k^T,$$

$${\rho }_k=\frac{1}{y_k^T{s}_k},V_k=I-{\rho }_k{y}_k{s}_k^T,$$

$$s_k=w_k-w_{k-1},y_k=\nabla F(w_k)-\nabla F(w_{k-1}).$$

$$H_{k+1}{y}_k=s_k,\text{and}{s}_k^T{y}_k>0,$$

$$w_{k+1}=w_k+p_k,$$

$$\underset{p}{\min }{m}_k(p)=F(w_k)+\nabla F(w_k{)}^{T}p+\frac{1}{2}{p}^T{B}_{k}p,$$

$$\text{s.t.}\Vert p\Vert \le {\Delta }_k,$$

$$B_{k+1}=B_k+\frac{(y_k-B_k{s}_k)(y_k-B_k{s}_k{)}^T}{(y_k-B_k{s}_k{)}^T{s}_k}.$$

$$F(w_k+{\alpha }_k{p}_k)\le F(w_k)-c_1{\alpha }_k\nabla F(w_k{)}^T{H}_k\nabla F(w_k)$$

$$\text{Cost}=\begin{array}{c}\text{Cost of gradient}\\ \text{computation}\end{array}+\begin{array}{c}\text{Cost of forming/}\\ \text{taking step.}\end{array}$$

$$\mu I⪯{\nabla }^{2}F(w)⪯LI,\text{for all }w\in {\mathbb{R}}^d.$$

$${\mu }_{1}I⪯H_k⪯{\mu }_{2}I,\text{for }k=0,1,2,\dots .$$

$$B_k^{(i)}=B_k^{(i-1)}-\frac{B_k^{(i-1)}{s}_{k,i}{s}_{k,i}^T{B}_k^{(i-1)}}{s_{k,i}^T{B}_k^{(i-1)}{s}_{k,i}}+\frac{y_{k,i}{y}_{k,i}^T}{y_{k,i}^T{s}_{k,i}}.$$

$$s_{k,j}=w_k-{\bar{w}}_j,y_{k,j}=\nabla F(w_k)-\nabla F({\bar{w}}_j)\text{Option I},$$

$$s_{k,j}=w_k-{\bar{w}}_j,y_{k,j}={\nabla }^{2}F(w_k)s_k\text{Option II}.$$

$$y_{k,j}^T{s}_{k,j}\ge \frac{1}{L}\Vert y_{k,j}{\Vert }^2$$

$$y_{k,j}^T{s}_{k,j}\le \frac{1}{\mu }\Vert y_{k,j}{\Vert }^2$$

$$\mu \le \frac{\Vert y_{k,j}{\Vert }^2}{y_{k,j}^T{s}_{k,j}}\le L.$$

$$\mu \Vert s_{k,j}{\Vert }^2\le y_{k,j}^T{s}_{k,j}=s_{k,j}^T{\nabla }^{2}F(w_k)s_{k,j}\le L\Vert s_{k,j}{\Vert }^2.$$

$$\frac{\Vert y_{k,j}{\Vert }^2}{y_{k,j}^T{s}_{k,j}}=\frac{s_{k,j}^T{\nabla }^{2}F(w_k{)}^2{s}_{k,j}}{s_{k,j}^T{\nabla }^{2}F(w_k)s_{k,j}},$$

$$\mu \le \frac{\Vert y_{k,j}{\Vert }^2}{y_{k,j}^T{s}_{k,j}}\le L.$$

$$Tr(B_{k+1})$$

$$\le Tr(B_k^{(0)})+\sum _{i=1}^m\frac{\Vert y_{k,i}{\Vert }^2}{y_{k,i}^T{s}_{k,i}}$$

$$\le Tr(B_k^{(0)})+mL\le C_1,$$

$$\det (B_{k+1})$$

$$=\det (B_k^{(0)})\prod _{i=1}^m\frac{y_{k,i}^T{s}_{k,i}}{s_{k,i}^T{s}_{k,i}}\frac{s_{k,i}^T{s}_{k,i}}{s_{k,i}^T{B}_k^{(i-1)}{s}_{k,i}}$$

$$\displaystyle\geq\det(B_{k}^{(0)})\Big{(}\frac{\mu}{C_{1}}\Big{)}^{m}\geq C_{2},$$

$$F(w_k)-F^{\star }$$

$$F(w_{k+1})$$

$$\le F(w_k)+\nabla F(w_k{)}^T(-\alpha H_k\nabla F(w_k))+\frac{L}{2}\Vert \alpha H_k\nabla F(w_k){\Vert }^2$$

$$\le F(w_k)-\alpha \nabla F(w_k{)}^T{H}_k\nabla F(w_k)+\frac{{\alpha }^2{\mu }_2^{2}L}{2}\Vert \nabla F(w_k){\Vert }^2$$

$$\le F(w_k)-\alpha {\mu }_1\Vert \nabla F(w_k){\Vert }^2+\frac{{\alpha }^2{\mu }_2^{2}L}{2}\Vert \nabla F(w_k){\Vert }^2$$

$$=F(w_k)-\alpha \left({\mu }_1-\alpha \frac{{\mu }_2^{2}L}{2}\right)\Vert \nabla F(w_k){\Vert }^2$$

$$\le F(w_k)-\alpha \frac{{\mu }_1}{2}\Vert \nabla F(w_k){\Vert }^2,$$

Full text

1 Introduction

In supervised machine learning, one seeks to minimize the empirical risk,

[TABLE]

where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is the composition of a prediction function (parametrized by $w$) and a loss function, and $(x^{i},y^{i})$, for $i=1,\dots,n$, denote the training examples (samples). Difficulties arise in minimizing the function $F$ for three main reasons: $(1)$ the number of samples $n$ is large; $(2)$ the number of variables $d$ is large; and, $(3)$ the objective function is nonconvex.

In the last decades, much effort has been devoted to the development of stochastic first-order methods that have a low per-iteration cost, enjoy optimal complexity, are easy to implement, and that have proven to be effective for many machine learning applications. At present, the preferred method for large-scale applications is the stochastic gradient (SG) method [55, 7], and its variance-reduced [33, 57, 21, 50] and adaptive variants [22, 37]. However, these methods have several issues: (1) they are highly sensitive to the choice of hyper-parameters (e.g., steplength and batch size) and tuning can be cumbersome; (2) they suffer from ill-conditioning; and, (3) they often offer limited opportunities for parallelism; see [2, 62, 39, 56, 8].

In order to alleviate these issues, stochastic Newton [11, 46, 5, 56, 61] and stochastic quasi-Newton [58, 12, 18, 49, 26, 3, 34, 4, 30, 31] methods have been proposed. These methods attempt to combine the speed of Newton’s method and the scalability of first-order methods by incorporating curvature information in a judicious manner, and have proven to work well for several machine learning tasks [2, 62].

With the advances in distributed and GPU computing, it is now possible to go beyond stochastic Newton and quasi-Newton methods and use large batches, or even the full dataset, to compute function, gradient and Hessian vector products in order to train machine learning models. In the large batch regime, one can take advantage of parallel and distributed computing and fully utilize the capabilities of GPUs. However, researchers have observed that well-tuned first-order methods (e.g., ADAM) are far more effective than full batch methods (e.g., LBFGS) for large-scale applications [28, 35].

Nevertheless, in this paper we focus on (full) batch methods that incorporate local second-order (curvature) information of the objective function. These methods mitigate the effects of ill-conditioning, avoid or diminish the need for hyper-parameter tuning, have enough concurrency to take advantage of parallel computing, and, due to requiring fewer iterations enjoy low communication costs in distributed computing environments. Specifically, we focus on quasi-Newton methods [52]; methods that construct curvature information using first-order (gradient) information. We propose two variants of classical quasi-Newton methods that sample a small number of random points at every iteration to build (inverse) Hessian approximations.

We are motivated by the results presented in Figure 1 that illustrate the performance (for 10 different starting points) of several stochastic and deterministic, first- and second-order methods on a toy neural network classification task, given budget; see Section 7 for details. As is clear from the results, first-order methods converge very slowly, and sometimes even fail to achieve 100% accuracy. Similarly, classical quasi-Newton methods are also slow or stagnate. On the other hand, methods that use the true Hessian are able to converge in very few iterations from all starting points. This seems to suggest that for some neural network training tasks second-order information is important, and that the curvature information captured by classical quasi-Newton methods may not be adequate or useful.

The key idea of our proposed methods is to leverage the fact that quasi-Newton methods can incorporate second-order information using only gradient information at a reasonable cost, but at the same time to enhance the (inverse) Hessian approximations by using more reliable (recent and local) information. The fundamental component of our methods, and what differentiates them from the classical variants, is the manner in which the curvature pairs are constructed. To this end, we propose to forget past curvature information and sample new curvature pairs at every iteration.

Contributions

Our contributions can be summarized as follows:

• •

We propose two novel quasi-Newton methods that use sampling to construct Hessian approximations. We analyze the convergence properties of both methods, and show that their theoretical guarantees match those of their classical limited memory counterparts.

• •

We derive expressions for the probability of accepting curvature pairs within our quasi-Newton frameworks. Moreover, we illustrate empirically the probability of accepting curvature pairs on a toy classification problem.

• •

We discuss the implementation costs of the sampled quasi-Newton methods and compare them to the classical variants, and illustrate the scaling properties of the methods compared to the SG method on distributed computing platforms on real large-scale network architectures.

• •

We illustrate the practical performance of the methods on a toy classification problem as well as on popular benchmarking binary classification and neural network training tasks, and show their advantages over the classical variants. We posit that this is the case since the (inverse) Hessian approximations constructed by our proposed methods capture better (more informative) curvature information. Moreover, the proposed methods are easily parallelizable and efficient in terms of iteration, epochs and communication.

The paper is organized as follows. We conclude this section with a literature review of quasi-Newton methods. We describe the classical (L)BFGS and (L)SR1 methods in Section 2, and in Section 3 we present our proposed sampled quasi-Newton variants. In Section 4, we discuss the computational cost of the proposed methods and show their scaling properties. We show the theoretical properties of our proposed methods in Section 5. In Section 6 we provide a theoretical and numerical analysis on the probabilities of accepting the sampled points within our proposed quasi-Newton frameworks. Numerical results on neural network training tasks are reported in Section 7. Finally, in Section 8 we provide some final remarks and discuss several avenues for future work.

Literature Review

Quasi-Newton methods, such as BFGS [9, 23, 25, 59] and SR1 [13, 16, 36] and their limited-memory variants LBFGS [51, 43] and LSR1 [10, 45], respectively, have been studied extensively in the deterministic nonlinear optimization literature. These methods incorporate curvature (second-order) information using only gradient (first-order) information, have good theoretical guarantees, and have proven to be effective in practice.

In the context of deep neural networks, both full batch and stochastic quasi-Newton methods seem to perform worse than (stochastic) first-order methods. Nevertheless, several stochastic quasi-Newton methods have been proposed; see e.g., [58, 12, 4]. What distinguishes these methods from one another is the way in which curvature pairs are constructed. Our methods borrow some of the ideas proposed in [12, 26, 44]. Specifically, we use Hessian vector products in lieu of gradient displacements.

Possibly the closest works to ours are Block BFGS [24] and its stochastic variant [26]. These methods construct multiple curvature pairs to update the quasi-Newton matrices. However, there are several key features that are different from our approach; in these works $(1)$ the Hessian approximation is not updated at every iteration, and $(2)$ they enforce that multiple secant equations hold simultaneously.

2 Quasi-Newton Methods

In this section, we review two classical quasi-Newton methods (BFGS and SR1) and their limited memory variants (LBFGS and LSR1). This will set the stage for our proposed sampled quasi-Newton methods.

2.1 BFGS and LBFGS

Let us begin by considering the BFGS method and then consider its limited memory version. At the $k$th iteration, the BFGS method computes a new iterate by the formula

[TABLE]

where $\alpha_{k}$ is the step length, $\nabla F(w_{k})$ is the gradient of (1.1) and $H_{k}$ is the inverse BFGS Hessian approximation that is updated at every iteration by means of the formula

[TABLE]

where the curvature pairs $(s_{k},y_{k})$ are defined as

[TABLE]

As is clear, the curvature pairs (2.2) are constructed sequentially (at every iteration), and as such the inverse Hessian approximation at the $k$th iteration $H_{k}$ depends on iterate (and gradient) information from past iterations.

The inverse BFGS Hessian approximations are constructed to satisfy two conditions:

[TABLE]

the secant and curvature conditions, respectively, as well as symmetry. Consequently, as a result, as long as the initial inverse Hessian approximation is positive definite, then all subsequent inverse BFGS Hessian approximations are also positive definite. Note, the new (inverse) Hessian approximation $H_{k+1}$ differs from the old approximation $H_{k}$ by a rank-2 matrix.

In the limited memory version, the matrix $H_{k}$ is defined at each iteration as the result of applying $m$ BFGS updates to a multiple of the identity matrix using the set of $m$ most recent curvature pairs $\{s_{i},y_{i}\}$ kept in storage. As a result, one need not store the dense inverse Hessian approximation, rather one can store two $m\times d$ matrices and compute the matrix-vector product in (2.1) via the two-loop recursion [52]. After the step has been computed, the oldest pair $(s_{j},y_{j})$ is discarded and the new curvature pair is stored.

2.2 SR1 and LSR1

Contrary to the BFGS updating formula, and as suggested by the name, the symmetric-rank-1 (SR1) updating formula allows one to satisfy the secant equation and maintain symmetry with a simpler rank-1 update. However, unlike BFGS, the SR1 update does not guarantee that the updated matrix maintains positive definiteness. As such, the SR1 method is usually implemented with a trust region; we introduce it in this way below.

At the $k$th iteration, the SR1 method computes a new iterate by the formula

[TABLE]

where $p_{k}$ is the minimizer of the following subproblem

[TABLE]

$\Delta_{k}$ is the trust region and $B_{k}$ is the SR1 Hessian approximation computed as

[TABLE]

Similar to LBFGS, in the limited memory version of SR1 the matrix $B_{k}$ is defined at each iteration as the result of applying $m$ SR1 updates to a multiple of the identity matrix, using a set of $m$ correction pairs $\{s_{i},y_{i}\}$ kept in storage.

3 Sampled Quasi-Newton Methods

In this section, we describe our two proposed sampled quasi-Newton methods; S-LBFGS and S-LSR1. The main idea of these methods, and what differentiates them from the classical variants, is the way in which curvature pairs are constructed. At every iteration, a small number ($m$) of points are sampled around the current iterate and used to construct a new set of curvature pairs. In other words, contrary to the sequential nature of classical quasi-Newton methods, our proposed methods forget all past curvature pairs and construct new curvature pairs from scratch via sampling.

Our motivation stems from the following observation: by constructing Hessian approximations via sampling, one is able to better capture curvature information of the objective function. In Figures 2 and 3, we show the spectrum of the true Hessian, and compare it to the spectra of different SR1 Hessian approximations at several points for two toy classification problems. As is clear from the results, the eigenvalues of the S-LSR1 Hessian approximations better match the eigenvalues of the true Hessian compared to the eigenvalues of the SR1 and LSR1 Hessian approximations. This is not surprising since S-LSR1 uses newly sampled local information, and unlike the classical variants does not rely on past information that could be significantly stale. Similar results were obtained for other problems; see Appendix A.2 for details.

This, of course, does not come for free. The classical variants construct curvature pairs as the optimization progresses at no additional cost, whereas the sampled quasi-Newton methods require the construction of $m$ new curvature pairs at every iteration. We discuss implementation issues and the computational cost of the sampled quasi-Newton methods in Sections 3.1 and 4.

We now discuss the way in which curvature pairs are constructed, and then formally introduce our proposed sampled quasi-Newton methods.

3.1 Sampling Curvature Pairs

As mentioned above, the key component of our proposed algorithms is the way in which curvature pairs are constructed. A pseudo-code of our proposed sampling strategy and construction of the curvature pairs in given in Algorithm 1. Let $S\in\mathbb{R}^{d\times m}$ and $Y\in\mathbb{R}^{d\times m}$ denote the matrices of all curvature pairs constructed during the $k$th iteration.

Both S-LBFGS and S-LSR1 use the subroutine described in Algorithm 1. At every iteration, given the current iterate and gradient, $m$ curvature pairs are constructed. The subroutine first samples points around the current iterate along a random directions $\sigma_{i}$ and sets the iterate displacement curvature pair ($s$), and then creates the gradient difference curvature pair ($y$) via gradient differences (Option I) or Hessian vector products (Option II). Note that the random directions $\sigma_{i}$ can be arbitrary; in the latter part of the paper (Sections 7 and 6), we make an explicit choise on the directions.

Our theory holds for both options; however, in our numerical experiments we present results with Option II only for the following reasons. Option I requires $m$ gradient evaluations ($m$ epochs), and thus requires accessing the data $m$ times. On the other hand, Option II only requires a single Hessian matrix product which can be computed very efficiently on a GPU, as the $y$ curvature pairs can be constructed simultaneously, i.e., $Y=\nabla^{2}F(w)S$, and thus only requires accessing the data once. Moreover, Option I requires choosing the sampling radius $r$, whereas Option II does not since it is scale invariant.

Before we proceed with our presentation of the S-LBFGS and S-LSR1 methods, we empirically compare the performance of a methods that uses Option I and Option II. As is clear from Figures 4 and 5, the performance of the method that uses Option I is highly dependant on the choice of the sampling radius ($r$). If this parameter is not chosen appropriately, the performance of the method can be slow. This is not the case when Option II is utilized, and one can attribute this to the fact that Option II is scale invariant. Moreover, the benefits of using Option II can clearly be observed in the plots in terms of epochs. Again, this is due to the fact each iteration using Option I requires accessing the data at $m$ times to construct the curvature pairs, whereas Option II required only a single access of the data to construct the curvature pairs.

3.2 Sampled LBFGS (S-LBFGS)

At the $k$th iteration, the S-LBFGS method computes a new iterate via (2.1), where the inverse Hessian approximation is constructed using the curvature pairs sampled by Algorithm 1. The S-LBFGS method is outlined in Algorithm 2.

Algorithm 2 is almost identical to the classical (L)BFGS algorithm [52]; however, it has two key differentiating elements: (1) the way in which curvature pairs are created; and, (2) the location in the algorithm where the curvature pairs are constructed. Both elements can be interpreted as features of S-LBFGS. First, using a similar argument as that for the S-LSR1 method (Figure 2), the inverse Hessian approximations constructed by this method better capture local curvature information of the objective function. Moreover, notice that the first set of curvature pairs is constructed before a single step is taken by the method (Line 2). This allows the method to take a quasi-Newton-type (well-scaled) step from the first iteration which is not the case for classical BFGS methods that usually take a gradient-type step in the first iteration, and in which imposing the correct scale can be an issue. This, possibly, is a more important implication of the method, as the first step taken by quasi-Newton methods can be of paramount importance.

In order to fully specify the S-LBFGS method, we need to describe how the steplength is selected (Algorithm 2, Step 4). We consider two variants of the method: $(1)$ constant steplength selection, and $(2)$ adaptive steplength selection. Our theory (Section 5, Theorems 5.2 and 5.5), explicitly defines the manner in which the steplength should be selected in order to ensure convergence. Of course, in practice, one can (potentially) use a larger steplength, and as such in this approach the steplength ($\alpha_{k}=\alpha$) is a tuneable parameter. We also consider an adaptive Armijo backtracking mechanism for selecting the steplength at every iteration. Given the current iterate $w_{k}$, the steplength is chosen to satisfy the following sufficient decrease condition

[TABLE]

where $c_{1}\in(0,1)$. The mechanism works as follows. Given an initial steplength (say $\alpha_{k}=1$), the function is evaluated at the trial point $w_{k}+\alpha_{k}p_{k}$ and condition (3.1) is checked. If the trial point satisfies (3.1), then the step is accepted. If the trial point does not satisfy (3.1), the steplength is reduced (e.g., $\alpha_{k}=\tau\alpha_{k}$ for $\tau\in(0,1)$). This process is repeated until a steplength that satisfies (3.1) is found. We should note that under reasonable assumptions on the function $F$ (see [52]) this procedure is well defined since the search direction uses the true gradient, $H_{k}$ is a positive definite matrix, and the true function is used in condition (3.1).

3.3 Sampled LSR1

At the $k$th iteration, the S-LSR1 method computes a new iterate via (2.3), where the Hessian approximation in (2.4) is constructed using the curvature pairs sampled by Algorithm 1. The S-LSR1 method is outlined in Algorithm 3.

The S-LSR1 method has the same key features as S-LBFGS that differentiates it from the classical SR1 methods. The subroutine adjustTR (Step 10, Algorithm 3) adjusts the trust-region based on the progress made by the method. For brevity we omit the details of this subroutine, and refer the reader to Appendix A.4 for the details.

4 Distributed Computing and Computational Cost

In this section, we show the scalability of the sampled quasi-Newton methods as compared to the SG method, and compare the computational cost to the classical variants.

4.1 Distributed Computing

Recently, there has been a huge effort to scale SG-type algorithms to solve Imagenet using hundreds of GPUs; see e.g., [27, 32, 1, 63]. In Figure 6 (left), we show how the batch size affects the number of images processed per second to compute the function, gradient and Hessian vector products on a NVIDIA Tesla P100 GPU for various deep neural networks111The structure of the deep neural network is taken from: https://github.com/tensorflow/models/tree/master/research/slim.; see Table 1.

As is clear, by using small batch sizes one is not able to fully utilize the power of GPUs. On the other hand, using larger batches in conjunction with SG-type algorithms does not necessarily reduce training time [20, 60]. Another observation that can be extracted from Figure 6 is that the cost of computing function values, gradients and Hessian vector products appears to be comparable for these networks.

In Figure 6 (bar plots), we compare the time to perform 1 epoch of the SG method (assuming we have 1M images) with the time to perform 1 iteration of S-LSR1. For SG, we show results for different batch sizes on each GPU222 Each GPU has 1 MPI process that is used for communicating updates. Note, we are running 4 MPI processes for each physical node, i.e., each node has 4 P100 GPUs: (1) batch size 16 (SGD 16); and, (2) batch size 32, 64 and 128 for vgg a, LeNet and alexnet v2, respectively (SGD Default). The reason that there is no significant benefit when using more GPUs for the SG method is that the cost is dominated by the communication. For S-LSR1, that is not the case; as we scale up the number of MPI processes, we get good performance improvements since there is much less communication involved. See Appendix A.5 for more details.

4.2 Cost, Storage and Parallelization

The cost per iteration of the different quasi-Newton methods can be deconstructed as follows:

[TABLE]

Note, motivated by the results in Figure 6, we assume that the cost computing a function value, gradient and Hessian vector product is comparable and is $\mathcal{O}(nd)$. The cost of computing the gradient is common for each method, whereas the search directions are computed differently for BFGS-type methods and SR1-type methods. More specifically, for BFGS methods we employ a line search and for SR1 method we use a trust region and solve the subproblem (2.4) using CG [52]. We denote the number of line search iterations and CG iterations as $\kappa_{ls}$ and $\kappa_{tr}$, respectively. Table 2 summarizes the computational cost and storage for the different quasi-Newton methods.

As is clear from Table 2, the proposed sampled quasi-Newton methods do not have a significantly higher cost per iteration than the classical limited memory variants of the methods. In the regime where $m\ll n,d$, the computational cost of the methods are $\mathcal{O}(nd)$. Moreover, the storage requirements for the sampled quasi-Newton methods is the same as that of limited-memory quasi-Newton methods. We should also note, that several computations that are required in our proposed methods are easily parallelizeable. These computations are the gradient evaluations, the function evaluations and the construction of the gradient displacement curvature pairs $y$.

5 Convergence Analysis

In this section, we present convergence analyses for the sampled quasi-Newton methods.

5.1 Sampled LBFGS

We derive convergence results for the sampled LBFGS method with fixed step sizes and adaptive step sizes for strongly convex and nonconvex functions.

5.1.1 Strongly Convex Functions

We make the following standard assumptions.

Assumption 1**.**

$F$* is twice continuously differentiable.*

Assumption 2**.**

There exist positive constants $\mu$ and $L$ such that

[TABLE]

First, we show that the inverse Hessian approximations $H_{k}$ generated by the sampled LBFGS method have eigenvalues that are uniformly bounded above and away from zero. The proof technique is an adaptation of that in [12, 3]; however, modifications are necessary since in our approach the inverse Hessian approximations are constructed using information only from the current iterate, and not constructed sequentially.

Lemma 5.1**.**

If Assumptions 1 and 2 hold, there exist constants $0<\mu_{1}\leq\mu_{2}$ such that the inverse Hessian approximations $\{H_{k}\}$ generated by Algorithm 2 satisfy,

[TABLE]

Proof.

Instead of analyzing the inverse Hessian approximation $H_{k}$, we study the direct Hessian approximation $B_{k}=H_{k}^{-1}$. In this case, the sampled LBFGS updating formula is given as follows. At the $k$th iteration, given a set of curvature pairs $(s_{k,j},y_{k,j})$, for $j=1,\dots,m$

Set $B_{k}^{(0)}=\frac{y_{k,l}^{T}y_{k,l}}{s_{k,l}^{T}y_{k,l}}I$, where $l$ is chosen uniformly at random from $\{1,\dots,m\}$. 2. 2.

For $i=1,\dots,m$ compute

[TABLE] 3. 3.

Set $B_{k+1}=B_{k}^{(m)}$.

In our algorithm (Algorithm 1), there are two options for constructing the curvature pairs $s_{k,j}$ and $y_{k,j}$. At the current iterate $w_{k}$ we sample points $\bar{w}_{j}$ for $j=1,\dots,m$ and set

[TABLE]

We now derive an upper and lower bound for $\frac{\|y_{k,j}\|^{2}}{y_{k,j}^{T}s_{k,j}}$, for all $j=1,\dots,m$, for both options.

Option I: A consequence of Assumption 2 is that the eigenvalues of the Hessian matrix are bounded above and away from zero. Utilizing this fact, the convexity of the objective function and the definitions (5.2), we have

[TABLE]

On the other hand, strong convexity of the functions, the consequence of Assumption 2 and definitions (5.2), provide a lower bound,

[TABLE]

Combining the upper and lower bounds (5.4) and (5.5)

[TABLE]

Option II: A consequence of Assumption 2 is that the eigenvalues of the Hessian matrix are bounded above and away from zero. Utilizing this fact and the definitions (5.2), we have

[TABLE]

We have that,

[TABLE]

and since $\nabla^{2}F(w_{k})$ is symmetric and positive definite, it has a square root and so

[TABLE]

The bounds on $\frac{\|y_{k,j}\|^{2}}{y_{k,j}^{T}s_{k,j}}$ prove that for any $l$ chosen uniformly at random from $\{1,\dots,m\}$ the eigenvalues of the matrices $B_{k}^{(0)}=\frac{y_{k,l}^{T}y_{k,l}}{s_{k,l}^{T}y_{k,l}}I$ at the start of the sampled LBFGS update cycles are bounded above and away from zero, for all $k$ and $l$. We now use a Trace-Determinant argument to show that the eigenvalues of $B_{k}$ are bounded above and away from zero.

Let $Tr(B)$ and $\det(B)$ denote the trace and determinant of matrix $B$, respectively. The trace of the matrix $B_{k+1}$ can be expressed as,

[TABLE]

for some positive constant $C_{1}$, where the inequalities above are due to (5.6), and the fact that the eigenvalues of the initial L-BFGS matrix $B_{k}^{(0)}$ are bounded above and away from zero.

Using a result due to Powell [53], the determinant of the matrix $B_{k+1}$ generated by the sampled LBFGS method can be expressed as,

[TABLE]

for some positive constant $C_{2}$, where the above inequalities are due to the fact that the largest eigenvalue of $B_{k}^{(i)}$ is less than $C_{1}$, Assumption 2, and the fact that $\frac{\mu}{C_{1}}<1$.

The trace (5.1.1) and determinant (5.1.1) inequalities derived above imply that largest eigenvalues of all matrices $B_{k}$ are bounded above, uniformly, and that the smallest eigenvalues of all matrices $B_{k}$ are bounded away from zero, uniformly. ∎

Constant Step Length

Utilizing Lemma 5.1, we show that the sampled LBFGS method with a constant step length converges linearly.

Theorem 5.2**.**

Suppose that Assumptions 1 and 2 hold, and let $F^{\star}=F(w^{\star})$, where $w^{\star}$ is the minimizer of $F$. Let $\{w_{k}\}$ be the iterates generated by Algorithm 2, where $0<\alpha_{k}=\alpha\leq\frac{\mu_{1}}{\mu_{2}^{2}L}$, and $w_{0}$ is the starting point. Then, for all $k\geq 0$,

[TABLE]

Proof.

We have that

[TABLE]

where the first inequality is due to Assumption 2, the second and third inequalities arise as a consequence of Lemma 5.1 and the last inequality is due to the choice of the steplength. By strong convexity, we have $2\mu(F(w)-F^{\star})\leq\|\nabla F(w)\|^{2}$, and thus

[TABLE]

Subtracting $F^{\star}$ from both sides,

[TABLE]

Recursive application of the above inequality yields the desired result. ∎

Theorem 5.2 shows that the S-LBFGS method converges to the optimal solution at a linear rate. This result is similar in nature to the result for LBFGS [43]. We should also mention that unlike first-order methods (e.g., Gradient Descent, $H_{k}=I$), the step length range and the rate of convergence of the S-LBFGS method depends on $\mu_{1}$ and $\mu_{2}$, the smallest and largest eigenvalues of the S-LBFGS Hessian approximation. In the worst-case, the presence of the matrix $H_{k}$ can make the results in Theorem 5.2 significantly worse than that of the first-order variant if the update has been unfortunate and generates ill-conditioned matrices. We should note, however, such worst-case behavior is almost never observed in practice for BFGS updating.

Adaptive Step Length

We now show a similar result for the case where the step length is chosen by an Armijo backtracking line search (3.1).

Theorem 5.3**.**

Suppose that Assumptions 1 and 2 hold. Let $\{w_{k}\}$ be the iterates generated by Algorithm 2, where $\alpha_{k}$ is the maximum value in $\{\tau^{-j}:j=0,1,\dots\}$ satisfying (3.1) with $0<c_{1}<1$, and $w_{0}$ is the starting point. Then for all $k\geq 0$,

[TABLE]

Proof.

Starting with (5.12) we have

[TABLE]

From the Armijo backtracking condition (3.1), we have

[TABLE]

Looking at (5.12) and (5.14), it is clear that the Armijo condition is satisfied when

[TABLE]

Thus, any $\alpha_{k}$ that satisfies (5.15) is guaranteed to satisfy the Armijo condition (3.1). Since we find $\alpha_{k}$ using a constant backtracking factor of $\tau<1$, we have that

[TABLE]

Therefore, from (5.12) and by (5.15) and (5.16) we have

[TABLE]

By strong convexity, we have $2\mu(F(w)-F^{\star})\leq\|\nabla F(w)\|^{2}$, and thus

[TABLE]

Subtracting $F^{\star}$ from both sides, and applying (5.18) recursively yields the desired result. ∎

Theorem 5.3 shows that the sampled LBFGS method with an adaptive backtracking line search converges to the optimal solution at a linear rate. We should note that this result is worse than the constant step length result (Theorem 5.2), i.e., the rate constant is larger. This is not surprising since this is a worst-case result; however, in practice, such an approach performs significantly better and does not require tuning the steplength parameter.

5.1.2 Nonconvex Functions

For nonconvex functions, the BFGS method is known fail [47, 19]. Even for LBFGS, which makes only a finite number of updates at each iteration, one cannot guarantee that the (inverse) Hessian approximations have eigenvalues that are uniformly bounded above and away from zero. To establish convergence of the BFGS method in the nonconvex setting several techniques have been proposed including: $(i)$ cautious updating [42]; $(ii)$ modified updating [41]; and $(iii)$ damping [54]. Here we employ a cautious strategy that is well suited to our particular algorithm; at the $k$th iteration, we update the (inverse) Hessian approximation using only the set of curvature pairs that satisfy

[TABLE]

where $\epsilon>0$ is a predetermined constant. If no curvature pairs satisfy (5.19), then the new (inverse) Hessian approximation is set to $H_{k}=I$. Using said mechanism we prove that the eigenvalues of the (inverse) Hessian approximations generated by the S-LBFGS method are bounded above and away from zero. For this analysis, we make the following assumptions in addition to Assumption 1.

Assumption 3**.**

The function $F(w)$ is bounded below by a scalar $\widehat{F}$.

Assumption 4**.**

The gradients of $F$ are $L$-Lipschitz continuous for all $w\in\mathbb{R}^{d}$.

Lemma 5.4**.**

Suppose that Assumptions 1 and 4 hold. Let $\{H_{k}\}$ be the inverse Hessian approximations generated by Algorithm 2, with the modification that the inverse approximation update is performed using only curvature pairs that satisfy (5.19), for some $\epsilon>0$, and $H_{k}=I$ if no curvature pairs satisfy (5.19). Then, there exist constants $0<\mu_{1}\leq\mu_{2}$ such that

[TABLE]

Proof.

Note, that in the nonconvex setting, there is a chance that no curvature pairs are selected in Algorithm 1. In this case, the inverse Hessian approximation is $H_{k}=I$, and thus $\mu_{1}=\mu_{2}=1$ and condition (5.20) is satisfied.

Similar to the proof of Lemma 5.1, we study the direct Hessian approximation $B_{k}=H_{k}^{-1}$. In our algorithm, there are two options for updating the curvature pairs $s_{k,j}$ and $y_{k,j}$:

[TABLE]

for $j=1,\dots,m$. Let $\tilde{m}_{k}\in\{1,...,m\}$ denote the number of curvature pairs that satisfy (5.19) at the $k$th iteration, where $m$ is the memory. At the $k$th iteration, given a set of curvature pairs $(s_{k,j},y_{k,j})$, for $j=1,\dots,\tilde{m}_{k}$ we update the Hessian approximation recursively (using the procedure described in the proof of Lemma 5.1, and set $B_{k+1}=B_{k}^{\tilde{m}_{k}}$.

In this setting, the skipping mechanism (5.19) provides both an upper and lower bound on the quantity $\frac{\|y_{k,j}\|^{2}}{y_{k,j}^{T}s_{k,j}}$, for both Options, which in turn ensures that the initial sampled LBFGS Hessian approximation is bounded above and away from zero.

The lower bound is attained by repeated application of Cauchy’s inequality to condition (5.19). We have from (5.19) that

[TABLE]

It follows that

[TABLE]

The upper bound is attained by the Lipschitz continuity of gradients,

[TABLE]

Combining (5.23) and (5.24), we have

[TABLE]

The bounds on $\frac{\|y_{k,j}\|^{2}}{y_{k,j}^{T}s_{k,j}}$ prove that for any $l$ chosen uniformly at random from $\{1,\dots,\tilde{m}_{k}\}$ the eigenvalues of the matrices $B_{k}^{(0)}=\frac{y_{k,l}^{T}y_{k,l}}{s_{k,l}^{T}y_{k,l}}I$ at the start of the sampled LBFGS update cycles are bounded above and away from zero, for all $k$ and $l$. The rest of the proof follows the same trace-determinant argument as in the proof of Lemma 5.1, the only difference being that the last inequality in 5.1.1 comes as a result of the cautious update strategy. ∎

Constant Step Length

Utilizing Lemma 5.4, we show that the sampled LBFGS with a cautious updating strategy and a constant step length converges.

Theorem 5.5**.**

Suppose that Assumptions 1, 3 and 4 hold. Let $\{w_{k}\}$ be the iterates generated by Algorithm 2, with the modification that the inverse Hessian approximation update is performed using only curvature pairs that satisfy (5.19), for some $\epsilon>0$, and $H_{k}=I$ if no curvature pairs satisfy (5.19), where $0<\alpha_{k}=\alpha\leq\frac{\mu_{1}}{\mu_{2}^{2}L}$, and $w_{0}$ is the starting point. Then,

[TABLE]

and, moreover, for any $T>1$,

[TABLE]

Proof.

We start with (5.13)

[TABLE]

Summing both sides of the above inequality from $k=0$ to $T-1$,

[TABLE]

The left-hand-side of the above inequality is a telescopic sum and thus,

[TABLE]

where the inequality is due to $\hat{F}\leq F(w_{T})$ (Assumption 3). Using the above, we have

[TABLE]

Taking limits we obtain,

[TABLE]

which implies (5.25). Dividing (5.26) by $T$ we conclude

[TABLE]

∎

Theorem 5.5 shows that, if a small enough step length is chosen, the S-LBFGS method converges to a stationary point.

Adaptive Step Length

We now show a similar result for the case where the step length is chosen by an Armijo backtracking line search (3.1).

Theorem 5.6**.**

Suppose that Assumptions 1, 3 and 4 hold. Let $\{w_{k}\}$ be the iterates generated by Algorithm 2, with the modification that the inverse Hessian approximation update is performed using only curvature pairs that satisfy (5.19), for some $\epsilon>0$, and $H_{k}=I$ if no curvature pairs satisfy (5.19), where $\alpha_{k}$ is the maximum value in $\{\tau^{-j}:j=0,1,\dots\}$ satisfying (3.1) with $0<c_{1}<1$, and where $w_{0}$ is the starting point. Then,

[TABLE]

and, moreover, for any $T>1$,

[TABLE]

Proof.

We start with (5.1.1)

[TABLE]

Summing both sides of the above inequality from $k=0$ to $T-1$,

[TABLE]

The left-hand-side of the above inequality is a telescopic sum and thus,

[TABLE]

where the inequality is due to $\hat{F}\leq F(w_{T})$ (Assumption 3). Using the above, we have

[TABLE]

Taking limits we obtain,

[TABLE]

which implies (5.29). Dividing (5.30) by $T$ we conclude

[TABLE]

∎

Theorem 5.7 shows that, the S-LBFGS method that employs an Armijo backtracking linesearch (3.1) converges to a stationary point.

Adaptive Step Length

We now show a similar result for the case where the step length is chosen by an Armijo backtracking line search (3.1).

Theorem 5.7**.**

Suppose that Assumptions 1, 3 and 4 hold. Let $\{w_{k}\}$ be the iterates generated by Algorithm 2, with the modification that the inverse Hessian approximation update is performed using only curvature pairs that satisfy (5.19), for some $\epsilon>0$, and $H_{k}=I$ if no curvature pairs satisfy (5.19), where $\alpha_{k}$ is the maximum value in $\{\tau^{-j}:j=0,1,\dots\}$ satisfying (3.1) with $0<c_{1}<1$, and where $w_{0}$ is the starting point. Then,

[TABLE]

and, moreover, for any $T>1$,

[TABLE]

Proof.

We start with (5.1.1)

[TABLE]

Summing both sides of the above inequality from $k=0$ to $T-1$,

[TABLE]

The left-hand-side of the above inequality is a telescopic sum and thus,

[TABLE]

where the inequality is due to $\hat{F}\leq F(w_{T})$ (Assumption 3). Using the above, we have

[TABLE]

Taking limits we obtain,

[TABLE]

which implies (5.29). Dividing (5.30) by $T$ we conclude

[TABLE]

∎

Theorem 5.7 shows that, the S-LBFGS method that employs an Armijo backtracking linesearch (3.1) converges to a stationary point.

5.2 Sampled LSR1

We derive convergence results for the sampled SR1 method for general nonconvex objective functions.

In order to establish convergence results one needs to ensure that the SR1 Hessian update equation (2.5) is well defined. To this end, we employ a cautious updating mechanism that is well suited to our particular algorithm. At the $k$th iteration, we update the Hessian approximation using only the set of curvature pairs that satisfy

[TABLE]

where $\epsilon>0$ is a predetermined constant. If no curvature pairs satisfy (A.4), then the new Hessian approximation is set to $B_{k}=I$. It is not trivial to test this condition in practice without explicitly constructing $d\times d$ matrices. We discuss this in detail in Section 7; see Appendix A.4 for more details.

For the analysis in this section, we make the following assumption in addition to 1, 3 and 4.

Assumption 5**.**

For all $k$,

[TABLE]

where $\xi\in(0,1)$ and $\beta_{k}=1+\|B_{k}\|$.

Assumption 5 ensures that at every iteration we solve the trust-region subproblem sufficiently accurately.

We prove that the Hessian approximations $B_{k}$ generated by the S-LSR1 method are uniformly bounded from above. The proof technique is an adaptation of that in [45]; however, modifications are necessary since the Hessian approximations are constructed using information only from the current iterate, and not constructed sequentially.

Lemma 5.8**.**

Suppose that Assumptions 1, 4 and 5 hold. Let $\{B_{k}\}$ be the Hessian approximations generated by Algorithm 3, with the modification that the approximation update is performed using only curvature pairs that satisfy (A.4), for some $\epsilon>0$, and $B_{k}=I$ if no curvature pairs satisfy (A.4). Then, there exists a constant $\nu_{2}>0$ such that

[TABLE]

Proof.

As in the proof of Lemma 5.4, note that there is a chance that no curvature pairs are selected in Algorithm 1. In this case, the Hessian approximation is $B_{k}=I$, and thus $\nu_{2}=1$ and condition (5.32) is satisfied.

We now consider the case where at least one curvature pair is selected by Algorithm 1. In this case, the sampled LSR1 updating formula is given as follows. Let $\tilde{m}_{k}\in\{1,...,m\}$ denote the number of curvature pairs that satisfy (A.4) at the $k$th iteration, where $m$ is the memory. At the $k$th iteration, given a set of curvature pairs $(s_{k,j},y_{k,j})$, for $j=1,\dots,\tilde{m}_{k}$

Set $B_{k}^{(0)}=\gamma_{k}I$, where $0\leq\gamma_{k}<\gamma$. 2. 2.

For $i=1,\dots,\tilde{m}_{k}$ compute

[TABLE] 3. 3.

Set $B_{k+1}=B_{k}^{(\tilde{m}_{k})}$.

In our algorithm (Algorithm 1), there are two options for constructing the curvature pairs $s_{k,j}$ and $y_{k,j}$. At the current iterate $w_{k}$ we sample points $\bar{w}_{j}$ for $j=1,\dots,m$ and set

[TABLE]

Given a set of $\tilde{m}_{k}$ curvature pairs that satisfy (A.4), we now prove an upper bound for $\|B_{k}\|$. We first prove the bound for a given iteration $k$ and for all updates to the Hessian approximation $i=0,1,\dots,\tilde{m}_{k}$ ($\|B_{k}^{i}\|$), and then get an upper bound for all $k$ ($\|B_{k}\|$).

For a given iteration $k$, we prove a bound on $\|B_{k}^{i}\|$ via induction, and show

[TABLE]

where $\bar{\gamma}_{k}$ is such that $\|\nabla^{2}F(w_{k})\|\leq\bar{\gamma}_{k}$, and whose existence follows from Assumption 4 ($\bar{\gamma}_{k}\leq L<\infty$). For $i=0$, the bound holds trivially since $B_{k}^{(0)}=\gamma_{k}I$. Now assume that (5.35) holds true for some $i\geq 0$. Note that all the curvature pairs that are used in the update of the Hessian approximation satisfy (A.4). By the definition of the SR1 updates, we have for some index $i+1$ that

[TABLE]

and thus

[TABLE]

where the first inequality is due to the application of the triangle inequality, the second inequality is due to condition (A.4), the fourth inequality is due to the application of the triangle inequality, and the fifth inequality is due to application of Cauchy’s inequality and in the last inequality we used that $\bar{\gamma_{k}}\geq\bar{\gamma}_{k,i+1}=\frac{\|y_{k,i+1}\|}{\|s_{k,i+1}\|}>0$. Substituting (5.35),

[TABLE]

which completes the inductive proof. Thus, for any $k$ we have an upper bound on the Hessian approximation. Therefore, since $B_{k+1}=B_{k}^{(\tilde{m}_{k})}$, the sampled SR1 Hessian approximation constructed at the $k$th iteration satisfies

[TABLE]

Now we generalize the result for all iterations $k$. For $k=0$, the bound holds trivially, since the first step of the sampled LSR1 method is a gradient method ($B_{0}=I$). For $k\geq 1$, we assume that $\gamma_{k}\leq\gamma<\infty$ and $\bar{\gamma}_{k}\leq\bar{\gamma}\leq L<\infty$ for all $k$, and thus

[TABLE]

for some $\nu_{2}>0$. This completes the proof. ∎

Utilizing Lemma 5.8, we show that the S-LSR1 with a cautious updating strategy converges. In order to prove the following result, we make use of well-known results for Trust-Region methods; see [17]. As such, the proof is identical to [17, Theorem 6.4.5]; to keep the paper self contained and due to the notation differences we include the proof.

Theorem 5.9**.**

Suppose that Assumptions 1, 3, 4 and 5 hold. Let $\{w_{k}\}$ be the iterates generated by Algorithm 3, with the modification that the Hessian approximation update is performed using only curvature pairs that satisfy (A.4), for some $\epsilon>0$, and $B_{k}=I$ if no curvature pairs satisfy (A.4). Then,

[TABLE]

Proof.

Assume, for the purpose of establishing a contradiction, that there is a subsequence of successful iterations (where $\rho_{k}>\eta_{1}$, Line 6, Algorithm 3), indexed by $t_{i}\subseteq\mathcal{S}$ where $\mathcal{S}=\{k\geq 0|\rho_{k}\geq\eta_{1}\}$, such that

[TABLE]

for some $\epsilon>0$ and for all $i$. Theorem 6.4.5 from [17] then ensures the existence for each $t_{i}$ of a first successful iteration $\ell(t_{i})>t_{i}$ such that

[TABLE]

Let $\ell_{i}=\ell(t_{i})$, we thus obtain that there is anotehr subsequence of $\mathcal{S}$ indexed by $\{\ell_{i}\}$ such that

[TABLE]

We now restrict our attention to the subsequence of successful iterations whose indices are in the set

[TABLE]

where $t_{i}$ and $\ell_{i}$ belong to the subsequences $\mathcal{S}$ and $\mathcal{K}$, respectively.

Using Assumption 5, the fact that $\mathcal{K}\subseteq\mathcal{S}$ and (5.37), we deduce that for $k\in\mathcal{K}$

[TABLE]

where we used the result of Lemma 5.8. Since the sequence $\{F(w_{k})\}$ is monotonically decreasing and bounded below (Assumption 3), it is convergent, and the left-hand-side of (5.38) must tend to zero as $k\rightarrow\infty$. Thus,

[TABLE]

As a consequence, the term containing $\Delta_{k}$ is the dominant term in the $\min$ (5.38) and we have, for $k\in\mathcal{K}$ sufficiently large,

[TABLE]

From this bound, we deduce that, for $i$ sufficiently large

[TABLE]

As a consequence of Assumption 3 and the monotonicity of the sequence $\{F(w_{k})\}$, we have that the right-hand-side of (5.41) must converge to zero, and thus $\|w_{t_{i}}-w_{\ell_{i}}\|\rightarrow 0$ as $i\rightarrow\infty$.

By continuity of the gradient (Assumption 1), we thus deduce that $\|\nabla F(w_{t_{i}})-\nabla F(w_{\ell_{i}})\|\rightarrow 0$. However, this is impossible because of the definitions of $\{t_{i}\}$ and $\{\ell_{i}\}$, which imply that $\|\nabla F(w_{t_{i}})-\nabla F(w_{\ell_{i}})\|\geq\delta$. Hence, no subsequence satisfying (5.36) can exist, and the theorem is proved. ∎

Theorem 5.9 shows that the sampled SR1 method converges to a stationary point. This result is similar in nature to that of the LSR1 method; see [45].

6 Probabilistic Bounds on Sampled Quasi-Newton Methods

Since our proposed methods randomly select $m$ curvature pairs ($\{s,y\}$) at every iteration and we require the pairs satisfy certain conditions ((5.19) and (A.4) for S-LBFGS and S-LSR1, respectively), a fair question to ask is how many pairs are accepted and used to construct Hessian approximations at every iteration. Alternatively, the question can be posed as what is the probability that a given random $\{s,y\}$ pair satisfies the required conditions and is used in the quasi-Newton Hessian approximations.

In this section, we present probabilistic bounds that illustrate the probability of accepting a given $\{s,y\}$ pair. To do this, we leverage the form of Option II for computing the $y$ vectors (given a vector $s$) and the fact that $s$ can be any random vector. We will assume throughout this section that $s$ is uniformly sampled on a unit sphere, i.e., $s\sim\mathcal{U}(\mathcal{S}(0,1))$. We also illustrate the probabilities of accepting pairs empirically for synthetic problems with different dimensions and acceptance tolerances, and on two toy classification problem.

6.1 Probabilistic Bounds for S-LBFGS

In this section, we present results that illustrate the probability that the pairs generated within the S-LBFGS method satisfy (5.19). We first derive an expression for the probability of accepting a pair $\{s,y\}$, and then provide some empirical evidence to show the probability of accepting pairs for different problems.

By Option II, (5.19) can be expressed as

[TABLE]

for any $w\in\mathbb{R}^{d}$. Notice that the middle term above is the Raleigh quotient of a random vector $s$ with respect to the Hessian matrix. Thus, for any $w\in\mathbb{R}^{d}$ and any given random vector $s$, we are interested in the following probability,

[TABLE]

The following theorem gives an expression for the probability of accepting the pair $\{s,y\}$.

Theorem 6.1**.**

Let $\lambda=(\lambda_{1},\lambda_{2},\dots,\lambda_{d})$ be the eigenvalues of the true Hessian at some point $w\in\mathbb{R}^{d}$ ($\nabla^{2}F(w)$), $s\in\mathbb{R}^{d}$ be a random vector uniformly distributed on a sphere, and $\epsilon>0$ be a prescribed tolerance. Then,

[TABLE]

Proof.

We refer interested readers to [6, Theorem 9] for the proof of this theorem. ∎

Although the result of Theorem 6.1 is interesting, a reasonable criticism is that it requires knowledge of the eigenvalues of the Hessian matrix, something that is prohibitively expensive to compute for many deep learning training problems. However, we show numerically that for certain neural network problems, the probability of accepting pairs is relatively high. Specifically, we present empirical results that illustrate the probability that a given pair $\{s,y\}$ is accepted for different problems and different $\epsilon$. The problems considered are summarized in Table 3, and the results are given in Figure 7.

The first two problems (Table 3) have synthetic eigenvalue distributions. The goal of these problems is to investigate the effect of $n$ and $\epsilon$ on the probability. As is clear, for the first problem, the probability of accepting a random $\{s,y\}$ pair is around $50\%$ and the probability decreases last as $\epsilon\rightarrow 1$, which is not surprising due to the eigenvalue structure. For the second problem, where the eigenvalues are smaller, the probability becomes almost zero for $\epsilon\leq 10^{-4}$. For the last two problems, we considered the structures of the toy classification problems (small and medium) for different points in parameter space. Note, points A, B and C are the same as those used in Figures 2 and 3. As is clear from the Figures 7c and 7d, the probability of accepting a random curvature pair is very high as long as $\epsilon$ is not too large. The main takeaway from these numerical results is that the probability of accepting curvature pairs is relatively large as long as the tolerance is not chosen to be too large (in practice $\epsilon\approx 10^{-4}-10^{-8}$).

6.2 Probabilistic Bounds on S-LSR1

In this section, we present results that illustrate the probability that the pairs generated within the S-LSR1 method satisfy (A.4). We first derive an expression for the probability of accepting a pair $\{s,y\}$, and then provide some empirical evidence to show the probability of accepting pairs for different problems.

By Option II, (A.4) can be expressed as

[TABLE]

for any $w\in\mathbb{R}^{d}$, where the matrix $B$ is some SR1 Hessian approximation. Clearly, the acceptance of a new pair $\{s,y\}$ depends on the matrix $B$. To be more precise, for some $w\in\mathbb{R}^{d}$, given a new pair $\{s_{j},y_{j}\}$ for $j=1,\dots,m$, (6.2) can be expressed as

[TABLE]

where $B^{(j)}$ is an SR1 Hessian approximation constructed using all the pairs $\{s_{i},y_{i}\}_{i<j}$, and $B^{(0)}$ is the initial SR1 Hessian approximation (potentially $B^{(0)}=0$). As is clear, the acceptance of the new pair $\{s_{j},y_{j}\}$ depends recursively on all previously accepted curvature pairs. Similar to the S-LBFGS case, the left-hand-side of (6.3) is a Raleigh quotient of a random vector $s_{j}$ with respect to the Hessian matrix and the matrix $B^{(j-1)}$.

Thus, for any $w\in\mathbb{R}^{d}$, $B^{(j-1)}\in\mathbb{R}^{d\times d}$ and any given random vector $s_{j}\in\mathbb{R}^{d}$, we are interested in the following probability,

[TABLE]

The following theorem gives an expression for the probability of accepting the pair $\{s_{j},y_{j}\}$.

Theorem 6.2**.**

Let $\tilde{\lambda}^{j-1}=(\tilde{\lambda}_{1}^{j-1},\tilde{\lambda}_{2}^{j-1},\dots,\tilde{\lambda}_{d}^{j-1})$ be the eigenvalues of the matrix $\nabla^{2}F(w)-B^{(j-1)}$, $s_{j}\in\mathbb{R}^{d}$ be a random vector uniformly distributed on a sphere, and $\epsilon>0$ be a prescribed tolerance. Then, for all $j\in\{1,\dots,m\}$,

[TABLE]

Proof.

The proof of this theorem is an adaptation of [6, Theorem 9]. Note that $\mathbb{P}[|X|\leq\eta]=\mathbb{P}[-\eta\leq X\leq\eta]=\mathbb{P}[X\leq\eta]-\mathbb{P}[X\leq-\eta]$, and $\mathbb{P}[|X|>\eta]=1-\mathbb{P}[|X|\leq\eta]$. ∎

As in the case for S-LBFGS, we now illustrate the probability of accepting pairs empirically. We conducted the same set of experiments as in Section 6.1; see Table 3 for details. The probability of accpeting pairs for the synthetic problems is larger for S-LSR1 than S-LBFGS. This is due to the fact that negative values of the Rayleigh quotient are also accepted, as long as they are large enough in magnitude. That being said, the relative performance when the eigenvalues are chosen to be smaller is similar to the S-LBFGS method. For the toy classification problems, as for the S-LBFGS method, the probability of accepting curvature pairs is close to 100% as long as $\epsilon$ is chosen to be small.

7 Numerical Experiments

In this section, we present numerical experiments on a toy classification problem as well as on popular benchmarking binary classification and neural network training tasks in order to illustrate the performance of our proposed sampled quasi-Newton methods.

7.1 Method Specifications and Details

Before we present the numerical results333All codes to reproduce results presented in this section are available at: https://github.com/OptMLGroup/SQN., we discuss the implementation details for all the methods. For ADAM [37], we tuned the steplength and batch size for each problem independently. For GD and BFGS-type methods, we computed the steplength using a backtracking Armijo line search [52]. For SR1-type methods, we solved the subproblems (2.4) using CG-Steihaug [52]. For BFGS and SR1, we constructed the full (inverse) Hessian approximations explicitly, whereas for the limited-memory we never constructed the full matrices. For limited-memory BFGS methods we used the two-loop recursion to get the search direction [52]. Implementing the limited memory SR1 methods is not trivial; we made use of the compact representations of the SR1 matrices [14] and computed the steps dynamically; see Appendix A.4 for details.

7.2 Toy Classification Problem

Consider the following simple classification problem, illustrated in Figure 9, consisting of two classes (red and blue) each with 50 data points. The goal of this classification task is to find a nonlinear decision boundary that separates the two classes. We trained three fully connected neural networks–small, medium and large–with sigmoid activation functions and 4 hidden layers. The details of the three networks are summarized in Table 4.

For this problem, we ran each method 100 times starting from different initial points and show the results for different budget levels. The results are summarized in Figure 10. In order to better visualize the relative performance of our proposed sampled quasi-Newton methods compared to the classical variants, we show accuracy vs. epochs plots in Figure 11. As is clear from the figures, the proposed methods outperform their classical variants as well as the first-order methods. See Appendix A.6 for more results.

For this problem, we ran each method 100 times starting from different initial points and show the results for different budget levels. The results are summarized in Figure 10. In order to better visualize the relative performance of our proposed sampled quasi-Newton methods compared to the classical variants, we show accuracy vs. epochs plots in Figure 11. As is clear from the figures, the proposed methods outperform their classical variants as well as the first-order methods. See Appendix A.6 for more results.

The toy classification problem is inherently complex. As is clear from the results, first-order methods (GD and ADAM) are not competitive with other reported methods, as they require a significant computational budget in order to achieve low classification error. It is worth noting that as we increase the size of the neural networks (over-parameterization), the performance of these methods becomes better. On the other hand, quasi-Newton methods have better performance in this complex, albeit small, problem, primarily due to the use of curvature information. Amongst the reported results, our sampled quasi-Newton methods significantly outperform the classical methods. We posit that this is the case due to the use of more recent and local curvature information in the updates.

7.3 Logistic Regression

Next we consider $\ell_{2}$-regularized logistic regression problems of the form

[TABLE]

where $(x_{i},y_{i})_{i=1}^{n}$ denote the training examples and $\lambda>0$ is the regularization parameter. We present results on two popular machine learning datasets (rcv1 and w8a; [15]); see Appendix A.7.1 for dataset details and more results. We compared the performance of the proposed sampled quasi-Newton methods with gradient descent (GD) and classical quasi-Newton methods (LSR1 and LBFGS). Figures 12 illustrates the performance of the methods in terms of optimality gap (training loss), training accuracy and testing accuracy. As is clear from Figure 12, the sample quasi-Newton methods are competitive with the classical variants in terms of all three metrics. One can also observed that in the initial stages of the optimization, it appears that the sampled quasi-Newton methods outperform their classical counterparts.

7.4 Nonlinear Least Squares

In this section, we consider nonlinear least squares problems [62] of the form

[TABLE]

where $(x_{i},y_{i})_{i=1}^{n}$ denote the training examples. We present results on the same datasets and compare against the same methods and using the same metrics as in Section 7.3. As is clear from Figure 13, the sampled quasi-Newton methods outperform their classical counterparts across the board. This is consistent with the results on other datasets; see Appendix A.7.2. We posit that this is due to the fact that sampling curvature pairs at every iterations allows for the method to incorporate more recent, local and reliable curvature information, a feature that can be indispensable in the nonconvex setting.

7.5 Neural Network Training: MNIST and CIFAR10

We illustrate the performance of the sampled quasi-Newton methods on standard benchmarking neural network training tasks: MNIST [40] and CIFAR10 [38]. The details of the problems are given in Table 5. For these problems we used sigmoid and softplus activation functions and softmax cross-entropy loss. For the memory-variant algorithms, we considered the memory from the set $m\in\{4,16,64,256\}$, and report the best performance with respect to the different memory sizes. The results of these experiments are given in Figures 14 and 15.

Overall, the sampled quasi-Newton methods outperform their classical variants. We should note that the goal of these experiments is not to perform better than ADAM, rather the performance of ADAM can be viewed as a benchmark. The reasons for this are two-fold. First, ADAM is a stochastic algorithm while the other reported methods are deterministic. Second, we report results for the best hyper-parameter settings for ADAM (well-tuned); see Appendix A.8, while the other methods do not require tuning or they are insensitive to the choice of hyper-parameters.

For the MNIST problems, the S-LSR1 method is able to achieve comparable accuracy to that of well-tuned ADAM, after a lot more epochs. That being said, in a distributed setting, the time to perform one iteration (one epoch) of S-LSR1 is significantly smaller than the time to perform one epoch of ADAM, and as such in terms of Wall Clock Time, the proposed method could be more efficient. With regards to the CIFAR10 problems, one can observe that our proposed sampled methods perform on par if not better than classical quasi-Newton methods. We posit that the reason that S-LSR1 has better performance than S-LBFGS is due to the possible utilization of negative curvature in the updates.

8 Final Remarks and Future Work

This paper describes two novel quasi-Newton methods; S-LBFGS and S-LSR1. Contrary to classical quasi-Newton methods, these methods forget past curvature information and sample new curvature information at every iteration. Numerical results show that the methods are efficient in practice, and the convergence guarantees of the methods match those of the classical variants.

Our algorithms can be extended to the stochastic setting where gradients and/or Hessians are computed inexactly. Moreover, the algorithms could be made adaptive following the ideas from [48, 29]. Furthermore, stronger theoretical (e.g., superlinear convergence) results could be proven for some variants of the sampled quasi-Newton methods. Finally, a large-scale numerical investigation of the method would test the limits of these methods.

Acknowledgements

This work was partially supported by the U.S. National Science Foundation, under award numbers NSF:CCF:1618717 and NSF:CCF:1740796, DARPA Lagrange award HR-001117S0039, and XSEDE Startup grant IRI180020.

Appendix A Additional Numerical Experiments and Method Details

In this section, we present additional numerical results and expand on some details about the methods.

A.1 Motivation Figure

In this section, we present more motivating plots showing the accuracy vs. iterations and accuracy vs. epochs for a toy classification problem. In the following experiments, we ran each method from 10 different initial points.

A.2 Eigenvalue Figures

In this section, we describe the procedure in which Figures 17 and 18 (Figures 2 & 3) were constructed. We plot the same figures below for ease of exposition.

To calculate the eigenvalues for SR1, LSR1 and S-LSR1 we used the following procedure.

We ran the SR1 method for T iterations on a toy classification problem. During the optimization, we computed the eigenvalues of the SR1 Hessian approximation at several points (e.g., A, B and C); black $\times$ marks on plots. 2. 2.

We stored all the curvature pairs $\{s_{k},y_{k}\}_{k=1}^{T}$ and the iterates $\{w_{k}\}_{k=1}^{T}$. 3. 3.

We constructed the true Hessian at all iterations and computed the eigenvalues of the true Hessian; dark blue $\bullet$ (positive eigenvalues) and light blue $\bullet$ (negative eigenvalues) marks on plots. 4. 4.

We then computed the limited-memory SR1 Hessian approximations at several points (e.g., A, B and C) using the $m$ most recent pairs and computed the eigenvalues of the approximations; orange $\blacktriangledown$ marks on plots. 5. 5.

Finally, we used the iterate information at points A, B and C, sampled $m$ points at random around those iterates with sampling radius $r$, constructed the sampled LSR1 Hessian approximations and computed the eigenvalues of the approximations; red $\bullet$ marks on plots.

Note: for Figure 17 we used $T=40$, $m=16$ and $r=0.01$, and for Figure 18 we used $T=70$, $m=32$ and $r=0.01$.

As is clear, the eigenvalues of the sampled LSR1 Hessian approximations better match the eigenvalues of the true Hessian. Similar results were obtained for other problems and for different parameters $m$ and $r$.

A.3 Trust-Region Management Subroutine

In this section we present, in detail, the Trust-Region management subroutine ($\Delta_{k+1}=\texttt{adjustTR}(\Delta_{k},\rho_{k})$) that is used in Algorithm 3. See [52] for further details.

A.4 Hessian-Free Implementation of Limited-Memory SR1 Methods

In this section, we discuss the practical implementation of limited-memory SR1 methods where we need not construct the Hessian approximation $B_{k}$ explicitly. For the purpose of this discuss we focus on the S-LSR1 method, but a similar approach can be used for the LSR1 method too. To do so, we utilize the compact representation of the Hessian approximation discussed in [14] which is equivalent to (2.5), i.e.,

[TABLE]

The compact representation can be expressed as follows:

[TABLE]

where $S_{k}=[s_{k,1},s_{k,2},\dots,s_{k,m}]\in\mathbb{R}^{d\times m}$ and $Y_{k}=[y_{k,1},y_{k,2},\dots,y_{k,m}]\in\mathbb{R}^{d\times m}$, and $B_{k}^{(0)}$ is a symmetric positive definite initial Hessian approximation, which for the purpose of this discussion we assume has the form $B_{k}^{(0)}=\gamma_{k}I$ ($0\leq\gamma_{k}<\gamma<\infty$). In (A.1), $D_{k}$ and $L_{k}$ are two $m\times m$ matrices that are defined as follows,

[TABLE]

[TABLE]

The curvature pairs in the matrices $S_{k}$ and $Y_{k}$ are pairs that satisfy the condition given in (A.4), i.e.,

[TABLE]

where $\epsilon>0$ is a predetermined constant.

In large-scale applications, it is not memory-efficient, or even possible for some applications, to store a $d\times d$ Hessian approximation matrix $B_{k+1}$. Instead, we can calculate the Hessian vector product $B_{k+1}v$, for some $v\in\mathbb{R}^{d}$, by leveraging the compact form of $B_{k}$ in (A.1) as follows:

[TABLE]

The above, $B_{k+1}v$, is very efficient in terms of memory, and even more importantly efficient to compute; the complexity of computing $B_{k+1}v$ is $\mathcal{O}(m^{2}d)$.

In Algorithm 3, we need to compute and use $B_{k+1}v$ in the following parts: (1) in checking the condition (A.4); (2) as part of the computation of solving the subproblem (2.4), i.e.,

[TABLE]

using the CG solver (see [52]); and, (3) in the calculation of $\rho_{k}$.

In the remainder of this section we describe the steps for checking whether the curvature pairs constructed by Algorithm 1 satisfy (A.4). This is by no means a trivial task; several researchers have proposed mechanisms for doing this [10, 45] by using spectral decomositions of $B_{k}$. We propose to do this in a dynamic manner leveraging (A.5).

The condition that we want to check (A.4) has the following form:

[TABLE]

for $i=1,\dots,m$, where $B_{k}^{(i)}$ are constructed recursively via,

[TABLE]

$S_{k}^{i}=[s_{k,1},s_{k,2},\dots,s_{k,i}]$, $Y_{k}^{i}=[y_{k,1},y_{k,2},\dots,y_{k,i}]$, and $D_{k}^{i}$ and $L_{k}^{i}$ are defined in equations (A.2) and (A.3), respectively, using $S_{k}^{i}$ and $Y_{k}^{i}$. Of course we want to check condition (A.6) without explicitly forming the matrices $B_{k}^{(i)}$, and instead construct $B_{k}^{(i)}s_{k,i+1}$ directly. To this end, by using (A.5), for any $i=0,\dots,m-1$ we have:

[TABLE]

where the matrices $S_{k}^{i}$, $Y_{k}^{i}$, $D_{k}^{i}$ and $L_{k}^{i}$ are defined as above.

The steps for checking (A.4) are as follows:

Consider the $k$th iteration of Algorithm 3, where the pairs $\bar{S}_{k}=[s_{k,1},\dots,s_{k,m}]$ and $\bar{Y}_{k}=[y_{k,1},\dots,y_{k,m}]$ and $B_{k}^{0}$ are constructed by Algorithm 1. Let $S_{k}=[\;]$ and $Y_{k}=[\;]$ be two empty matrices. 2. 2.

For any $i=1,\dots,m$, consider the pair $(s_{k,i},y_{k,i})$ and compute $B_{k}^{(i-1)}s_{k,i}$ by using (A.8) and the updated lists $S_{k}$ and $Y_{k}$. Note that for $i=1$ the matrices $S_{k}$ and $Y_{k}$ are empty and $B_{k}^{(0)}$ is the initial Hessian approximation, thus condition (A.6) can be checked directly.

1. (a)

If condition (A.6) is satisfied for this pair, then $S_{k}=[S_{k}\;s_{k,i}]$ add $Y_{k}=[Y_{k}\;y_{k,i}]$ 2. (b)

Else, discard the pair $(s_{k,i},y_{k,i})$

Using the mechanism described above, we recursively check condition (A.6), and construct well defined pairs $S_{k}$ and $Y_{k}$ which are used for the calculation of $p_{k}$ and $\rho_{k}$. As mentioned above, we can use the same idea for implementing the LSR1 method. Finally, we should note that the idea described above was successfully implemented in a distributed variant on the S-LSR1 method in [30].

A.5 Cost of communication

In this section, we show experiments conducted on a HPC cluster using a Cray Aries High Speed Network. The bandwidth ranges depending on the distance between nodes. We compiled the C++ code with the provided cray compiler.

In Figure 19, we show how the duration (seconds) of Broadcast and Reduce increases when vectors of longer length are processed.

In Figure 20, we show how long it takes (seconds) to perform Broadcast and Reduce operations for vectors of a given length if performed on different numbers of MPI processes. We have performed each operation 100 times and are showing the average time and 95% confidence intervals.

A.6 Toy Example

In this section, we present additional numerical results for the toy classification problem described in Section 7.2. In the following experiments, we ran each method from 100 different initial points.

A.6.1 Performance of Methods on small, medium and large toy classification problems - Box-plots

The following box-plots show the accuracy achieved by different methods for different budgets (epochs) and iterations.

A.6.2 Performance of Methods on small, medium and large toy classification problems

In this section, we present more experiments showing accuracy vs. iterations and accuracy vs. epochs for different methods on toy classification problem.

A.6.3 Comparison of BFGS-type methods

In this section, we present more experiments showing the accuracy achieved in terms of iterations and epochs for BFGS-type methods on toy classification problem.

A.6.4 Comparison of SR1-type methods

In this section, we present more experiments showing the accuracy achieved in terms of iterations and epochs for SR1-type methods on toy classification problem.

A.7 Binary Classification

In this section, we present additional numerical results for the binary classification tasks Sections 7.3 & 7.4. We show results on the following datasets444Datasets are available at https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/.