A Multi-Batch L-BFGS Method for Machine Learning

Introduction

It is common in machine learning to encounter optimization problems involving millions of parameters and very large datasets. To deal with the computational demands imposed by such applications, high performance implementations of stochastic gradient and batch quasi-Newton methods have been developed . In this paper we study a batch approach based on the L-BFGS method that strives to reach the right balance between efficient learning and productive parallelism.

In supervised learning, one seeks to minimize empirical risk,

At present, the preferred optimization method is the stochastic gradient descent (SGD) method , and its variants , which are implemented either in an asynchronous manner (e.g. when using a parameter server in a distributed setting) or following a synchronous mini-batch approach that exploits parallelism in the gradient evaluation . A drawback of the asynchronous approach is that it cannot use large batches, as this would cause updates to become too dense and compromise the stability and scalability of the method . As a result, the algorithm spends more time in communication as compared to computation. On the other hand, using a synchronous mini-batch approach one can achieve a near-linear decrease in the number of SGD iterations as the mini-batch size is increased, up to a certain point after which the increase in computation is not offset by the faster convergence .

An alternative to SGD is a batch method, such as L-BFGS, which is able to reach high training accuracy and allows one to perform more computation per node, so as to achieve a better balance with communication costs . Batch methods are, however, not as efficient learning algorithms as SGD in a sequential setting . To benefit from the strength of both methods some high performance systems employ SGD at the start and later switch to a batch method .

Multi-Batch Method. In this paper, we follow a different approach consisting of a single method that selects a sizeable subset (batch) of the training data to compute a step, and changes this batch at each iteration to improve the learning abilities of the method. We call this a multi-batch approach to differentiate it from the mini-batch approach used in conjunction with SGD, which employs a very small subset of the training data. When using large batches it is natural to employ a quasi-Newton method, as incorporating second-order information imposes little computational overhead and improves the stability and speed of the method. We focus here on the L-BFGS method, which employs gradient information to update an estimate of the Hessian and computes a step in $O(d)$ flops, where $d$ is the number of variables. The multi-batch approach can, however, cause difficulties to L-BFGS because this method employs gradient differences to update Hessian approximations. When the gradients used in these differences are based on different data points, the updating procedure can be unstable. Similar difficulties arise in a parallel implementation of the standard L-BFGS method, if some of the computational nodes devoted to the evaluation of the function and gradient are unable to return results on time — as this again amounts to using different data points to evaluate the function and gradient at the beginning and the end of the iteration. The goal of this paper is to show that stable quasi-Newton updating can be achieved in both settings without incurring extra computational cost, or special synchronization. The key is to perform quasi-Newton updating based on the overlap between consecutive batches. The only restriction is that this overlap should not be too small, something that can be achieved in most situations.

Contributions. We describe a novel implementation of the batch L-BFGS method that is robust in the absence of sample consistency; i.e., when different samples are used to evaluate the objective function and its gradient at consecutive iterations. The numerical experiments show that the method proposed in this paper — which we call the multi-batch L-BFGS method — achieves a good balance between computation and communication costs. We also analyze the convergence properties of the new method (using a fixed step length strategy) on both convex and nonconvex problems.

The Multi-Batch Quasi-Newton Method

In a pure batch approach, one applies a gradient based method, such as L-BFGS , to the deterministic optimization problem (1.1). When the number $n$ of training examples is large, it is natural to parallelize the evaluation of $F$ and $\nabla F$ by assigning the computation of the component functions $f_{i}$ to different processors. If this is done on a distributed platform, it is possible for some of the computational nodes to be slower than the rest. In this case, the contribution of the slow (or unresponsive) computational nodes could be ignored given the stochastic nature of the objective function. This leads, however, to an inconsistency in the objective function and gradient at the beginning and at the end of the iteration, which can be detrimental to quasi-Newton methods. Thus, we seek to find a fault-tolerant variant of the batch L-BFGS method that is capable of dealing with slow or unresponsive computational nodes.

A similar challenge arises in a multi-batch implementation of the L-BFGS method in which the entire training set $T=\{(x^{i},y^{i})_{i=1}^{n}\}$ is not employed at every iteration, but rather, a subset of the data is used to compute the gradient. Specifically, we consider a method in which the dataset is randomly divided into a number of batches — say 10, 50, or 100 — and the minimization is performed with respect to a different batch at every iteration. At the $k$-th iteration, the algorithm chooses a batch $S_{k}\subset\{1,\ldots,n\}$, computes

and takes a step along the direction $-H_{k}g_{k}^{S_{k}}$, where $H_{k}$ is an approximation to $\nabla^{2}F(w_{k})^{-1}$. Allowing the sample $S_{k}$ to change freely at every iteration gives this approach flexibility of implementation and is beneficial to the learning process, as we show in Section 4. (We refer to $S_{k}$ as the sample of training points, even though $S_{k}$ only indexes those points.)

The case of unresponsive computational nodes and the multi-batch method are similar. The main difference is that node failures create unpredictable changes to the samples $S_{k}$, whereas a multi-batch method has control over sample generation. In either case, the algorithm employs a stochastic approximation to the gradient and can no longer be considered deterministic. We must, however, distinguish our setting from that of the classical SGD method, which employs small mini-batches and noisy gradient approximations. Our algorithm operates with much larger batches so that distributing the function evaluation is beneficial and the compute time of $g_{k}^{S_{k}}$ is not overwhelmed by communication costs. This gives rise to gradients with relatively small variance and justifies the use of a second-order method such as L-BFGS.

Robust Quasi-Newton Updating. The difficulties created by the use of a different sample $S_{k}$ at each iteration can be circumvented if consecutive samples $S_{k}$ and $S_{k+1}$ overlap, so that $O_{k}=S_{k}\cap S_{k+1}\neq\emptyset.$ One can then perform stable quasi-Newton updating by computing gradient differences based on this overlap, i.e., by defining

in the notation given in (2.2). The correction pair $(y_{k},s_{k})$ can then be used in the BFGS update. When the overlap set $O_{k}$ is not too small, $y_{k}$ is a useful approximation of the curvature of the objective function $F$ along the most recent displacement, and will lead to a productive quasi-Newton step. This observation is based on an important property of Newton-like methods, namely that there is much more freedom in choosing a Hessian approximation than in computing the gradient . Thus, a smaller sample $O_{k}$ can be employed for updating the inverse Hessian approximation $H_{k}$ than for computing the batch gradient $g_{k}^{S_{k}}$ in the search direction $-H_{k}g_{k}^{S_{k}}$. In summary, by ensuring that unresponsive nodes do not constitute the vast majority of all working nodes in a fault-tolerant parallel implementation, or by exerting a small degree of control over the creation of the samples $S_{k}$ in the multi-batch method, one can design a robust method that naturally builds upon the fundamental properties of BFGS updating.

We should mention in passing that a commonly used strategy for ensuring stability of quasi-Newton updating in machine learning is to enforce gradient consistency , i.e., to use the same sample $S_{k}$ to compute gradient evaluations at the beginning and the end of the iteration. Another popular remedy is to use the same batch $S_{k}$ for multiple iterations , alleviating the gradient inconsistency problem at the price of slower convergence. In this paper, we assume that achieving such sample consistency is not possible (in the fault-tolerant case) or desirable (in a multi-batch framework), and wish to design a new variant of L-BFGS that imposes minimal restrictions in the sample changes.

At the $k$-th iteration, the multi-batch BFGS algorithm chooses a set $S_{k}\subset\{1,\ldots,n\}$ and computes a new iterate

where $\alpha_{k}$ is the step length, $g_{k}^{S_{k}}$ is the batch gradient (2.2) and $H_{k}$ is the inverse BFGS Hessian matrix approximation that is updated at every iteration by means of the formula

To compute the correction vectors $(s_{k},y_{k})$, we determine the overlap set $O_{k}=S_{k}\cap S_{k+1}$ consisting of the samples that are common at the $k$-th and $k+1$-st iterations. We define

and compute the correction vectors as in (2.3). In this paper we assume that $\alpha_{k}$ is constant.

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 a set of $m$ correction pairs $\{s_{i},y_{i}\}$ kept in storage. The memory parameter $m$ is typically in the range 2 to 20. When computing the matrix-vector product in (2.4) it is not necessary to form that matrix $H_{k}$ since one can obtain this product via the two-loop recursion , using the $m$ most recent correction pairs $\{s_{i},y_{i}\}$. After the step has been computed, the oldest pair $(s_{j},y_{j})$ is discarded and the new curvature pair is stored.

A pseudo-code of the proposed method is given below, and depends on several parameters. The parameter $r$ denotes the fraction of samples in the dataset used to define the gradient, i.e., $r=\frac{\left|S\right|}{n}$. The parameter $o$ denotes the length of overlap between consecutive samples, and is defined as a fraction of the number of samples in a given batch $S$, i.e., $o=\frac{\left|O\right|}{\left|S\right|}$.

2 Sample Generation

We now discuss how the sample $S_{k+1}$ is created at each iteration (Line 8 in Algorithm 1).

Let $\mathcal{J}_{k}\subset\{1,2,...,B\}$ and $\mathcal{J}_{k+1}\subset\{1,2,...,B\}$ be the set of indices of all nodes that returned a gradient at the $k$-th and $k+1$-st iterations, respectively. Using this notation $S_{k}=\cup_{j\in\mathcal{J}_{k}}\mathcal{B}_{j}$ and $S_{k+1}=\cup_{j\in\mathcal{J}_{k+1}}\mathcal{B}_{j}$, and we define $O_{k}=\cup_{j\in\mathcal{J}_{k}\cap\mathcal{J}_{k+1}}\mathcal{B}_{j}$. The simplest implementation in this setting preallocates the data on each compute node, requiring minimal data communication, i.e., only one data transfer. In this case the samples $S_{k}$ will be independent if node failures occur randomly. On the other hand, if the same set of nodes fail, then sample creation will be biased, which is harmful both in theory and practice. One way to ensure independent sampling is to shuffle and redistribute the data to all nodes after a certain number of iterations.

Multi-batch Sampling. We propose two strategies for the multi-batch setting.

Figure 1b illustrates the sample creation process in the first strategy. The dataset is shuffled and batches are generated by collecting subsets of the training set, in order. Every set (except $S_{0}$) is of the form $S_{k}=\{O_{k-1},N_{k},O_{k}\}$, where $O_{k-1}$ and $O_{k}$ are the overlapping samples with batches $S_{k-1}$ and $S_{k+1}$ respectively, and $N_{k}$ are the samples that are unique to batch $S_{k}$. After each pass through the dataset, the samples are reshuffled, and the procedure described above is repeated. In our implementation samples are drawn without replacement, guaranteeing that after every pass (epoch) all samples are used. This strategy has the advantage that it requires no extra computation in the evaluation of $g_{k}^{O_{k}}$ and $g_{k+1}^{O_{k}}$, but the samples $\{S_{k}\}$ are not independent.

The second sampling strategy is simpler and requires less control. At every iteration $k$, a batch $S_{k}$ is created by randomly selecting $\left|S_{k}\right|$ elements from $\{1,\ldots n\}$. The overlapping set $O_{k}$ is then formed by randomly selecting $\left|O_{k}\right|$ elements from $S_{k}$ (subsampling). This strategy is slightly more expensive since $g_{k+1}^{O_{k}}$ requires extra computation, but if the overlap is small this cost is not significant.

Convergence Analysis

In this section, we analyze the convergence properties of the multi-batch L-BFGS method (Algorithm 1) when applied to the minimization of strongly convex and nonconvex objective functions, using a fixed step length strategy. We assume that the goal is to minimize the empirical risk $F$ given in (1.1), but note that a similar analysis could be used to study the minimization of the expected risk.

Due to the stochastic nature of the multi-batch approach, every iteration of Algorithm 1 employs a gradient that contains errors that do not converge to zero. Therefore, by using a fixed step length strategy one cannot establish convergence to the optimal solution $w^{\star}$, but only convergence to a neighborhood of $w^{\star}$ . Nevertheless, this result is of interest as it reflects the common practice of using a fixed step length and decreasing it only if the desired testing error has not been achieved. It also illustrates the tradeoffs that arise between the size of the batch and the step length.

In our analysis, we make the following assumptions about the objective function and the algorithm.

$F$ is twice continuously differentiable.

Note that Assumption $A.2$ implies that the entire Hessian $\nabla^{2}F(w)$ also satisfies

for some constants $\lambda,\Lambda>0$. Assuming that every sub-sampled function $F^{O}(w)$ is strongly convex is not unreasonable as a regularization term is commonly added in practice when that is not the case.

We begin by showing that the inverse Hessian approximations $H_{k}$ generated by the multi-batch L-BFGS method have eigenvalues that are uniformly bounded above and away from zero. The proof technique used is an adaptation of that in .

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

Utilizing Lemma 3.1, we show that the multi-batch L-BFGS method with a constant step length converges to a neighborhood of the optimal solution.

Suppose that Assumptions A.1-A.4 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 1 with $\alpha_{k}=\alpha\in(0,\frac{1}{2\mu_{1}\lambda})$, starting from $w_{0}$. Then for all $k\geq 0$,

The bound provided by this theorem has two components: (i) a term decaying linearly to zero, and (ii) a term identifying the neighborhood of convergence. Note that a larger step length yields a more favorable constant in the linearly decaying term, at the cost of an increase in the size of the neighborhood of convergence. We will consider again these tradeoffs in Section 4, where we also note that larger batches increase the opportunities for parallelism and improve the limiting accuracy in the solution, but slow down the learning abilities of the algorithm.

One can establish convergence of the multi-batch L-BFGS method to the optimal solution $w^{\star}$ by employing a sequence of step lengths $\{\alpha_{k}\}$ that converge to zero according to the schedule proposed by Robbins and Monro . However, that provides only a sublinear rate of convergence, which is of little interest in our context where large batches are employed and some type of linear convergence is expected. In this light, Theorem 3.2 is more relevant to practice.

2 Nonconvex case

The BFGS method is known to fail on noconvex problems . Even for L-BFGS, which makes only a finite number of updates at each iteration, one cannot guarantee that the Hessian approximations have eigenvalues that are uniformly bounded above and away from zero. To establish convergence of the BFGS method in the nonconvex case cautious updating procedures have been proposed . Here we employ a cautious strategy that is well suited to our particular algorithm; we skip the update, i.e., set $H_{k+1}=H_{k}$, if the curvature condition

is not satisfied, where $\epsilon>0$ is a predetermined constant. Using said mechanism we show that the eigenvalues of the Hessian matrix approximations generated by the multi-batch L-BFGS method are bounded above and away from zero (Lemma 3.3). The analysis presented in this section is based on the following assumptions.

$F$ is twice continuously differentiable.

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

Suppose that Assumptions B.1-B.2 hold and let $\epsilon>0$ be given. Let $\{H_{k}\}$ be the Hessian approximations generated by Algorithm 1, with the modification that $H_{k+1}=H_{k}$ whenever (3.5) is not satisfied. Then, there exist constants $0<\mu_{1}\leq\mu_{2}$ such that

We can now follow the analysis in [4, Chapter 4] to establish the following result about the behavior of the gradient norm for the multi-batch L-BFGS method with a cautious update strategy.

Suppose that Assumptions B.1-B.5 above hold, and let $\epsilon>0$ be given. Let $\{w_{k}\}$ be the iterates generated by Algorithm 1, with $\alpha_{k}=\alpha\in(0,\frac{\mu_{1}}{\mu_{2}^{2}\eta\Lambda})$, starting from $w_{0}$, and with the modification that $H_{k+1}=H_{k}$ whenever (3.5) is not satisfied. Then,

This result bounds the average norm of the gradient of $F$ after the first $L-1$ iterations, and shows that the iterates spend increasingly more time in regions where the objective function has a small gradient.

Numerical Results

In this Section, we present numerical results that evaluate the proposed robust multi-batch L-BFGS scheme (Algorithm 1) on logistic regression problems. Figure 2 shows the performance on the webspam datasetLIBSVM: https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary.html. , where we compare it against three methods: (i) multi-batch L-BFGS without enforcing sample consistency (L-BFGS), where gradient differences are computed using different samples, i.e., $y_{k}=g_{k+1}^{S_{k+1}}-g_{k}^{S_{k}}$; (ii) multi-batch gradient descent (Gradient Descent), which is obtained by setting $H_{k}=I$ in Algorithm 1; and, (iii) serial SGD, where at every iteration one sample is used to compute the gradient. We run each method with 10 different random seeds, and, where applicable, report results for different batch ($r$) and overlap ($o$) sizes. The proposed method is more stable than the standard L-BFGS method; this is especially noticeable when $r$ is small. On the other hand, serial SGD achieves similar accuracy as the robust L-BFGS method and at a similar rate (e.g., $r=1\%$), at the cost of $n$ communications per epochs versus $\frac{1}{r(1-o)}$ communications per epoch. Figure 2 also indicates that the robust L-BFGS method is not too sensitive to the size of overlap. Similar behavior was observed on other datasets, in regimes where $r\cdot o$ was not too small; see Appendix B.1. We mention in passing that the L-BFGS step was computed using the a vector-free implementation proposed in .

We also explore the performance of the robust multi-batch L-BFGS method in the presence of node failures (faults), and compare it to the multi-batch variant that does not enforce sample consistency (L-BFGS). Figure 3 illustrates the performance of the methods on the webspam dataset, for various probabilities of node failures $p\in\{0.1,0.3,0.5\}$, and suggests that the robust L-BFGS variant is more stable; see Appendix B.2 for further results.

Lastly, we study the strong and weak scaling properties of the robust L-BFGS method on artificial data (Figure 4). We measure the time needed to compute a gradient (Gradient) and the associated communication (Gradient+C), as well as, the time needed to compute the L-BFGS direction (L-BFGS) and the associated communication (L-BFGS+C), for various batch sizes ($r$). The figure on the left shows strong scaling of multi-batch LBFGS on a $d=10^{4}$ dimensional problem with $n=10^{7}$ samples. The size of input data is 24GB, and we vary the number of MPI processes, $K\in\{1,2,\dots,128\}$. The time it takes to compute the gradient decreases with $K$, however, for small values of $r$, the communication time exceeds the compute time. The figure on the right shows weak scaling on a problem of similar size, but with varying sparsity. Each sample has $10\cdot K$ non-zero elements, thus for any $K$ the size of local problem is roughly $1.5$GB (for $K=128$ size of data 192GB). We observe almost constant time for the gradient computation while the cost of computing the L-BFGS direction decreases with $K$; however, if communication is considered, the overall time needed to compute the L-BFGS direction increases slightly. For more details see Appendix C.

Conclusion

This paper describes a novel variant of the L-BFGS method that is robust and efficient in two settings. The first occurs in the presence of node failures in a distributed computing implementation; the second arises when one wishes to employ a different batch at each iteration in order to accelerate learning. The proposed method avoids the pitfalls of using inconsistent gradient differences by performing quasi-Newton updating based on the overlap between consecutive samples. Numerical results show that the method is efficient in practice, and a convergence analysis illustrates its theoretical properties.

The first two authors were supported by the Office of Naval Research award N000141410313, the Department of Energy grant DE-FG02-87ER25047 and the National Science Foundation grant DMS-1620022. Martin Takáč was supported by National Science Foundation grant CCF-1618717.

References

Appendix A Proofs and Technical Results

We first restate the Assumptions that we use in the Convergence Analysis section (Section 3). Assumption $A$ and $B$ are used in the strongly convex and nonconvex cases, respectively.

$F$ is twice continuously differentiable.

There exist positive constants $\hat{\lambda}$ and $\hat{\Lambda}$ such that

Note that Assumption $A.2$ implies that the entire Hessian $\nabla^{2}F(w)$ also satisfies

$F$ is twice continuously differentiable.

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

There exist constants $\gamma\geq 0$ and $\eta>0$ such that

A.2 Proof of Lemma 3.1 (Strongly Convex Case)

The following Lemma shows that the eigenvalues of the matrices generated by the multi-batch L-BFGS method are bounded above and away from zero if $F$ is strongly convex.

If Assumptions A.1-A.2 above hold, there exist constants $0<\mu_{1}\leq\mu_{2}$ such that the Hessian approximations $\{H_{k}\}$ generated by the multi-batch L-BFGS method (Algorithm 1) satisfy

Instead of analyzing the inverse Hessian approximation $H_{k}$, we study the direct Hessian approximation $B_{k}=H_{k}^{-1}$. In this case, the limited memory quasi-Newton updating formula is given as follows

The curvature pairs $s_{k}$ and $y_{k}$ are updated via the following formulae

A consequence of Assumption $A.2$ is that the eigenvalues of any sub-sampled Hessian ($\left|O\right|$ samples) are bounded above and away from zero. Utilizing this fact, the convexity of component functions and the definitions (A.12), we have

On the other hand, strong convexity of the sub-sampled functions, the consequence of Assumption $A.2$ and definitions (A.12), provide a lower bound,

Combining the upper and lower bounds (A.13) and (A.14)

The above proves that the eigenvalues of the matrices $B_{k}^{(0)}=\frac{y_{k}^{T}y_{k}}{s_{k}^{T}y_{k}}I$ at the start of the L-BFGS update cycles are bounded above and away from zero, for all $k$. We now use a Trace-Determinant argument to show that the eigenvalues of $B_{k}$ are bounded above and away from zero.

for some positive constant $C_{1}$, where the inequalities above are due to (A.15), 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 , the determinant of the matrix $B_{k+1}$ generated by the multi-batch L-BFGS method can be expressed as,

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}$ and Assumption $A.2$.

The trace (A.2) and determinant (A.2) 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. ∎

A.3 Proof of Theorem 3.2 (Strongly Convex Case)

Utilizing the result from Lemma 3.1, we now prove a linear convergence result to a neighborhood of the optimal solution, for the case where Assumptions $A$ hold.

Suppose that Assumptions A.1-A.4 above hold, and let $F^{\star}=F(w^{\star})$, where $w^{\star}$ is the minimizer of $F$. Let $\{w_{k}\}$ be the iterates generated by the multi-batch L-BFGS method (Algorithm 1) with

starting from $w_{0}$. Then for all $k\geq 0$,

where the first inequality arises due to (A.9), and the second inequality arises as a consequence of Lemma 3.1.

Taking the expectation (over $S_{k}$) of equation (A.3)

where in the first inequality we make use of Assumption $A.5$, and the second inequality arises due to Lemma 3.1 and Assumption $A.4$.

Since $F$ is $\lambda$-strongly convex, we can use the following relationship between the norm of the gradient squared, and the distance of the $k$-th iterate from the optimal solution.

where the expectation is over all batches $S_{0},S_{1},...,S_{k-1}$ and all history starting with $w_{0}$. Thus (A.3) can be expressed as,

from which we deduce that in order to reduce the value with respect to the previous function value, the step length needs to be in the range

Subtracting $\frac{\alpha\mu_{2}^{2}\gamma^{2}\Lambda}{4\mu_{1}\lambda}$ from either side of (A.22) yields

Finally using the definition of $\phi_{k}$ (A.21) with the above expression yields the desired result,

A.4 Proof of Lemma 3.3 (Nonconvex Case)

The following Lemma shows that the eigenvalues of the matrices generated by the multi-batch L-BFGS method are bounded above and away from zero (nonconvex case).

Suppose that Assumptions B.1-B.2 hold and let $\epsilon>0$ be given. Let $\{H_{k}\}$ be the Hessian approximations generated by the multi-batch L-BFGS method (Algorithm 1), with the modification that the Hessian approximation $H_{k}$ update is performed only when

else $H_{k+1}=H_{k}$. Then, there exist constants $0<\mu_{1}\leq\mu_{2}$ such that

Similar to the proof of Lemma 3.1, we study the direct Hessian approximation $B_{k}=H_{k}^{-1}$.

The curvature pairs $s_{k}$ and $y_{k}$ are updated via the following formulae

The skipping mechanism (A.25) provides both an upper and lower bound on the quantity $\frac{\|y_{k}\|^{2}}{y_{k}^{T}s_{k}}$, which in turn ensures that the initial L-BFGS Hessian approximation is bounded above and away from zero. The lower bound is attained by repeated application of Cauchy’s inequality to condition (A.25). We have from (A.25) that

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

Re-arranging the above expression yields the desired upper bound,

The above proves that the eigenvalues of the matrices $B_{k}^{(0)}=\frac{y_{k}^{T}y_{k}}{s_{k}^{T}y_{k}}I$ at the start of the L-BFGS update cycles are bounded above and away from zero, for all $k$. The rest of the proof follows the same trace-determinant argument as in the proof of Lemma 3.1, the only difference being that the last inequality in A.2 comes as a result of the cautious update strategy. ∎

A.5 Proof of Theorem 3.4 (Nonconvex Case)

Utilizing the result from Lemma 3.3, we can now establish the following result about the behavior of the gradient norm for the multi-batch L-BFGS method with a cautious update strategy.

Suppose that Assumptions B.1-B.5 above hold. Let $\{w_{k}\}$ be the iterates generated by the multi-batch L-BFGS method (Algorithm 1) with

where $w_{0}$ is the starting point. Also, suppose that if

for some ${\epsilon}>0$, the inverse L-BFGS Hessian approximation is skipped, $H_{k+1}=H_{k}$. Then, for all $k\geq 0$,

where the second inequality holds due to Assumption $B.4$, and the fourth inequality is obtained by using the upper bound on the step length. Taking an expectation over all batches $S_{0},S_{1},...,S_{k-1}$ and all history starting with $w_{0}$ yields

The left-hand-side of the above inequality is a telescoping sum

Substituting the above expression into (A.5) and re-arranging terms

Dividing the above equation by $L$ completes the proof. ∎

Appendix B Extended Numerical Experiments - Real Datasets

In this Section, we present further numerical results, on the datasets listed in Table 1, in both the multi-batch and fault-tolerant settings. Note, that some of the datasets are too small, and thus, there is no reason to run them on a distributed platform; however, we include them as they are part of the standard benchmarking datasets.

Notation. Let $n$ denote the number of training samples in a given dataset, $d$ the dimension of the parameter vector $w$, and $K$ the number of MPI processes used. The parameter $r$ denotes the fraction of samples in the dataset used to define the gradient, i.e., $r=\frac{\left|S\right|}{n}$. The parameter $o$ denotes the length of overlap between consecutive samples, and is defined as a fraction of the number of samples in a given batch $S$, i.e., $o=\frac{\left|O\right|}{\left|S\right|}$.

We focus on logistic regression classification; the objective function is given by

where $(x^{i},y^{i})_{i=1}^{n}$ denote the training examples and $\sigma=\frac{1}{n}$ is the regularization parameter.

For the experiments in this section (Figures 5-13), we run four methods:

(Robust L-BFGS) robust multi-batch L-BFGS (Algorithm 1),

(L-BFGS) multi-batch L-BFGS without enforcing sample consistency; gradient differences are computed using different samples, i.e., $y_{k}=g_{k+1}^{S_{k+1}}-g_{k}^{S_{k}}$,

(Gradient Descent) multi-batch gradient descent; obtained by setting $H_{k}=I$ in Algorithm 1,

(SGD) serial SGD; at every iteration one sample is used to compute the gradient.

In Figures 5-13 we show the evolution of $\|\nabla F(w)\|$ for different step lengths $\alpha$, and for various batch ($\left|S\right|=r\cdot n$) and overlap ($\left|O\right|=o\cdot\left|S\right|$) sizes. Each Figure (5-13) consists of 10 plots that illustrate the performance of the methods with the following parameters:

Top 3 plots: $\alpha=1$, $o=20\%$ and $r=1\%,5\%,10\%$

Middle 3 plots: $\alpha=0.1$, $o=20\%$ and $r=1\%,5\%,10\%$

Bottom 4 plots: $\alpha=1$, $r=1\%$ and $o=5\%,10\%,20\%,30\%$

As is expected for quasi-Newton methods, robust L-BFGS performs best with a step-size $\alpha=1$, for the most part.

B.2 Fault-tolerant L-BFGS Implementation

If we run a distributed algorithm, for example on a shared computer cluster, then we may experience delays. Such delays can be caused by other processes running on the same compute node, node failures and for other reasons. As a result, given a computational (time) budget, these delays may cause nodes to fail to return a value. To illustrate this behavior, and to motivate the robust fault-tolerant L-BFGS method, we run a simple benchmark MPI code on two different environments:

Amazon EC2 – Amazon EC2 is a cloud system provided by Amazon. It is expected that if load balancing is done properly, the execution time will have small noise; however, the network and communication can still be an issue. (4 MPI processes)

Shared Cluster – In our shared cluster, multiple jobs run on each node, with some jobs being more demanding than others. Even though each node has 16 cores, the amount of resources each job can utilize changes over time. In terms of communication, we have a GigaBit network. (11 MPI processes, running on 11 nodes)

We run a simple code on the cloud/cluster, with MPI communication. We generate two matrices $A,B\in R^{n\times n}$, then synchronize all MPI processes and compute $C=A\cdot B$ using the GSL C-BLAS library. The time is measured and recorded as computational time. After the matrix product is computed, the result is sent to the master/root node using asynchronous communication, and the time required is recorded. The process is repeated 3000 times.

The results of the experiment described above are captured in Figure 14. As expected, on the Amazon EC2 cloud, the matrix-matrix multiplication takes roughly the same time for all replications and the noise in communication is relatively small. In this example the cost of communication is negligible when compared to the cost of computation. On our shared cluster, one cannot guarantee that all resources are exclusively used for a specific process, and thus, the computation and communication time is considerably more stochastic and unbalanced. For some cases the difference between the minimum and maximum computation (communication) time varies by an order of magnitude or more. Hence, on such a platform a fault-tolerant algorithm that only uses information from nodes that return an update within a preallocated budget is a natural choice.

In Figures 15-19 we show a comparison of the proposed robust multi-batch L-BFGS method and the multi-batch L-BFGS method that does not enforce sample consistency (L-BFGS). In these experiments, $p$ denotes the probability that a single node (MPI process) will not return a gradient evaluated on local data within a given time budget. We illustrate the performance of the methods for $\alpha=0.1$ and $p\in\{0.1,0.2,0.3,0.4,0.5\}$. We observe that the robust implementation is not affected much by the failure probability $p$.

Appendix C Scaling of Robust Multi-Batch L-BFGS Implementation

In this Section, we study the strong and weak scaling properties of the robust multi-batch L-BFGS method on an artificial dataset. For various values of $r$ and $K$, we measure the time needed to compute a gradient (Gradient) and the time needed to compute and communicate the gradient (Gradient+C), as well as, the time needed to compute the L-BFGS direction (L-BFGS) and the associated communication overhead (L-BFGS+C).

Figure 20 depicts the strong scaling properties of our proposed algorithm. We generate a dataset with $n=10^{7}$ samples and $d=10^{4}$ dimensions, where each sample has 160 randomly chosen non-zero elements (dataset size 24GB). We run our code for different values of $r$ (different batch sizes $S_{k}$), with $K=1,2,\dots,128$ number of MPI processes.

One can observe that the compute time for the gradient and the L-BFGS direction decreases as $K$ is increased. However, when communication time is considered, the combined cost increases slightly as $K$ is increased. Notice that for large $K$, even when $r=10\%$ (i.e., $10\%$ of all samples processed in one iteration, $\sim$18MB of data), the amount of local work is not sufficient to overcome the communication cost.

C.2 Weak Scaling - Fixed Problem Dimension, Increasing Data Size

In order to illustrate the weak scaling properties of the algorithm, we generate a data-matrix $X\in R^{10^{7}\times 10^{4}}$, and run it on a shared cluster with $K=1,2,4,8,\dots,128$ MPI processes. For a given number of MPI processes ($K$), each sample contains $10\cdot K$ non-zero elements. Effectively, the dimension of the problem is fixed, but sparsity of the data is decreased as more MPI processes are used. The size of the input data is 1.5 $\cdot K$ GB (i.e., 1.5GB per MPI process).

The compute time for the gradient is almost constant, this is because the amount of work per MPI process (rank) is almost identical; see Figure 21. On the other hand, because we are using a Vector-Free L-BFGS implementation for computing the L-BFGS direction, the amount of time needed for each node to compute the L-BFGS direction is decreasing as $K$ is increased. However, increasing $K$ does lead to larger communication overhead, which can be observed in Figure 21. For $K=128$ (192GB of data) and $r=10\%$, almost 20GB of data are processed per iteration in less than 0.1 seconds, which implies that one epoch would take around 1 second.

C.3 Increasing Problem Dimension, Fixed Data Size and K𝐾K

In this experiment, we investigate the effect of a change in the dimension $d$ of the problem on the performance of the algorithm. We fix the size of data ($29GB$) and the number of MPI processes ($K=8$). We generate data with $n=10^{7}$ samples, where each sample has 200 non-zero elements. Figure 22 shows that increasing the dimension $d$ has a mild effect on the computation time of the gradient, while the effect on the time needed to compute the L-BFGS direction is more apparent. However, if communication time is taken into consideration, the time required for the gradient computation and the L-BFGS direction computation increase as $d$ is increased.