A Robust Multi-Batch L-BFGS Method for Machine Learning
Albert S. Berahas Martin Takáč
math.OC cs.LG stat.ML












A Robust Multi-Batch L-BFGS Method for Machine Learning[This work substantially extends <cit.> published at the Neural Information Processing Systems (NeurIPS) conference in 2016.]
    
Albert S. Berahasa^∗^∗Corresponding author. Email: <albertberahas@u.northwestern.edu>
and Martin Takáča
aDepartment of Industrial and Systems Engineering, Lehigh University, Bethlehem, PA, USA
v5.0 released July 2015

    December 30, 2023
==============================================================================================================================================================================================================================




This paper describes an implementation of the L-BFGS method designed to deal with two adversarial situations. The first occurs in distributed computing environments where some of the computational nodes devoted to the evaluation of the function and gradient are unable to return results on time. A similar challenge occurs in a multi-batch approach in which the data points used to compute function and gradients are purposely changed at each iteration to accelerate the learning process. Difficulties arise because L-BFGS employs gradient differences to update the Hessian approximations, and when these gradients are computed using different data points the updating process can be unstable. This paper shows how to perform stable quasi-Newton updating in the multi-batch setting, studies the convergence properties for both convex and nonconvex functions, and illustrates the behavior of the algorithm in a distributed computing platform on binary classification logistic regression and neural network training problems that arise in machine learning.



L-BFGS, multi-batch, fault-tolerant, sampling, consistency, overlap.


90C30; 90C06; 90C53. 




§ INTRODUCTION




It is common in machine learning to encounter optimization problems involving tens of millions of training examples and millions of variables.  To deal with the demands of time, storage and processing power imposed by such applications, high performance implementations of stochastic gradient and batch quasi-Newton methods have been developed; see e.g.,  <cit.>. In this paper we study a batch approach based on the L-BFGS method <cit.> that strives to reach the right balance between efficient learning and productive parallelism.

At present, due to its fast learning properties and low per-iteration cost, the preferred method for very large scale applications is the stochastic gradient  (SG) method <cit.>, and its variance-reduced and accelerated variants <cit.>. These methods are implemented either in an asynchronous manner (e.g., using a parameter server in a distributed setting) or following a synchronous mini-batch approach that exploits parallelism in the gradient evaluations <cit.>. 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 <cit.>. 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 SG 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 <cit.>.

An alternative to SG-type methods are batch methods, such as L-BFGS <cit.>, because they parallelize well and are able to achieve high training accuracy. Batch methods allow for more computation per node, so as to achieve a better balance with the communication costs <cit.>; however, batch methods are not as efficient learning algorithms as SG methods in a sequential setting <cit.>. To benefit from both types of methods, some high performance machine learning systems implement both types of methods  <cit.>, and algorithms that transition from the stochastic to the batch regime <cit.> have also received attention recently.



The goal of this paper is to propose a single method that selects a sizeable subset (batch) of the training data to compute a step and changes this batch at every iteration to improve the learning abilities of the method. In order to differentiate it from the mini-batch approach used in conjunction with the SG method, which employs a very small subset of the training data, we call this the multi-batch approach. In this regime 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. However, the multi-batch approach can cause difficulties to quasi-Newton methods as these methods employ gradient differences to update the Hessian approximations.



More specifically, in this paper we study how to design a robust multi-batch implementation of the limited-memory version of the classical BFGS method <cit.>—which we call the multi-batch L-BFGS method—in the presence of two adverse situations <cit.>. The first occurs in parallel implementations when some of the computational nodes devoted to the evaluation of the function and gradient are unable to return results on time, i.e., in the presence of faults.  This amounts to using different data points to evaluate the function and gradient at the beginning and the end of the iteration, which can be  harmful to quasi-Newton methods since they employ gradient differences  to update Hessian approximations.  A similar challenge occurs in a multi-batch approach in which the data points used to compute the function and gradient are purposely changed at each iteration (or every several iterations) to  accelerate the learning process. The main objective of this paper is to show that stable quasi-Newton updating can be achieved in these 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 insignificant, something that can be expected, or easily enforced, in most situations.

Recently, several stochastic quasi-Newton (SQN) methods have been proposed; see e.g., <cit.>. The methods enumerated above differ in three major aspects: (i) the update rules for the curvature (correction) pairs and the Hessian approximation, (ii) the frequency of updating, and (iii) the required extra computational cost and synchronization required. Our method is different from these methods predominantly due to the fact that it does not modify the BFGS update equations or the form of the curvature pairs, and does not require extra (gradient) computations.  Additionally, our method is designed to work in a distributed settings with faults, in which faults occur randomly and sample consistency cannot be assumed, and as such several SQN methods are not suitable.


We analyze the convergence properties of the multi-batch L-BFGS method using a fixed step length strategy, as well as a diminishing step length strategy, on both strongly convex and nonconvex problems. This is appropriate in our setting, as using a fixed step length approach is popular in practice, and facilitates the study of the stability of quasi-Newton updating in a distributed setting. For strongly convex functions, we show that the algorithm converges, at a linear rate, to an approximate solution whose accuracy depends on the variance of the gradients and the step length. In the nonconvex setting, we show that if cautious BFGS updating is employed, the expected value of the average norm-squared of the gradient is bounded. 

We present numerical experiments on a plethora of problems that arise in machine learning and deep learning. We first illustrate the robustness of our proposed approach on binary classification logistic regression problems on a distributed computing platform with faults and in the serial multi-batch setting. The results indicate that the proposed method achieves a good balance between computation and communication costs. Moreover, we present results on neural network training tasks that illustrate that when larger batch-size is used, our algorithm is competitive with the state-of-the-art. Finally, we demonstrate the strong and weak scaling properties of the proposed method.




The paper is organized as follows. In Section <ref> we describe the multi-batch L-BFGS method in detail. In Section <ref> we provide convergence analyses for the proposed method for strongly convex and nonconvex functions. Numerical results that illustrate the practical performance and robustness of the multi-batch L-BFGS method are reported in Section <ref>. Finally, in Section <ref> we provide some concluding remarks.




§ A MULTI-BATCH QUASI-NEWTON METHOD




Ideally, in supervised learning, one seeks to minimize  expected risk, defined as

    R(w)= ∫_Ω f(w; x, y) d P(x, y) = 𝔼[f(w;x,y)],

 where (x,y) are input-output pairs, f : ℝ^d →ℝ is the composition of a prediction function (parametrized by w) and a loss function, and Ω is the space of input-output pairs endowed with a probability distribution P(x,y). Since the distribution P is typically not known, one  approximates (<ref>) by the empirical risk
 
    F(w) = 1/n∑_i=1^nf(w;x^i,y^i) 1/n∑_i=1^nf_i(w),

where (x^i, y^i), for i=1, …, n, denote the training examples, also referred to as data points or samples. The training problem consists of finding an optimal choice of the parameters w ∈ℝ^d with respect to F, i.e., to compute a solution of the problem

    min_w∈ℝ^dF(w) =1/n∑_i=1^nf_i(w).




In a pure batch approach, one applies a gradient-based method to the deterministic optimization problem (<ref>). In this regime, a popular method is L-BFGS <cit.>. When n is large, it is natural to parallelize the computation of F and ∇ F by assigning the evaluation of component functions f_i, or subsets of the component functions, to different processors. If this is done on a distributed computing platform, it is possible for some of the computational nodes, dedicated to a portion of the evaluation of the objective function and the gradient, to be slower than the rest. In this case, the contribution of the slow (or unresponsive) computational nodes could potentially be ignored given the stochastic nature of the true objective function (<ref>). However, this leads 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, as mentioned above. Hence, we seek to develop a fault-tolerant version 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 only a subset of the data is used to compute the gradient at every iteration. We consider a method in which the dataset is randomly divided into  a number of batches and the minimization is performed with respect to a different batch at every iteration. Specifically, at the k-th iteration the algorithm chooses S_k ⊂{1, …, n}, computes 

    F^S_k(w_k)=1/|S_k|∑_i∈ S_kf_i(w_k),     g_k^S_k=∇ F^S_k(w_k) = 1/|S_k|∑_i∈ S_k∇ f_i(w_k) ,

and takes a step along the direction - H_k g_k^S_k, where H_k is an approximation to ∇^2 F(w_k)^-1. Allowing the sample S_k to change freely at every iteration gives this approach flexibility and is beneficial to the learning process. Note, 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 regime are similar in nature, i.e., the samples S_k used change from one iteration to the next. The main difference is that node failures create unpredictable changes to the samples, whereas a multi-batch method has control over the 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 SG method, which employs small mini-batches. Our algorithm operates with much larger batches so that distributing the computation of the function and gradient is beneficial, and the compute time 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.

The robust implementation of the L-BFGS method, proposed in <cit.>, is based on the following observation:  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 have an overlap, so that

    O_k= S_k∩ S_k+1≠∅.

 One can then perform stable quasi-Newton updating by computing gradient differences based on this overlap, i.e.,  by defining 

    y_k+1=g_k+1^O_k-g_k^O_k,      s_k+1 = w_k+1-w_k,

in the notation given in (<ref>), and using this correction pair (y_k, s_k) 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 along the most recent displacement, and leads 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 <cit.>. More specifically, 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 used to define the search direction - H_k g_k^S_k. In summary, by ensuring that unresponsive nodes do not constitute the vast majority of all compute nodes in a fault-tolerant parallel implementation, or by exerting a small degree of control in the creation of the samples S_k in the multi-batch regime, one can design a robust method that naturally builds upon the fundamental properties of BFGS updating.

We should mention that a  commonly used fix  for ensuring stability of quasi-Newton updating in machine learning is to enforce gradient consistency <cit.>, i.e., to use the same sample S_k to compute gradient evaluations at the beginning and the end of the iteration, at the cost of double gradient evaluations.  Another popular remedy is to use the same batch S_k for multiple iterations <cit.>, alleviating the gradient inconsistency problem at the price of slower convergence. In this paper, we assume that such sample consistency is not possible (the fault-tolerant case) or desirable (the multi-batch regime), and wish to design and analyze an implementation of L-BFGS that imposes minimal restrictions in the changes of the sample.



 §.§ Specification of the Method

Let us begin by considering a robust implementation of the multi-batch BFGS method and  then consider its limited memory version. At  the k-th iteration, the multi-batch BFGS algorithm chooses a set S_k ⊂{1, …, n} and computes a new iterate by the formula

    w_k+1=w_k-α_kH_k g_k^S_k ,

where α_k is the step length, g_k^S_k is the batch gradient (<ref>) and H_k is the inverse BFGS
Hessian matrix approximation that is updated at every iteration by means
of the formula

    H_k+1=V_k^TH_kV_k+ρ_ks_ks_k^T,     ρ_k=1/y_k^Ts_k,      V_k=1-ρ_ky_ks_k^T .

 To compute the correction vectors (s_k, y_k), we determine the overlap set O_k = S_k ∩ S_k+1 consisting of the samples that are common at the k-th and k+1-st iterations. We define  
 
    g_k^O_k=∇ F^O_k(w_k)=1/|O_k|∑_i∈ O_k∇ f_i(w_k),

and compute the correction pairs as in (<ref>). This completely specifies the algorithm, except for the choice of step length α_k; in this paper we consider constant and diminishing step lengths.

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 search direction (matrix-vector product) in (<ref>) it is not necessary to form the dense matrix H_k since one can obtain this product via the two-loop recursion <cit.>, using the m most recent correction pairs. Employing this mechanism, the search direction can be computed in 𝒪(d) floating operations, where d is the number of variables. 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 multi-batch limited-memory BFGS algorithm is given in Algorithm <ref>, and depends on several parameters.  The parameter r denotes the fraction of samples in the dataset used to define the gradient, i.e., r = | S|/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 = | O|/| S|.





 §.§ Sample Generation



The fault-tolerant and multi-batch settings differ in the way the samples S_k and O_k are formed (Lines <ref> & <ref>, Algorithm <ref>). In the former, sampling is done automatically as a by-product of the nodes that fail to return a computation (gradient evaluation). In the latter, the samples S_k and O_k used at every iteration are purposefully changed in order to accelerate the learning process, thus sampling is user controlled. In either setting, independent sampling can be achieved; a necessary condition to establish convergence results. We first describe the fault-tolerant setting, and then propose two sampling strategies that can be employed in the multi-batch setting. Let T= { (x^i, y^i), for i=1, …, n} denote the training set.











  
*Fault-Tolerant
Consider a distributed implementation in which slave nodes read the current iterate w_k from the master node, compute a local gradient on a subset of the dataset, and send it back to the master node for aggregation in the calculation (<ref>).  Given a time (computational) budget, it is possible for some nodes to fail to return a result. The schematic in Figure <ref> illustrates the gradient calculation across two iterations, k and k+1, in the presence of faults.  Here,  B is the total number of slave nodes, ℬ_i for i=1,...,B denote the batches of data that each slave node i receives (T = ∪_i ℬ_i),  and 
∇̃f(w) is the gradient calculation using all nodes that responded within the preallocated time.

Let 𝒥_k⊂{1,2,...,B} and 𝒥_k+1⊂{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 = ∪_j∈𝒥_kℬ_j and S_k+1 = ∪_j∈𝒥_k+1ℬ_j, and we define O_k = ∪_j ∈𝒥_k∩𝒥_k+1ℬ_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 are independent if node failures occur randomly. On the other hand, if the same set of nodes fail, then the sample creation will be biased, which is harmful both in theory and in practice.  One way to ensure independent sampling is to shuffle and redistribute the data to all nodes after every iteration or after a certain number of iterations.



  
*Multi-Batch Sampling
In the multi-batch setting several strategies can be employed, with the only restriction that consecutive batches S_k and S_k+1 should, to a certain degree, overlap. We propose two sampling strategies: (i) overlaps O_k are forced in the sample creation process, (ii) the overlapping set O_k is subsampled from the batch S_k. In practice the two strategies perform on par, however, there is a subtle difference. In the second strategy the batches are sampled independently, something that is not true for the strategy in which overlapping samples are forced. The independent sampling strategy of course does not come for free as this strategy incurs an increase in computational cost per iteration. However, as mentioned above, the overlapping set O_k need not be very large, and thus the increase in cost is negligible as compared to the rest of the computation. We now describe the two approaches in more detail.

Figure <ref> 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 epoch (pass over the whole dataset) 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 | S_k | elements from {1,… n}. The set O_k is then formed by randomly selecting | O_k | elements from S_k (subsampling). Note, in this sampling strategy the samples O_k need not be in the set S_k+1. 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 <ref>) when applied to the minimization of strongly convex  and nonconvex objective functions, using a fixed step length strategy, as well as a diminishing step length strategy. We assume that the goal is to minimize the empirical risk F (<ref>), but note that a similar analysis could be used to study the minimization of the expected risk (<ref>).  



 §.§ Strongly Convex case

Due to the stochastic nature of the multi-batch approach, every iteration of Algorithm <ref> 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^⋆, but only convergence to a neighborhood of w^⋆ <cit.>. 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.




  
*Assumptions A
 
 

  * F is twice continuously differentiable.

  * There exist positive constants λ̂ and Λ̂ such that
λ̂ I ≼∇^2F^O(w) ≼Λ̂ I
for all w ∈ℝ^d and all sets O ⊂{1,2,…,n} of length |O| = o· r· n.

  * There exist constants γ≥ 0 and η≥ 1 such that 𝔼_S[ ∇  F^S(w) ^2 ] ≤γ^2 + η∇ F(w)^2
for all w ∈ℝ^d and all sets S⊂{1,2,…,n}  of length |S|=r · n.

  * The samples S are drawn independently and ∇ F^S(w) is an unbiased estimator of the true gradient ∇ F(w) for all w ∈ℝ^d, i.e.,
𝔼_S[ ∇ F^S(w)] = ∇ F(w).

Note that Assumption A.2 implies that the entire Hessian ∇^2F(w) also satisfies 

    λ I ≼∇^2F(w) ≼Λ I,    ∀ w ∈ℝ^d,

for some constants λ,  Λ>0. Assuming that every subsampled 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 <cit.>.

	
If Assumptions A.1 & A.2 hold, there exist constants 0<μ_1≤μ_2 such that the inverse Hessian approximations {H_k} generated by Algorithm <ref> satisfy

    μ_1 I ≼ H_k ≼μ_2 I,    for  k=0,1,2,…





Instead of analyzing the inverse Hessian approximation H_k, we study the Hessian approximation B_k = H_k^-1. In this case, the limited memory quasi-Newton updating formula is given as follows:


  * Set B_k^(0)=y_k^Ty_k/s_k^Ty_kI and m̃ = min{k,m}; where m is the memory in L-BFGS.

  * For i=0,...,m̃-1 set j=k-m̃+1+i and compute

    B_k^(i+1)=B_k^(i)-B_k^(i)s_js_j^TB_k^(i)/s_j^TB_k^(i)s_j + y_jy_j^T/y_j^Ts_j.


  * Set B_k+1 = B_k^(m̃).


The curvature pairs s_k and y_k are updated via the following formulae

    y_k+1=g_k+1^O_k-g_k^O_k,      s_k+1 = w_k+1-w_k.

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

    y_k^Ts_k ≥1/Λ̂y_k^2       ⇒  y_k^2 /y_k^Ts_k≤Λ̂.

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

    y_k^Ts_k ≤1/λ̂y_k^2       ⇒  y_k^2 /y_k^Ts_k≥λ̂.
   







Combining the upper and lower bounds (<ref>) and (<ref>)  

    λ̂≤y_k^2/y_k^Ts_k≤Λ̂.


The above proves that the eigenvalues of the matrices B_k^(0)=y_k^Ty_k/s_k^Ty_kI 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. Let tr(B) and (B) denote the trace and determinant of matrix B, respectively, and set j_i = k-m̃+i. The trace of the matrix B_k+1 can be expressed as, 

    tr(B_k+1)    = tr(B_k^(0)) - tr∑_i=1^m̃(B_k^(i-1)s_j_is_j_i^TB_k^(i-1)/s_j_i^TB_k^(i-1)s_j_i) + tr∑_i=1^m̃y_j_iy_j_i^T/y_j_i^Ts_j_i
       ≤ tr(B_k^(0)) + ∑_i=1^m̃y_j_i^2/y_j_i^Ts_j_i
    ≤ tr(B_k^(0)) + m̃Λ̂
     ≤ C_1,

for some positive constant C_1, where the inequalities above are due to (<ref>), 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 <cit.>, the determinant of the matrix  B_k+1 generated by the multi-batch L-BFGS method can be expressed as, 

    (B_k+1)    = (B_k^(0)) ∏_i=1^m̃y_j_i^Ts_j_i/s_j_i^TB_k^(i-1)s_j_i
    =  (B_k^(0)) ∏_i=1^m̃y_j_i^Ts_j_i/s_j_i^Ts_j_is_j_i^Ts_j_i/s_j_i^TB_k^(i-1)s_j_i
       ≥(B_k^(0)) ( λ̂/C_1)^m̃
     ≥ C_2,

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 (<ref>) and determinant (<ref>) inequalities derived above imply that the largest eigenvalues of all matrices B_k are bounded from above, uniformly, and that the smallest eigenvalues of all matrices B_k are bounded away from zero, uniformly. 


Before we present the main theorem for the multi-batch L-BFGS method that employs constant step lengths, we state one more intermediate Lemma that bounds the distance between the function value at any point w∈ℝ^d and the optimal function value with respect to the norm of the gradient squared.

 	
Let Assumptions A.1 & A.2 hold, and let F^⋆ = F(w^⋆), where w^⋆ is the minimizer of F. Then, for all w∈ℝ^d,

    2λ(F(w)-F^⋆) ≤∇ F(w)^2.



 As a result of Assumptions A.1, A.2 and (<ref>), for all x,y ∈ℝ^d

    F(x) ≤ F(y) + ∇ F(y)^T (x-y) + 1/2λ∇ F(y) - ∇ F(x)^2;

see <cit.>.
Let x=w and y=w^⋆

    F(w)    ≤ F^⋆ + ∇ F(w^⋆)^T (w-w^⋆) + 1/2λ∇ F(w) - ∇ F(w^⋆)^2 
       ≤ F^⋆ + 1/2λ∇ F(w) ^2.

Re-arranging the above expression yields the desired result.


Utilizing Lemmas <ref> and <ref>, we show that the multi-batch L-BFGS method with a constant step length converges linearly to a neighborhood of the optimal solution.




Suppose that Assumptions A.1-A.4 hold, and let F^⋆ = F(w^⋆), where w^⋆ is the minimizer of F. Let {w_k} be the iterates generated by Algorithm <ref>, where α_k = α satisfies

    0 <  α≤μ_1/μ_2^2 ηλ,

and w_0 is the starting point.
Then for all k≥ 0,

    𝔼[ F(w_k) - F^⋆]    ≤( 1-αμ_1 λ)^k [ F(w_0) - F^⋆] + [ 1-(1-αμ_1 λ)^k]αμ_2^2 γ ^2 Λ/2 μ_1 λ
       αμ_2^2 γ ^2 Λ/2 μ_1 λ.




We have that

    F(w_k+1)     = F(w_k -α H_k ∇ F^S_k(w_k)) 
       ≤ F(w_k) + ∇ F(w_k)^T (-α H_k ∇ F^S_k(w_k)) + Λ/2α H_k ∇ F^S_k(w_k)^2 
       ≤ F(w_k) - α∇ F(w_k)^T  H_k ∇ F^S_k(w_k) + α^2 μ_2^2 Λ/2∇ F^S_k(w_k)^2,

where the first inequality arises due to (<ref>), and the second inequality arises as a consequence of Lemma <ref>. Taking the expectation (over S_k) of equation (<ref>)

    𝔼_S_k[ F(w_k+1)]    ≤ F(w_k) - α∇ F(w_k)^T  H_k ∇ F(w_k) + α^2 μ_2^2 Λ/2𝔼_S_k[ ∇ F^S_k(w_k)^2 ] 
       ≤ F(w_k) - αμ_1 ∇ F(w_k) ^2  + α^2 μ_2^2 Λ/2( γ^2 + η∇ F(w)^2) 
        = F(w_k) - α( μ_1 - αμ_2^2 ηΛ/2)∇ F(w_k) ^2  + α^2 μ_2^2 γ^2Λ/2
       ≤ F(w_k) - αμ_1/2∇ F(w_k) ^2  + α^2 μ_2^2 γ^2Λ/2

where the first inequality makes use of Assumption A.4, the second inequality arises due to Lemma <ref> and Assumption A.3 and the third inequality is due to the step length (<ref>). Since F is λ-strongly convex, we can substitute the result of Lemma <ref> in (<ref>),

    𝔼_S_k [F(w_k+1)]    ≤ F(w_k) - αμ_1/2∇ F(w_k) ^2  + α^2 μ_2^2 γ^2 Λ/2
       ≤ F(w_k) - αμ_1 λ  [F(w_k) - F^⋆]  + α^2 μ_2^2 γ^2 Λ/2.

Let 

    ϕ_k = 𝔼[F(w_k) - F^⋆],

where the expectation is over all batches S_0,S_1,...,S_k-1 and all history starting with w_0. Equation (<ref>) can be expressed as,

    ϕ_k+1≤ (1 - αμ_1 λ ) ϕ_k + α^2 μ_2^2 γ^2 Λ/2.

Since the step length is chosen according to (<ref>) we deduce that 0 ≤ (1 - αμ_1 λ ) < 1.

Subtracting αμ_2^2γ^2Λ/2μ_1 λ from either side of (<ref>) yields

    ϕ_k+1 - αμ_2^2γ^2Λ/2μ_1 λ   ≤ (1 - αμ_1 λ ) ϕ_k + α^2 μ_2^2 γ^2 Λ/2 - αμ_2^2γ^2Λ/2μ_1 λ
        = (1 - αμ_1 λ ) [ ϕ_k  - αμ_2^2γ^2Λ/2μ_1 λ].

Recursive application of (<ref>) yields

    ϕ_k - αμ_2^2γ^2Λ/2μ_1 λ   ≤  (1 - αμ_1 λ )^k [ ϕ_0  - αμ_2^2γ^2Λ/2μ_1 λ],

and thus, 

    ϕ_k≤ (1 - αμ_1 λ )^k ϕ_0  + [ 1-(1-αμ_1 λ)^k]αμ_2^2 γ ^2 Λ/2 μ_1 λ.

Finally using the definition of ϕ_k (<ref>) with the above expression yields the desired result

    𝔼[ F(w_k) - F^⋆] ≤( 1-αμ_1 λ)^k [ F(w_0) - F^⋆] + [ 1-(1-αμ_1 λ)^k]αμ_2^2 γ ^2 Λ/2 μ_1 λ.




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, 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 consider these tradeoffs further in Section <ref>, where we also note that larger batch sizes increase the opportunities for parallelism and improve the limiting accuracy in the solution, but slow down the learning abilities of the algorithm. We should also mention that unlike the first-order variant of the algorithm (H_k = I), the step length range prescribed by the multi-batch L-BFGS method depends on μ_1 and μ_2, the smallest and largest eigenvalues of the L-BFGS Hessian approximation. In the worst-case, the presence of the matrix H_k can make the limit in Theorem <ref> 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.

 One can establish convergence of the multi-batch L-BFGS method to the optimal solution w^⋆ by employing a  sequence of step lengths {α_k } that converge to zero according to the schedule  proposed by <cit.>. However, that provides only a sub-linear 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 <ref> is more relevant to practice; nonetheless, we state the theorem here for completeness, and, for brevity, refer the reader to <cit.> for more details and the proof.
 
	
Suppose that Assumptions A.1-A.4 hold, and let F^⋆ = F(w^⋆), where w^⋆ is the minimizer of F. Let {w_k} be the iterates generated by Algorithm <ref> with 

    α_k = β/k+1    and     β > μ_1/μ_2^2 ηλ,

starting from w_0. Then for all k≥ 0,

    𝔼[ F(w_k) - F^⋆] ≤Q(β)/k+1,

where Q(β) = max{μ_2^2β^2γ^2Λ/2(2μ_1λβ-1), F(w_0) - F^⋆}.
 

Theorem <ref> shows that, for strongly convex functions, the multi-batch L-BFGS method with an appropriate schedule of diminishing step lengths converges to the optimal solution at a sub-linear rate. We should mention that another way to establish convergence to the optimal solution for the multi-batch L-BFGS method is to employ variance reduced gradients <cit.>. In this setting, one can establish linear convergence to the optimal solution using constant step lengths. We defer the analysis of the multi-batch L-BFGS method that employs variance reduced gradients to a different study <cit.>. 



 §.§ Nonconvex case

The BFGS method is known to fail on nonconvex problems <cit.>. 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 (L-)BFGS method in the nonconvex setting several techniques have been proposed including cautious updating <cit.>, modified updating <cit.> and damping <cit.>. Here we employ a cautious strategy that is well suited to our particular algorithm; we skip the Hessian update, i.e., set H_k+1 = H_k,  if the curvature condition  

    y_k^Ts_k ≥ϵ s_k ^2

 is not satisfied, where ϵ>0 is a predetermined constant. On the other hand, sufficient curvature is guaranteed when the updates are not skipped. 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 <ref>). The analysis presented in this section is based on the following assumptions.
 
  

  
*Assumptions B
 



  * F is twice continuously differentiable.

  * The gradients of F are Λ-Lipschitz continuous for all w ∈ℝ^d; the gradients of F^S are Λ_S-Lipschitz continuous for all w ∈ℝ^d and all sets S ⊂{1,2,…,n} of length |S|=r · n; and, the gradients of F^O are Λ_O-Lipschitz continuous for all w ∈ℝ^d and all sets O ⊂{1,2,…,n} of length |O| = o· r· n.

  * The function F(w) is bounded below by a scalar F .

  * There exist constants γ≥ 0 and η≥ 1 such that 

𝔼_S[ ∇  F^S(w) ^2 ] ≤γ^2 + η∇ F(w)^2


for all w ∈ℝ^d and all sets S⊂{1,2,…,n}  of length |S|=r · n. 

  * The samples S are drawn independently and ∇ F^S(w) is an unbiased estimator of the true gradient ∇ F(w) for all w ∈ℝ^d, i.e.,

𝔼_S [ ∇ F^S(w) ] = ∇ F(w).




Similar to the strongly convex case, we first show that the eigenvalues of the L-BFGS Hessian approximations are bounded above and away from zero.

		
Suppose that Assumptions B.1 & B.2  hold. Let {H_k } be the inverse Hessian approximations generated by Algorithm <ref>, with the modification that the inverse Hessian approximation update is performed only when  (<ref>) is satisfied,



for some ϵ >0, else H_k+1 = H_k. Then, there exist constants 0<μ_1≤μ_2 such that 

    μ_1 I ≼ H_k ≼μ_2 I,    for  k=0,1,2,…




Similar to the proof of Lemma <ref>, 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

    y_k+1=g_k+1^O_k-g_k^O_k,      s_k+1 = w_k+1-w_k.

The skipping mechanism (<ref>) provides both an upper and lower bound on the quantity y_k^2 /y_k^Ts_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 (<ref>). We have from (<ref>) that

    ϵ s_k ^2    ≤ y_k^Ts_k ≤ y_k  s_k   ⇒   s_k ≤1/ϵ y_k ,
 
from which it follows that

    s_k^Ty_k ≤ s_k  y_k ≤1/ϵ y_k ^2   ⇒   y_k ^2/s_k^Ty_k≥ϵ.

The upper bound is attained by the Lipschitz continuity of sample gradients (Assumption B.2),

    y_k^Ts_k  ≥ϵ s_k ^2  
    	≥ϵ y_k ^2/Λ_O^2  ⇒   y_k ^2/s_k^Ty_k≤Λ_O^2/ϵ.
 
Combining (<ref>) and (<ref>),

    ϵ≤y_k^2 /y_k^Ts_k≤Λ_O^2/ϵ.


The above proves that the eigenvalues of the matrices B_k^(0)=y_k^Ty_k/s_k^Ty_kI 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 <ref>, the only difference being that the last inequality in (<ref>) comes as a result of the cautious update strategy.


We now follow the analysis in <cit.> 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 hold. Let {w_k} be the iterates generated by Algorithm <ref>, with the modification that the inverse Hessian approximation update is performed only when  (<ref>) is satisfied,



for some ϵ >0, else H_k+1 = H_k,  where α_k = α satisfies

    0 <  α≤μ_1/μ_2^2 ηλ,

and w_0 is the starting point. Then, for all k≥0,

    𝔼[1/τ∑_k=0^τ-1∇ F(w_k) ^2 ]    ≤αμ_2^2 γ^2 Λ/μ_1  + 2[ F(w_0) - F]/αμ_1 τ
       αμ_2^2 γ^2 Λ/μ_1 .




Starting with (<ref>) and taking 







the expectation over all  batches S_0,S_1,...,S_k-1 and all history starting with w_0 yields

    ϕ_k+1-ϕ_k ≤ - αμ_1/2𝔼∇ F(w_k) ^2  + α^2 μ_2^2 γ^2 Λ/2,

where ϕ_k = 𝔼[F(w_k)]. Summing (<ref>) over the first τ iterations

    ∑_k=0^τ-1 [ϕ_k+1-ϕ_k]    ≤ - αμ_1/2∑_k=0^τ-1𝔼∇ F(w_k) ^2 + ∑_k=0^τ-1α^2 μ_2^2 γ^2 Λ/2
       = - αμ_1/2𝔼[∑_k=0^τ-1∇ F(w_k) ^2 ] +  α^2 μ_2^2 γ^2 Λτ/2.

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

    ∑_k=0^τ-1[ϕ_k+1-ϕ_k ]     = ϕ_τ-ϕ_0 
    	= 𝔼[F(w_τ)] -F(w_0) 
    	≥F -F(w_0).

Substituting the above expression into (<ref>) and rearranging terms

    𝔼[∑_k=0^τ-1∇ F(w_k) ^2 ] ≤αμ_2^2 γ^2 Λτ/μ_1  + 2[ F(w_0) - F]/αμ_1 .

Dividing the above equation by τ completes the proof.


This result bounds the average norm of the gradient of F  after the first τ-1 iterations, and shows that, in expectation, the iterates spend increasingly more time in regions where the objective function has a small gradient. Under appropriate conditions, we can establish a convergence rate for the multi-batch L-BFGS method with cautious updates to a stationary point of F, similar to the results proven for the SG method <cit.>. For completeness we state and prove the result.

	 
Suppose that Assumptions B.1-B.5 hold. Let {w_k} be the iterates generated by Algorithm <ref>, with the modification that the inverse Hessian approximation update is performed only when  (<ref>) is satisfied,  



for some ϵ >0, else H_k+1 = H_k. Let

    α_k = α = c/√(τ),    c = √(2(F(w_0) - F̂)/μ_2^2 γ^2 Λ),   δ(α) = μ_1 - αμ_2^2  ηΛ/2,

where 
τ > 
c^2μ_2^4 η^2 Λ^2/4μ_1^2, 
and w_0 is the starting point. Then, 

    min_0 ≤ k ≤τ-1𝔼[∇ F(w_k)^2] ≤√(2(F(w_0) - F̂)μ_2^2 γ^2 Λ/δ(α)^2 τ).



 
Starting with (<ref>), we have

    𝔼_S_k[ F(w_k+1)]    ≤ F(w_k) - α(μ_1 - αμ_2^2  ηΛ/2) ∇ F(w_k) ^2 + α^2 μ_2^2 γ^2Λ/2
        = F(w_k) - αδ(α) ∇ F(w_k) ^2 + α^2 μ_2^2 γ^2Λ/2,

where δ(α) = μ_1 - αμ_2^2  ηΛ/2. We require that this quantity is greater than zero, δ(α)>0; this discussion is deferred to the end of the proof.

Taking an expectation over all  batches S_0,S_1,...,S_k-1 and all history starting with w_0, and rearranging (<ref>) yields

    𝔼[∇ F(w_k) ^2] ≤1/αδ(α)𝔼[F(w_k)-F(w_k+1)] + αμ_2^2 γ^2 Λ/2 δ(α).

Summing over k=0,...,τ-1 and dividing by τ

    min_0 ≤ k ≤τ-1𝔼[∇ F(w_k)^2]    ≤1/τ∑_k=0^τ-1𝔼[∇ F(w_k) ^2] 
       ≤1/αδ(α) τ𝔼[F(w_0)-F(w_τ)] + αμ_2^2 γ^2 Λ/2 δ(α)
       ≤1/αδ(α) τ[F(w_0)-F̂] + αμ_2^2 γ^2 Λ/2 δ(α)
       ≤1/δ(c/√(τ)) c √(τ)[F(w_0)-F̂] + c μ_2^2 γ^2 Λ/2 δ(c/√(τ))√(τ).

The first inequality holds because the minimum value is less than the average value, and the third inequality holds because F̂≤ F(x_τ) (Assumption B.3). The last expression comes as a result of using the definition of the step length, α = c/√(τ). Setting 

    c = √(2(F(w_0) - F̂)/μ_2^2 γ^2 Λ),

yields the desired result.

We now comment on the quantity δ(α) that first appears in (<ref>), and that is required to be positive. To ensure that δ(α)>0, the step length must satisfy, α < 2μ_1/μ_2^2 ηΛ. Since the explicit form of the step length is α = c/√(τ), where c is (<ref>), we require that 

    α =  c/√(τ) < 2μ_1/μ_2^2 ηΛ.

In order to ensure that (<ref>) holds, we impose that

    τ > c^2μ_2^4 η^2 Λ^2/4μ_1^2 = (F(w_0) - F̂)μ_2^2η^2 Λ/2γ^2μ_1^2.



The result of Theorem <ref> establishes a sub-linear rate of convergence, to a stationary point of F, for the multi-batch L-BFGS method on nonconvex objective functions. The result is somewhat strange as it requires a priori knowledge of τ, the total number of iteration. In practice, one would use α_k = 1/√(k), which would result in a 𝒪(1/√(k)) convergence rate.





§ NUMERICAL RESULTS




We present numerical experiments on several problems that arise in machine learning, such as logistic regression binary classification and neural network training, in order to evaluate the performance of the proposed multi-batch L-BFGS method. The experiments verify that the proposed method is robust, competitive and achieves a good balance between computation and communication in the distributed setting. In Section <ref>, we evaluate the performance of the multi-batch L-BFGS method on binary classification tasks in both the multi-batch and fault-tolerant settings. In Section <ref>, we demonstrate the performance of the multi-batch L-BFGS method on neural network training tasks,

 and compare against some of the state-of-the-art methods. Finally, in Section <ref>, we illustrate the strong and weak scaling properties of the multi-batch L-BFGS method.



 §.§ Logistic Regression


In this section, we focus on logistic regression problems; the optimization problem can be stated as:

    min_w ∈ℝ^d F(w) =  1/n∑_i=1^nlog(1+e^-y^i(w^Tx^i))
      + σ/2w^2,

where (x^i, y^i)_i=1^n  denote the training examples and σ = 1/n is the regularization parameter. 

We present numerical results that evaluate the performance of the proposed robust multi-batch L-BFGS scheme  (Algorithm <ref>) in both the multi-batch (Figure <ref>) and fault tolerant (Figure <ref>) settings, on the  dataset[LIBSVM: <https://www.csie.ntu.edu.tw/ cjlin/libsvmtools/datasets/binary.html>]. We compare our proposed method (Robust L-BFGS) 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^S_k+1_k+1 - g^S_k_k; (ii) multi-batch gradient descent (Gradient Descent), which is obtained by setting H_k = I in Algorithm 1; and, (iii) serial SGD (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. In Figures <ref> and <ref> we show the evolution of the norm of the gradient in terms of epochs.



In the multi-batch setting, 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, at the cost of n communications per epoch versus 1/r(1-o) 

communications per epoch. Figure <ref> also indicates that the robust L-BFGS method is not too sensitive to the size of the overlap. Similar behavior was observed on other datasets, in regimes where r· o was not too small; see <cit.>.



Figure  <ref> shows 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 the presence of faults. 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 α=0.1 and p∈{0.1, 0.3, 0.5}. We observe that the robust implementation is not affected much by the failure probability p. Similar behavior was observed on other datasets; see <cit.>.





 

 §.§ Neural Networks


In this section, we study the performance of the multi-batch L-BFGS method[Code available at: <https://github.com/OptMLGroup/Multi-Batch_L-BFGS>.] on Neural Network tasks, on the MNIST and CIFAR10/CIFAR100 datasets[MNIST available at: <http://yann.lecun.com/exdb/mnist/>. CIFAR10/CIFAR100 available at: <https://www.cs.toronto.edu/ kriz/cifar.html>.]. Table <ref> summarizes the network architectures that we used; see <cit.> 
for more details. The first problem is convex; all other problems are nonconvex.








We implemented our algorithm in PyTorch <cit.> and compare against popular and readily available algorithms: (i) SGD <cit.>, and (ii)  Adam <cit.>. We denote our proposed method as LBFGS in the figures in this section. Note, we implemented our method with the cautious updating strategy, and For each method, we conducted a grid search to find the best learning rate α∈{2^0,2^-1,…,  2^-10}, and also investigated the effect of different batch sizes | S | ∈{50, 100, 200, 500, 1000, 2000, 4000 }; see <cit.> 
 for detailed experiments with all batch sizes. For the multi-batch L-BFGS method we also investigated the effect of history length m ∈{1,2,5,10,20 }. The overlap used in our proposed method was 20% of the batch, o = 0.20.



The authors in <cit.> observed that the widely used Barzilai-Borwein-type scaling s_k^T y_k/y_k^T y_k I of the initial Hessian approximation may lead to quasi-Newton updates that are not stable when small batch sizes are employed, especially for deep neural training tasks, and as such propose an Agadrad-like scaling of the initial BFGS matrix. To obviate this instability, we implement a variant of the multi-batch L-BFGS method (LBFGS2) in which we scale the initial Hessian approximations as α I. We ran experiments with both scaling strategies and the overall results were similar. Therefore, in the figures in this section we only show results for the latter strategy. 



Figure <ref> illustrates the evolution of the running maximum of the training accuracy (dashed lines) and testing accuracy (solid lines) for the best parameter settings for each method over the first 20 epochs of training. By best parameter settings we mean the run that achieved the highest training accuracy within 20 epochs. One can make several observations. First, it appears that the multi-batch L-BFGS method is competitive with the first-order methods on all training problems, except for CIFAR10 LeNet, in terms of both training and testing accuracy. Second, for half of the problems (three out of six), the best runs of the multi-batch L-BFGS method use larger batch sizes than the first-order methods. Of course, this benefit is not as clear on the neural network training problems as it is in the logistic regression problems.  Third, the benefits of incorporating second-order information are not as apparent in these nonconvex problems as compared to the problems presented in Section <ref>. We attribute this to two things: (1) nonconvex problems are hard, and (2) quasi-Newton methods are better are capturing curvature information for convex problems.





We now investigate the effect of the batch size. In Figure <ref> we show the evolution of the running maximum of the training/testing accuracy for different batch sizes | S | ∈{50, 500, 4000 } for three of the problems. For a complete set of results, see <cit.>. 
Overall, one can observe that for small batch sizes, the multi-batch L-BFGS variants perform worse than the first order methods. However, when large batches are employed (a regime that is favorable foe GPU computing), the multi-batch L-BFGS method performs on par with the other methods. Moreover, it appears that on several problems the performance of the multi-batch L-BFGS method is less affected by the size of the batch, i.e., the variability in the final training and testing error (after 20 epochs) in terms of batch size is smaller for the multi-batch L-BFGS method than for the stochastic first-order methods; see also <cit.>.




In order to understand why the multi-batch L-BFGS method does not perform well for small batches, we looked at two diagnostic measures: (i) the angle between the true gradient curvature vector y_d and subsampled gradient curvature vector y_s (⟨ y_s,y_d ⟩/ y_s y_d); and (ii) the ratio of subsampled gradient curvature vector to true gradient curvature vector (y_s/y_d). These measures indicate how informative the curvature information captured by the multi-batch L-BFGS method really is. Values close to 1 (dashed red lines) are ideal for both measures. We chose 3 different points (the starting point, a point after 3 epochs of Adam, and a point after 10 epochs of Adam). From those points we took a gradient descent step with sufficiently small step length, and computed the true gradient curvature vector (y_d). We also computed 100 different stochastic variants of the gradient curvature vector (y_s) using different batch sizes (| S | ∈{50, 100, 200, 500, 1000, 2000, 4000 }), and calculated the values of the two metrics. We illustrate the results in Figure <ref>; see <cit.> 
for more results. Several observations can be drawn from this figure. First, not surprisingly, the metrics improve (get closer to 1) as the batch size increases. Second, for the convex case, the metrics perform as expected both close and far from the solution; as a result (sufficiently) good curvature information is captured and the method performs well. On the other hand, for the nonconvex problems, the metrics indicate that, especially for small batch sizes, the curvature information captured can be terrible.




r0.5


  

    
    < g r a p h i c s >

  

Relative slow-down of computational time to compute gradients for different batch-sizes compared to the computational time to compute gradients with batch-size 1.



We should note that using a large batch size is not a bottleneck for current high-performance computing hardware, on the contrary, using small batch sizes leads to under-utilization of computational hardware and in fact hinders the ability to parallelize the methods. Figure <ref> illustrates this phenomenon; we show the average computational time of a gradient evaluation (over 1,000 evaluations) for different batch sizes ranging from 1 to 4096 relative to the computational time of the computation of the gradient using a single data point. The computation was performed on an NVIDIA Tesla K80 GPU using PyTorch. It is clear that for the Multiclass Linear classification task, the compute time of a single gradient is roughly the same as the compute time of a gradient based on a batch of size 1024, whereas for the larger training tasks, the compute time of the gradients appear constant up to batch sizes of roughly 128. We should  note however that their is a risk of decreased generalization when increasing the batch size, unless other strategies such as modifying the step size or regularization are used; see e.g., <cit.>.





































































 








































































































































 §.§ Scaling of the 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 artificial data. For various values of batch size (r) and nodes (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). The function of which we are computing the gradient is logistic regression. The L-BFGS direction is computed using the Vector-Free L-BFGS implementation <cit.>. We should note that the time to compute the gradient, which of course is required for computing the L-BFGS direction, is not included in L-BFGS and L-BFGS+C. We report the extra time to compute the L-BFGS step, after having computed the gradient. Thus, the goal of this section is to show that the time needed to compute the L-BFGS direction is insignificant compared to the cost of computing the gradient, which is needed in any case to run first-order methods. 





  §.§.§ Strong Scaling


r0.4


  

    
    < g r a p h i c s >

  

  

Strong scaling of robust multi-batch L-BFGS on a problem with artificial data; 
  n=10^7 and d=10^4. Each sample has 160 non-zero elements (dataset size 24GB). 

Figure <ref> depicts the strong scaling properties of the multi-batch L-BFGS method, for different batch sizes (r) and nodes (K=1,2,...,128). For this task, 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). One can observe that as the number of nodes (K) is increased, the compute times for the gradient and the L-BFGS direction decrease. 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, ∼18MB of data), the amount of local work is not sufficient to overcome the communication cost. 


   

   


  §.§.§ Weak Scaling –  Fixed Problem Dimension, Increasing Data Size


r0.4


  

    
    < g r a p h i c s >

  

  

Weak scaling of robust multi-batch L-BFGS on a problem with artificial data;  
  n=10^7 and d=10^4. Each sample has 10· K non-zero elements (size of local problem 1.5GB). 

In order to illustrate the weak scaling properties of the algorithm, we generate a data-matrix X ∈ℝ^n	× d (n=10^7, d = 10^4), and compute the gradient and the L-BFGS direction on a shared cluster with different number of MPI processes (K=1,2,...,128). Each sample has 10· K non-zero elements, thus for any K the size of local problem is roughly 1.5GB (for K=128 size of data 192GB). Effectively, the dataset size (n) is held fixed, but the sparsity of the data decreases as more MPI processes are used. 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 <ref>. On the other hand, because we are using a Vector-Free L-BFGS implementation <cit.> for computing the L-BFGS direction, the amount of time needed for each node to compute the L-BFGS direction decreases as K is increased. However, increasing K does lead to larger communication overhead, and as such the overall time needed to compute and communicate the L-BFGS direction increases slightly as K is increased. 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. 

 




  §.§.§ Increasing Problem Dimension, Fixed Data Size and K


r0.4


  

    
    < g r a p h i c s >

  

    

Scaling of robust multi-batch L-BFGS on a problem with artificial data; n=10^7, increasing d and K=8 MPI processes. Each sample had 200 non-zero elements (dataset size 29GB). 

In this experiment, we investigate the effect of a change in the dimension (d) of the problem on the computation of the gradient and the L-BFGS direction. We fix the size of data (29GB) and the number of MPI processes (K=8), and generate data with n=10^7 samples, where each sample has 200 non-zero elements. Figure <ref> 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. We should note that the results presented in Figure <ref> are not surprising; there is
minimal change in performance (in terms of the gradient computation) as dimension increases, since the number of nonzero elements is fixed and sparse matrix operations are emplyed.



























  

  
    


§ FINAL REMARKS




In this paper, we assumed that sample consistency is not possible (fault-tolerant setting) or desirable (multi-batch setting), and described a novel and robust variant of the L-BFGS method designed to deal with two adversarial situations. The success of the algorithm relies on the fact that gradient differences need not be computed on the full batch, rather a small subset can be used alleviating the need for double function evaluations while still maintaining useful curvature information. The method enforces a small degree of control in the sampling process and avoids the pitfalls of using inconsistent gradient differences by performing quasi-Newton updating on the overlap between consecutive samples. 

Our numerical results indicate that provided the overlap is not too small, the proposed method is efficient in practice on machine learning tasks such as binary classification logistic regression and neural network training. The experiments presented in this paper show that the empirical performance of the method matches that predicted by the theory for both strongly convex and nonconvex functions. Specifically, in the strongly convex case the multi-batch L-BFGS method with a constant step length converges to a neighborhood of the solution at a linear rate, and in the nonconvex case the iterates produced by the multi-batch L-BFGS method converge to a neighborhood of a stationary point.

Of course, the development, both theoretical and practical, of stochastic quasi-Newton methods is far from complete, and there are many interesting directions that can and should be investigated. Theoretical analysis that would suggest the batch size and overlap size would be of great interest in practice. Moreover, an investigation of the multi-batch L-BFGS method that employs variance reduced gradients in lieu of the stochastic gradients could have both theoretical and practical advantages. Finally, a stochastic line search that could work in conjunction with the multi-batch L-BFGS method would be novel both algorithmically and theoretically, and would most probably make the method even more competitive in practice.

  


  §.§.§ Acknowledgements
 This work was partially supported by DARPA Lagrange award HR-001117S0039 and U.S. National Science Foundation, under award numbers NSF:CCF:1618717, NSF:CMMI:1663256 and NSF:CCF:1740796.




gOMS






§ EXTENDED NUMERICAL RESULTS - REAL DATASETS - LOGISTIC REGRESSION



In this section, we present further numerical results on binary classification logistic regression problems, on the datasets listed in Table <ref>, 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.




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

    min_w ∈ℝ^d F(w) =  1/n∑_i=1^nlog(1+e^-y^i(w^Tx^i))
      + σ/2w^2,

where (x^i, y^i)_i=1^n  denote the training examples and σ = 1/n is the regularization parameter.



 §.§ Extended Numerical Results - Multi-Batch Setting


For the experiments in this section (Figures <ref>-<ref>), we ran four methods:

	
  * (Robust L-BFGS) robust multi-batch L-BFGS (Algorithm <ref>),
	
  * (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 <ref>,
	
  * (SGD) serial SGD; at every iteration one sample is used to compute the gradient.

In Figures <ref>-<ref> we 
show the evolution of ∇ F(w) for different step lengths α, and for various batch (| S| = r· n) and overlap (| O| = o ·| S|) sizes. Each Figure (<ref>-<ref>) consists of 10 plots that illustrate the performance of the methods with the following parameters:

	
  * Top 3 plots: α=1, o=20% and r=1%,5%,10%
	
  * Middle 3 plots: α=0.1, o=20% and r=1%,5%,10%
	
  * Bottom 4 plots: α=1, r=1% and o=5%,10%,20%,30%




































 §.§ Extended Numerical Results - Fault-Tolerant Setting



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/or 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 – On 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 ∈  R^n × n, then  synchronize all MPI processes and compute C=A· 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 <ref>. 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. In some cases, the difference between the minimum and maximum computation and communication times vary 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  <ref>-<ref> we present 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 α=0.1 and p∈{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.











 























§ EXTENDED NUMERICAL RESULTS - NEURAL NETWORKS



In this section, we present complete numerical results for the Neural Network training tasks on the MNIST[Available at: <http://yann.lecun.com/exdb/mnist/>.] and CIFAR10/CIFAR100[Available at: <https://www.cs.toronto.edu/ kriz/cifar.html>.] datasets. More specifically, we show the training and testing accuracy of the methods for all batch sizes (| S | ∈{50, 100, 200, 500, 1000, 2000, 4000 }). The tasks are summarized in Table <ref>, and the details are explained below.




	
  * MNIST Convex: This problem was a convex problem. It is a simple neural network with no hidden layers and soft-max cross-entropy loss function.
	
  * MNIST DNN (SoftPlus):
	This network consists of two convolutional layers and two fully connected layer. The first convolutional layer has 32 filters for each 5 × 5 patch, and the second convolutional layer has 64 filters for each 5 × 5 patch. Every convolutional layer is followed by a 2 × 2 max pooling layer. The first fully connected layer has 1024 neurons and the second fully connected layer produces a 10 dimensional output. We used SoftPlus activation functions.	
	The loss function is soft-max cross-entropy loss.
	
  * MNIST DNN (ReLU):
	Same as above, but with ReLU activation functions.
	
  * CIFAR10 LeNet:
	This network consists of two convolutional layers and three fully connected layer. The first convolutional layer has 6 filters for each 5 × 5 patch, and the second convolutional layer has 16 filters for each 5 × 5 patch. Every convolutional layer is followed by a 2 × 2 max pooling layer. The first fully connected layer has 120 neurons and the second fully connected layer
has 84 neurons. The last fully connected layer produces a 10 dimensional output. The activation functions used is ReLU  and the loss function used is soft-max cross-entropy loss. 
	
  * CIFAR10 VGG11 and CIFAR100 VGG11:
This is standard VGG11 network which contains 8 convolution layers, each followed by batch-normalization and ReLU activation functions. There are also 5 max-pooling layers and one average pooling layer. The output of desired dimension (10 or 100) is achieve by fully connected layer.  	
	The loss function is soft-max cross-entropy loss. 

 
For the experiments in this section (Figures <ref> and <ref>), we ran the following methods:

	
  * (LFBGS) multi-batch L-BFGS with the cautious strategy: Algorithm <ref> with standard initial scaling (γ_k I, where γ_k = s_k-1^T y_k-1/y_k-1^T y_k-1) or α I initial scaling,
	
  * (Adam) <cit.>,
	
  * (SGD) <cit.>.


Figures <ref> and <ref>  
show the evolution of training and testing accuracy for different batch sizes for the different neural network training tasks. Figures <ref> and <ref> show the two diagnostic metrics for different batch sizes.
































