Title: A FAST QUASI-NEWTON-TYPE METHOD FOR LARGE-SCALE STOCHASTIC OPTIMISATION

Abstract: During recent years there has been an increased interest in stochastic adaptations of limited memory quasi-Newton methods, which compared to pure gradientbased routines can improve the convergence by incorporating second order information. In this work we propose a direct least-squares approach conceptually similar to the limited memory quasi-Newton methods, but that computes the search direction in a slightly different way. This is achieved in a fast and numerically robust manner by maintaining a Cholesky factor of low dimension. This is combined with a stochastic line search relying upon fulfilment of the Wolfe condition in a backtracking manner, where the step length is adaptively modified with respect to the optimisation progress. We support our new algorithm by providing several theoretical results guaranteeing its performance. The performance is demonstrated on real-world benchmark problems which shows improved results in comparison with already established methods.

1. INTRODUCTION:

In learning algorithms we often face the classical and hard problem of stochastic optimisation where we need to minimise some non-convex cost function f (x) min

EQUATION

when we only have access to noisy evaluations of the function and its gradients. We take a particular interest in situations where the number of data and/or the number of unknowns d is very large.

The importance of this problem has been increasing for quite some time now. The reason is simple: many important applied problems ask for its solution, including most of the supervised machine learning algorithms when applied to large-scale settings. There are two important situations where the non-convex stochastic optimisation problem arise. Firstly, for large-scale problems it is often prohibitive to evaluate the cost function and its gradient on the entire dataset. Instead, it is divided into several mini-batches via subsampling, making the problem stochastic. This situation arise in most applications of deep learning. Secondly, when randomised algorithms are used to approximately compute the cost function and its gradients the result is always stochastic.

Our contributions are: 1. A new stochastic line search algorithm allowing for adaptive step-lengths akin to what is done in the corresponding state-of-the-art deterministic optimisation algorithms. This is enabled via a stochastic formulation of the first Wolfe condition. 2. We provide a new and efficient way of incorporating second-order (curvature) information into the stochastic optimiser via a direct least-squares approach conceptually similar to the popular limited memory quasi-Newton methods. 3. To facilitate a fast and numerically robust implementation we have derived tailored updating of a small dimension Cholesky factor given the new measurement pair (with dimension equal to the memory length). 4. We support our new developments by establishing several theoretical properties of the resulting algorithm. The performance is also demonstrated on real-world benchmark problems which shows improved convergence properties over current state-of-the-art methods.

2. RELATED WORK:

Due to its importance, the stochastic optimisation problem is rather well studied by now. The first stochastic optimisation algorithm was introduced by Robbins & Monro (1951) . It makes use of first-order information only, motivating the name stochastic gradient (SG), which is the contemporary term (Bottou et al., 2018) for these algorithms, originally referred to as stochastic approximation.

Interestingly most SG algorithms are not descent methods since the stochastic nature of the update can easily produce a new iterate corresponding to an increase in the cost function. Instead, they are Markov chain methods in that their update rule defines a Markov chain.

The basic first-order SG algorithms have recently been significantly improved by the introduction of various noise reduction techniques, see e.g. (Johnson & Zhang, 2013; Schmidt et al., 2013; Konečný & Richtárik, 2017; Defazio et al., 2014) .

The well-known drawback of all first-order methods is the lack of curvature information. Analogously to the deterministic setting, there is a lot to be gained in extracting and using second-order information that is maintained in the form of the Hessian matrix. The standard quasi-Newton method is the BFGS method, named after its inventors (Broyden, 1967; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970) . In its basic form, this algorithm does not scale to the large-scale settings we are interested in. The idea of only making use of the most recent iterates and gradients in forming the inverse Hessian approximation was suggested by Nocedal (1980) and Liu & Nocedal (1989) . The resulting L-BFGS method is computationally cheaper with a significantly reduced memory footprint. Due to its simplicity and good performance, this has become one of the most commonly used second-order methods for large-scale problems. Our developments makes use of the same trick underlying L-BFGS, but it is carefully tailored to the stochastic setting.

Over the past decade we have witnessed increasing capabilities of these so-called stochastic quasi-Newton methods, the category to which our developments belong. The work by Schraudolph et al. (2007) developed modifications of BFGS and its limited memory version. There has also been a series of papers approximating the inverse Hessian with a diagonal matrix, see e.g. Bordes et al. (2009) and Duchi et al. (2011) . The idea of exploiting regularisation together with BFGS was successfully introduced by Mokhtari & Ribeiro (2014) . Later Mokhtari & Ribeiro (2015) also developed a stochastic L-BFGS algorithm without regularisation. The idea of replacing the stochastic gradient difference in the BFGS update with a subsampled Hessian-vector product was recently introduced by Byrd et al. (2016) , and Wang et al. (2017) derived a damped L-BFGS method.

Over the past five years we have also seen quite a lot of fruitful activity in combining the stochastic quasi-Newton algorithms with various first-order noise reduction methods (Moritz et al., 2016; Gower et al., 2016) . A thorough and forward-looking overview of SG and its use within a modern machine learning context is provided by Bottou et al. (2018) . It also includes interesting accounts of possible improvements along the lines of first-order noise reduction and second-order methods.

It is interesting-and perhaps somewhat surprising-to note that it is only very recently that stochastic line search algorithms have started to become available. One nice example is the approach proposed by Mahsereci & Hennig (2017) which uses the framework of Gaussian processes and Bayesian optimisation. The step length is chosen that best satisfies a probabilistic measure combining reduction in the cost function with satisfaction of the Wolfe conditions. Conceptually more similar to our procedure is the line search proposed by Bollapragada et al. (2018) , which is tailored for problems that are using sampled mini-batches, as is common practice within deep learning.

3. STOCHASTIC LINE SEARCH:

In deterministic line search algorithms we first compute a search direction p k and then decide how far to move along that direction according to

EQUATION

where α k > 0 is referred to as the step length. The search direction is of the form

EQUATION

where H k denotes an approximation of the inverse Hessian matrix. The question of how far to move in this direction can be framed as the following scalar minimisation problem

EQUATION

The most common way of dealing with this problem is to settle for a possibly sup-optimal solution that guarantees at least a sufficient decrease. Such at solution can be obtained by selecting a step length α k that satisfies the following inequality

EQUATION

where c ∈ (0, 1). The above condition is known as the first Wolfe condition (Wolfe, 1969; 1971) or the Armijo condition (Armijo, 1966) .

While deterministic line search algorithms have been well established for a long time, their stochastic counterparts are still to a large extent missing. Consider the case when the measurements of the function and its gradient are given by

EQUATION

where g k ∇f (x)| x=x k , and e k ∈ R and v k ∈ R d×1 denote noise on the function and gradient evaluations, respectively. Furthermore we assume that

EQUATION

)

The challenge is that since f k and g k are random variables it is not obvious how to select a step length α k that-in some sense-guarantees a descent direction.

We explore the idea of requiring equation 5 to be fulfilled in expectation when the exact quantities are replaced with their stochastic counterparts,

EQUATION

where p k = -H k g k . This is certainly one way in which we can reason about the Wolfe condition in the stochastic setting we are interested in. Although satisfaction of equation 8 does not leave any guarantees when considering a single measurement, it still serves as an important property that could be exploited to provide robustness for the entire optimisation procedure. To motivate our proposed algorithm we hence start by establishing the following results.

lemma 1 (Stochastic Wolfe condition 1) Assume that i) the gradient estimates are unbiased , ii) the cost function estimates are possibly biased , and iii) a descent direction is ensured in expectation E p T k g k < 0. Then (for small α:s)

EQUATION

Proof 1 See Appendix A.1.

lemma 2 There exists a step length

α k > 0 such that E f (x k + α k p k ) -f (x k ) < 0. Proof 2 See Appendix A.2

Relying upon these results we propose a line search based on the idea of repeatedly decreasing the step length α k until the stochastic version of the first Wolfe condition is satisfied. Pseudo-code for this procedure is given in Algorithm 1. An input to this algorithm is the search direction p k , which can be computed using any preferred method. The step length is initially set to be the minimum of the "natural" step length 1 and the iteration dependent value ξ/k. In this way the initial step length is kept at 1 until k > ξ, a point after which it is decreased at the rate 1/k. Then we check whether the new point x k + α p k satisfies the stochastic version of the first Wolfe condition. If this is not the case, we decrease α k with the scale factor ρ. This is repeated until the condition is met, unless we hit an upper bound max{0, τ -k} on the number of backtracking iterations, where τ > k is a positive integer. With this restriction the decrease of the step length is limited, and when k ≥ τ we use the initial step length no matter if the Wolfe condition is satisfied or not. The purpose of ξ is to facilitate convergence by reducing the step length as the minima is approached. The motivation of τ comes from arguing that the backtracking loop does not contribute as much when the step length is small, and hence it is reasonable to provide a limit on its reduction.

Algorithm 1 Stochastic backtracking line search:

Require: Iteration index k, spatial point x k , search direction p k , scale factor ρ ∈ (0, 1), reduction limit ξ ≥ 1, backtracking limit τ > 0.

1: Set the initial step length α k = min{1, ξ/k} 2: Set i = 1 3: while f (x k + α p k ) > f (x k ) + cα k g T k p k and i ≤ max{0, τ -k} do 4:

Reduce the step length α k ← ρα k 5:

Set i ← i + 1 6: end while 7: Update

x k+1 = x k + α k p k

4. COMPUTING THE SEARCH DIRECTION:

In this section we address the problem of computing a search direction based on having a limited memory available for storing previous gradients and associated iterates. The approach we adopt is similar to limited memory quasi-Newton methods, but here we employ a direct least-squares estimate of the inverse Hessian matrix rather than more well-known methods such as damped L-BFGS and L-SR1. In contrast to traditional quasi-Newton methods, this approach does not strictly impose neither symmetry or fulfilment of the secant condition. It may appear peculiar to relax these requirements. However, in this setting is not obvious that enforced symmetry necessarily produces a better search direction. Furthermore, the secant condition relies upon the rather strong approximation that the Hessian matrix is constant between two subsequent iterations. Treating the condition less strictly might be helpful when that approximation is poor, perhaps especially in a stochastic environment. We construct a limited-memory inverse Hessian approximation in Section 4.1 and show how to update this representation in Section 4.2. In Section 4.3 we discuss a particular design choice associated with the update and Section 4.4 provides a means to ensure that a descent direction is calculated. The complete algorithm combining the procedure presented in this section with the line search routine in Algorithm 1 is given in Appendix C.

4.1. QUASI-NEWTON INVERSE HESSIAN APPROXIMATIONS:

According to the secant condition (see e.g. Fletcher (1987) ), the inverse Hessian matrix H k should satisfy

EQUATION

where

y k = g k -g k-1 and s k = x k -x k-1 .

Since there are generally more unknown values in H k than can be determined from y k and s k alone, quasi-Newton methods update H k from a previous estimate H k-1 by solving regularised problems of the type

EQUATION

where

X 2 F,W = XW 2 F = trace(W T X T XW

) and the choice of weighting matrix W results in different algorithms (see Hennig (2015) for an interesting perspective on this). Examples of the most common quasi-Newton methods are given in Appendix B.

We employ a similar approach and determine H k as the solution to the following regularised leastsquares problem

EQUATION

where Y k and S k hold a limited number of past y k 's and s k 's according to

EQUATION

Here, m << d is the memory limit and y k = g k -g k-1 is an estimate of y k . The regulator matrix Hk acts as a prior on H and can be modified at each iteration k. The parameter λ > 0 is used to control the relative cost of the two terms in equation 12. It can be verified that the solution to the above least-squares problem ( 12) is given by

EQUATION

)

where I denotes the identity matrix. The above inverse Hessian estimate can be used to generate a search direction in the standard manner by scaling the negative gradient, that is

EQUATION

However, for large-scale problems this is not practical since it involves the inverse of a large matrix.

To ameliorate this difficulty, we adopt the standard approach by storing only a minimal (limited memory) representation of the inverse Hessian estimate H k . To describe this, note that the dimensions of the matrices involved are

EQUATION

)

We can employ the Sherman-Morrison-Woodbury formula to arrive at the following equivalent expression for H k

EQUATION

Importantly, the matrix inverse λI + Y T k Y k -1 is now by construction a positive definite matrix of size m × m. Therefore, we construct and maintain a Cholesky factor of λI + Y T k Y k since this leads to efficient solutions. In particular, if we express this matrix via a Cholesky decomposition

EQUATION

where R k ∈ R m×m is an upper triangular matrix, then the search direction p k = -H k g k can be computed via

EQUATION

Note that the above expressions for p k involve first computing w k , which itself involves a computationally efficient forward-backward substitution (recalling that R k is an m × m upper triangular matrix and the memory length m is typically 10-50). Furthermore, Hk is typically diagonal so that Hk z k is also efficient to compute. The remaining operations involve four matrix-vector products and two vector additions. Therefore, for problems where d >> m then the matrix-vector products will dominate the computational cost.

Constructing R k can be achieved in several ways. The so-called normal-equation method constructs the (upper triangular) part of λI + Y T k Y k and then employs a Cholesky routine, which produces R k in O(n m(m+1) 2 +m 3 /3) operations. Alternatively, we can compute R k by applying Givens rotations or Householder reflections to the matrix

M k = √ λI Y T k T . This costs O(2m 2 ((n + m) -m/3))

operations, and is therefore more expensive, but typically offers better numerical accuracy (Golub & Van Loan, 2012) .

4.2. FAST AND ROBUST INCLUSION OF NEW MEASUREMENTS:

In order to maximise the speed, we have developed a method for updating a Cholesky factor given the new measurement pair (s k+1 , y k+1 ). Suppose we start with a Cholesky factor R k at iteration k such that

EQUATION

Assume, without loss of generality, that Y k and S k are ordered in the following manner

EQUATION

where Y 1 , Y 2 , S 1 and S 2 are defined as

EQUATION

EQUATION

and is an appropriate integer so that Y k and S k have m columns. The above ordering arises from "wrapping-around" the index when storing the measurements. We create the new Y k+1 and S k+1 by replacing the oldest column entries, y k-m+1 and s k-m+1 , with the latest measurements y k+1 and s k+1 , respectively, so that

EQUATION

)

The aim is to generate a new Cholesky factor R k+1 such that

EQUATION

To this end, let the upper triangular matrix R k be written conformally with the columns of Y k as

EQUATION

so that R 1 and R 2 have the same number of columns as Y 1 and Y 2 , respectively. Furthermore, r 1 is a column vector, r 2 is a scalar and r 3 is a row vector. Therefore,

EQUATION

)

By observing a common structure for the update λI

EQUATION

where r 4 , r 5 and r 6 are determined by

EQUATION

)

The final term R 6 can be obtained by noticing that

EQUATION

which implies

EQUATION

Therefore R 6 can be obtained in a computationally very efficient manner by down-dating and updating the Cholesky factor R 4 with the rank-1 matrices r T 6 r 6 and r T 3 r 3 , respectively (see e.g. Section 12.5.3 in Golub & Van Loan (2012) ).

4.3. SELECTING Hk:

There is no magic way of selecting the prior matrix Hk . However, in practise it has proved very useful to employ a simple strategy of Hk γ k I, where the positive scalar γ k > 0 is adaptively chosen in each iteration. As a crude measure of progress we adopt the following rule

EQUATION

where κ ≥ 1 is a scale parameter, and q corresponds to the number of backtracking loops in the line search; the values κ = 1.3 and q = 3 were found to work well in practise. The intuition behind equation 31 is that if no modification of the step length α k is made, we can allow for a more "aggressive" regularisation. Note that a low γ k is favouring small elements in H k . Since p k = -H g k , this limits the magnitude of p k and the change x k+1 -x k is kept down. Hence, it is good practice to set the initial scaling γ 0 relatively small and then let it scale up as the optimisation progresses. Furthermore, we should point out that a diagonal Hk comes with an efficiency benefit, since the product Hk z k in equation 19 then is obtained as the element-wise product between two vectors.

4.4. ENSURING A DESCENT DIRECTION:

In deterministic quasi-Newton methods, the search direction p k must be chosen to ensure a descent direction such that p T k g k < 0, since this guarantees reduction in the cost function for sufficiently small step lengths α k . Since p k = -Hg k , we have that p T k g k = -g T k H k g k which is always negative if the approximation H k of the inverse Hessian is positive definite. Otherwise, we can modify the search direction by subtracting a multiple of the gradient p k ← p k -β k g k . This is motivated by noticing that

EQUATION

which always can be made negative by selecting β k large enough, i.e. if

EQUATION

)

In the stochastic setting, the condition above does not strictly enforce a descent direction. Hence the search direction p k as determined by equation 15 is not a descent direction in general. However, ensuring that the condition is fulfilled in expectation is necessary since this is one of the assumptions made in lemma 1. Hence, we now establish the following result.

EQUATION

then

EQUATION

Proof 3 See Appendix A.3.

In the practical setting we can not use equation 34 as it is, since we do not have access to g k and p k , nor the noise variance σ 2 g . Instead we suggest a pragmatic approach in which these quantities are replaced by their estimates g k , p k and σ 2 g . The noise estimate σ 2 g could either be regarded a design parameter or empirically calculated from one or more sets of repeatedly collected gradient measurements. Nevertheless, picking β k sufficiently large ensures fulfilment of equation 32 and equation 35 simultaneously.

5. NUMERICAL EXPERIMENTS:

Let us now put our new developments to the test on a suite of problems from different categories to exhibit different properties and challenges. In Section 5.1 we study a commonly used benchmark, namely the collection of logistic classification problems described by Chang & Lin (2011) in the form of their library for support vector machines (LIBSVM). In Section 5.2 we consider an optimisation problem arising from the use of deep learning to solve the classical machine learning benchmark MNIST1 , where the task is to classify images of handwritten digits. Also, we test our method on training a neural network on the CIFAR-10 dataset (Krizhevsky, 2009) .

In our experiments we compare against relevant state-of-the-art methods. All experiments were run on a MacBook Pro 2.8GHz laptop with 16GB of RAM using Matlab 2018b. All routines where programmed in C and compiled via Matlab's mex command and linked against Matlab's Level-1,2 BLAS libraries. More details about the experiments are available in Appendix D. The source code used to produce the results will be made freely available.

5.1. LOGISTIC LOSS AND A 2-NORM REGULARISER:

The task here is to solve eight different empirical risk minimisation problems using a logistic loss function with an L2 regulariser (two are shown here and all eight are profiled in Appendix D). The data is taken from Chang & Lin (2011) . These problems are commonly used for profiling optimisation algorithms of the kind introduced in this paper, facilitating comparison with existing state-of-the-art algorithms. More specifically, we have used a similar set-up as Gower et al. (2016) , which inspired this study. The chosen algorithm parameters for each case is detailed in Appendix D.

We compared our limited memory least-squares approach (denoted as LMLS) against two existing methods from the literature, namely, the limited memory stochastic BFGS method after Bollapragada et al. ( 2018 

5.2. DEEP LEARNING:

Deep convolutional neural networks (CNNs) with multiple layers of convolution, pooling and nonlinear activation functions are delivering state-of-the-art results on many tasks in computer vision.

We are here borrowing the stochastic optimisation problems arising in using such a deep CNN to solve the MNIST and CIFAR-10 benchmarks. The particular CNN structure used for the MNIST example employs 5 × 5 convolution kernels, pooling layers and a fully connected layer at the end. We made use of the publicly available code provided by Zhang (2016) , which contains all the implementation details. For the CIFAR-10 example, the network includes 13-layers with more than 150,000 weights, see Appendix D for details. The MATLAB toolbox MatConvNet (Vedaldi & Lenc, 2015) was used in the implementation. Figures 1c and 1d show the average cost versus time for 10 Monte-Carlo trials with four different algorithms: 1. the method developed here (LMLS), 2. a stochastic limited memory BFGS method after Bollapragada et al. ( 2018) (denoted LBFGS), 3. Adam developed by Kingma & Ba (2015) , and 4. stochastic gradient (denoted SG). Note that all algorithms make use of the same gradient code.

6. CONCLUSION AND FUTURE WORK:

In this work we have presented a least-squares based limited memory optimisation routine that benefits from second order information by approximating the inverse Hessian matrix. The procedure is conceptually similar to quasi-Newton methods, although we do not explicitly enforce symmetry or satisfaction of the secant condition. By regularising with respect to an inverse Hessian prior, we allow for an adaptive aggressiveness in the search direction update. We have shown that the computations can be made robust and efficient using tailored Cholesky decompositions, with a cost that scales linearly in the problem dimension. Our method is designed for stochastic problems through a line search that repeatedly decreases the step length so as to satisfy the first Wolfe condition in expectation. Theoretical results have been established that support the proposed procedure. The method shows improved convergence properties over existing algorithms when applied to benchmark problems of various size and complexity.

A PROOFS OF THE LEMMAS A.1 THE FIRST STOCHASTIC WOLFE CONDITION:

In the interest of a simple notation we note that we can drop the sub-index k from all variables since all reasoning in this proof is for iteration k. We also drop x from the arguments and write f rather than f (x), etc. However, we do make the randomness explicit, since this is crucial for arriving at the correct answer. example we write f z to indicate that the random variable z was used in computing the estimate of the function value at the iterate x k . This means that the noise contaminating the measurement corresponds to the realisation of z.

Consider the stochastic version of the Wolfe condition

EQUATION

where the expected value is w.r.t. the randomness used in computing the required estimators as indicated by z and z . Using Taylor series we can for small step lengths α express f z (x + α p z ) according to

EQUATION

)

Here it is crucial to note that the randomness in the above estimate stem from two different sources. The search direction p z enters in the argument of the Taylor expansion which is why it has to be computed before actually performing the Taylor expansion, implying that the randomness is due to z and not z for the search direction. Based on equation 36 we have the following expression for the stochastic Wolfe condition

EQUATION

Let b denote a possible bias in the cost function estimator, then we have that E z f z = f + b and

E z f z = f + b.

Hence, the first two terms in equation 38 cancel and we have

c ≤ E z,z p T z g z E z,z [ p T z g z ] = E z p T z E z [ g z ] E z [ p T z g z ] = p T g E z [ p T z g z ]

.

(39)

The denominator can be written as

EQUATION

Using the definition of covariance we have that

EQUATION

where we in the last equality made use of the assumption that the gradients are unbiased and that Cov g = σ g I. Summarising we have that

EQUATION

and hence we get

EQUATION

)

Recalling the assurance given in lemma 3 with the associated proof in Appendix A.3, we can always modify p to guarantee that this bound is positive.

A.2 REDUCTION IN COST FUNCTION:

In this section we show that there exists an α > 0 such that

EQUATION

We do this by studying the difference

EQUATION

in which the first term can be expressed using a Taylor series according to

EQUATION

)

For small of α we can discard the term O(α p z 2 ) from this expression. Hence we have that

EQUATION

Following from Section 4 we can ensure p T g to be negative, and hence equation 44 holds as long as α is chosen small enough for the Taylor expansion to be a valid approximation.

A.3 ENSURING A DESCENT DIRECTION:

Note that

EQUATION

Using equation 42 we have

EQUATION

and it directly follows from this that

EQUATION

Hence

EQUATION

and we see that

EQUATION

B QUASI-NEWTON METHODS

The historically most popular quasi-Newton method is given by the BFGS algorithm (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970 )

EQUATION

A closely related version is obtained if the optimisation in equation 11 is done with respect to the Hessian matrix rather than to its inverse. This results in the so-called DFP formula

EQUATION

)

Both of these are rank 2 updates that ensure H k to be positive definite, an attractive feature in the sense that it guarantees a descent direction. However, it may cause problems in regions where the true inverse Hessian is indefinite. Another well-known alternative is the symmetric rank 1 (SR1) method

EQUATION

)

Apart from BFGS and DFP above, this is a rank 1 update that in general does not preserve positive definiteness. For this reason it has received a particular interest within the so-called trust-region framework, where its ability of producing indefinite inverse Hessian approximations is being regarded as a major strength (Nocedal & Wright, 2006) .

A main drawback of quasi-Newton methods is that they scale poorly to large problems, since the number of elements in H k equals the square of the input dimension d. This has resulted in developments of so-called limited memory algorithms, which rely upon the idea of forming the search direction p k directly without the need of storing the full matrix H k . To that end a number m << d of past differences y and s are being stored in memory. A notable member of this family is the L-BFGS algorithm (Nocedal, 1980) , which has been widely within large-scale optimisation. During recent years, not at least due to the growing number of deep learning applications, there has been an increasing interest in adapting deterministic limited memory methods to the stochastic setting (Wang et al., 2017; Moritz et al., 2016; Gower et al., 2016; Schraudolph et al., 2007; Bordes et al., 2009; Mokhtari & Ribeiro, 2015; Byrd et al., 2016; Bollapragada et al., 2018) .

C RESULTING ALGORITHM:

We summarise our ideas from Section 3 and 4 in Algorithm 2. The if-statement on line 5 is included to provide a safe-guard against numerical instability.

D EXPERIMENT DETAILS:

Details for the datasets used in Section 5 are listed in Table 1 , including the parameter choices we made in our algorithm. Here, b denotes the mini-batch size. The results of all logistic regression experiments are collected in Figure 2 (including those already shown in Figure 1 ), and the neural network results in Figures 1d and 1c are provided in a more readable size in Figure 3a and 3b. Detailed information of the network structure in the CIFAR experiment is provided in the printout shown in Figure 4 .

The adaptive procedure of selecting the inverse Hessian prior Hk did not have much impact in most of the logistic regression examples, and thus it was kept constant Hk = γ 0 I except for in the covtype and URL problems. In the CIFAR and MNIST problems, however, the procedure was found to be of significant importance. Require: An initial estimate x 1 , a maximum number of iterations k max , memory limit m, regularisation parameter λ, scale factor ρ ∈ (0, 1), reduction limit ξ ≥ 1, backtracking limit τ > 0 1: Set k = 1 2: while k < k max do 3:

Obtain a measurement of the cost function and its gradient

f k = f (x k ) + e k , g k = g k + v k .

4:

Set Hk = γ k I where

γ k = κγ k-1 , if α k-1 = 1, γ k-1 /κ, if α k-1 < 1/ρ q , γ k-1 ,

otherwise. Select p k as

p k = -Hk z k -λ -1 S k (Y T k z k ), z k = g k -Y k w k , w k = R -1 k R -T k Y T

k g k , with details provided in Section 4.2. Set k ← k + 1 24: end while LMLS refers to limited memory least-squares approach developed in this paper, SG refers to basic stochastic gradient and Adam refers to Kingma & Ba (2015) .

:

yann.lecun.com/exdb/mnist/

Figures and Tables:

Figure 1: Performance on two classification tasks using a logistic loss with a two-norm regulariser ((a) and (b)), and two deep convolutional neural networks (CNNs) used for recognising images of handwritten digits from the MNIST data (c), and classification of the images in the CIFAR-10 data (d). Lighter shaded lines indicate individual runs, whereas the darker shaded line indicates the average.

) (denoted as LBFGS) and the stochastic variance reduced gradient descent (SVRG) by Johnson & Zhang (2013) (denoted SVRG). Figures 1a and 1b show the cost versus time for 50 Monte-Carlo experiments.

Y k and S k by replacing the oldest vector-pairs in Y k-1 and S k-1 with y k and s k Y k and S k by adding y k and s k to Y k-1 and S kY k = Y k-1 and S k = S k

15:Setp k ← p k -β k g k with β k where while f (x k + α p k ) > f (x k ) + cα k g T k p k and i ≤ max{0, τ -k} do 19:Reduce the step length α k ← ρα k x k+1 = x k + α k p k 23:

Figure 2: Performance on classification tasks using a logistic loss with a two-norm regulariser.

Figure 3: Solving the optimisation problem used in training a state-of-the-art deep convolutional neural network (CNN) used for recognising images in the (a) MNIST and (b) CIFAR-10 dataset.LMLS refers to limited memory least-squares approach developed in this paper, SG refers to basic stochastic gradient and Adam refers toKingma & Ba (2015).

List of the datasets used in the experiments, where n denotes the size of the dataset (column 2) and d denotes the number of variables in the optimisation problem (column 3). The remaining columns list our design parameters, including the mini-batch size b.

Table content: <table><tr><td>Problem</td><td>n</td><td>d</td><td>b</td><td>m</td><td>λ</td><td>ξ</td><td>τ</td><td>γ 0</td></tr><tr><td>gisette</td><td>6 000</td><td>5 000</td><td>770</td><td>20</td><td>10 -4</td><td>50</td><td>10</td><td>10</td></tr><tr><td>covtype</td><td>581 012</td><td>54</td><td>7 620</td><td>55</td><td>10 -4</td><td>100</td><td>5</td><td>100</td></tr><tr><td>HIGGS</td><td>11 000 000</td><td>28</td><td>66 340</td><td>28</td><td>10 -4</td><td>20</td><td>5</td><td>100</td></tr><tr><td>SUSY</td><td>3 548 466</td><td>18</td><td>5 000</td><td>18</td><td>10 -4</td><td>50</td><td>10</td><td>100</td></tr><tr><td>epsilon</td><td>400 000</td><td>2 000</td><td>1 000</td><td>20</td><td>10 -4</td><td>150</td><td>1</td><td>100</td></tr><tr><td>rcv1</td><td>20 242</td><td>47 236</td><td>710</td><td>10</td><td>10 -4</td><td>50</td><td>10</td><td>600</td></tr><tr><td>URL</td><td>2 396 130</td><td>3 231 961</td><td>1548</td><td>5</td><td>10 -4</td><td>200</td><td>50</td><td>100</td></tr><tr><td>spam</td><td>82 970</td><td>823 470</td><td>2048</td><td>2</td><td>10 -4</td><td>150</td><td>10</td><td>200</td></tr><tr><td>MNIST</td><td>60 000</td><td>3898</td><td>1000</td><td>50</td><td>8 • 10 -4</td><td>150</td><td>50</td><td>1</td></tr><tr><td>CIFAR</td><td>50 000</td><td>150 000</td><td>200</td><td>20</td><td>10 -2</td><td>100</td><td>10</td><td>0.5</td></tr></table>

Peer review (iclr2019):
 This paper presents a new quasi-Newton method for stochastic optimization that solves a regularized least-squares problem to approximate curvature information that relaxes both the symmetry and secant conditions typically ensured in quasi-Newton methods. In addition to this, the authors propose a stochastic Armijo backtracking line search to determine the steplength that utilizes an initial steplength of 1 but switches to a diminishing steplength in later iterations. In order to make this approach computationally tractable, the authors propose updating and maintaining a Cholesky decomposition of a crucial matrix in the Hessian approximation. Although it is a good attempt at developing a new method, the paper ultimately lacks a convincing explanation (both theoretical and empirical) supporting their ideas, as I will critique below. 1. Stochastic Line Search Determining a steplength in the stochastic setting is a difficult problem, and I appreciate the authors’ attempt to attack this problem by looking at stochastic line searches. However, the paper lacks much detail and rigorous reasoning in the description and proof for the stochastic line search. First, the theory gives conditions that the Armijo condition holds in expectation. Proving anything about stochastic line searches is particularly difficult, so I’m on board with proving a result in expectation and doing something different in practice. However, much of the detail on how this is implemented in practice is lacking. How are the samples chosen for the line search? If we go along with the proposed theory, then when the function is reevaluated in the line search, a new sample is used. If this is the case, can one guarantee that the practical Armijo condition will hold? How often does the line search fail? How does the choice of the samples affect the cost of evaluating the line search? The theory also suggests that the particular choice of c is dependent on each iteration, particularly the inner product between the true search direction and the true gradient at iteration k. Does this allow for a fixed c to be used? How is c chosen? Is it fixed or adaptive? What happens as the true gradient approaches 0? The algorithm also places a limit on the number of backtracks permitted that decreases as the iteration count increases. What does the algorithm do when the line search fails? Does one simply take the step although the Armijo condition does not hold? In deterministic optimization, BFGS typically needs a smaller steplength in the beginning as the algorithm learns the scale of the problem, then eventually accepts the unit steplength to obtain fast local convergence. The line search proposed here uses an initial steplength of so that in early iterations, a steplength of 1 is used and in later iterations the algorithm uses a steplength. When this is combined with the diminishing maximum number of backtracking iterations, this will eventually yield an algorithm with a steplength of . Why is this preferred? Are the other algorithms in the numerical experiments tuned similarly? The theory also asks for a descent direction to be ensured in expectation. However, it is not the case that , so it is not correct to claim that a descent direction is ensured in expectation. Rather, the condition is requiring the angle between the negative stochastic gradient direction and search direction to be acute in expectation. All the proofs also depend on a linear Taylor approximation that is not well-explained, and I’m wary of proofs that utilize approximations in this way. Indeed, the precise statement is that , where . How does this affect the proof? Lastly, I would propose for the authors to change the name of their condition to the “Armijo condition” rather than using the term “1st Wolfe condition” since the Wolfe condition is typically associated with the curvature condition (p_k’ g_new >= c_2 p_k’ g_k), hence referring to a very different line search. 2. Design of the Quasi-Newton Matrix The authors develop an approach for designing the quasi-Newton matrix that does not strictly impose symmetry or the secant condition. The authors claim that this done because “it is not obvious that enforced symmetry necessarily produces a better search direction” and “treating the [secant] condition less strictly might be helpful when [the Hessian] approximation is poor”. This explanation seems insufficient to me to explain why relaxing these conditions via a regularized least-squares approach would yield a better algorithm, particularly in the noisy or stochastic setting. The lack of symmetry seems particularly strange; one would expect the true Hessian in the stochastic setting to still be symmetric, and one would still expect the secant condition to hold if the “true” gradients were accessible. It is also unclear how this approach takes advantage of the stochastic structure that exists within the problem. Additionally, the quasi-Newton matrix is defined based on the solution of a regularized least squares problem with a regularization parameter lambda. It seems to me that the key to the approximation is the balance between the two terms in the objective. How is lambda chosen? What is the effect of lambda as a tuned parameter, and how does it affect the quality of the Hessian approximation? It is unclear to me how this could be chosen in a more systematic way. The matrix also does not ensure positive definiteness, hence requiring a multiple of the gradient direction to be added to the search direction. In this case, the key parameter beta must be chosen carefully. What is a typical value of beta that is used for each of these problems? One would hope that beta is small, but if it is large, it may suggest that the search direction is primarily dominated by the stochastic gradient direction and hence the quasi-Newton matrix is not useful. The interplay of these different parameters needs to be investigated carefully. Lastly, since (L-)BFGS use a weighted Frobenius norm, I am curious why the authors decided to use a non-weighted Frobenius norm to define the matrix. How does changing the norm affect the Hessian approximation? All of these questions place the onus on the numerical experiments to see if these relaxations will ultimately yield a better algorithm. 3. Numerical Experiments As written, although the range of problems is broad and the numerical experiments show much promise, I do not believe that I could replicate the experiments conducted in the paper. In particular, how is SVRG and L-BFGS tuned? How is the steplength chosen? What (initial) batch sizes are used? Is the progressive batching mechanism used? (If the progressive batching mechanism is not used, then the authors should refer to the original multi-batch paper by Berahas, et al. [1] which do not increase the batch size and use a constant steplength.) In addition, a more fair comparison would include the stochastic quasi-Newton method in [2] that also utilize diminishing steplengths, which use Hessian-vector products in place of gradient differences. Multi-batch L-BFGS will only converge if the batch size is increased or the steplength diminished, and it’s not clear if either of these are done in the paper. Typos/Grammatical Errors: - Pg. 1: Commas are needed in some sentences, i.e. “Firstly, for large scale problems, it is…”; “…compute the cost function and its gradients, the result is…” - Pg. 2: “Interestingly, most SG algorithms…” - Pg 3: Remove “at least a” in second line - Pg. 3: suboptimal, not sup-optimal - Pg. 3: “Such a solution”, not “Such at solution” - Pg. 3: Capitalize Lemma - Pg. 4: fulfillment, not fulfilment - Pg. 7: Capitalize Lemma - Pg. 11: Before (42), Cov .hat{g} = .sigma_g^2 I - Pg. 11: Capitalize Lemma Summary: In summary, although the ideas appear to provide better numerical performance, it is difficult to evaluate if the ideas proposed in this paper actually yield a better algorithm. Many algorithmic details are left unanswered, and the paper lacks mathematical or empirical evidence to support their claims. More experimental and theoretical work is needed before the manuscript can be considered for publication. References: [1] Berahas, Albert S., Jorge Nocedal, and Martin Takác. "A multi-batch l-bfgs method for machine learning." Advances in Neural Information Processing Systems. 2016. [2] Byrd, Richard H., et al. "A stochastic quasi-Newton method for large-scale optimization." SIAM Journal on Optimization 26.2 (2016): 1008-1031. [3] Schraudolph, Nicol N., Jin Yu, and Simon Günter. "A stochastic quasi-Newton method for online convex optimization." Artificial Intelligence and Statistics. 2007.