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):
 The authors present an interesting variation of the standard QN methods. Their main point of departure from LBFGS/SR1 is in constructing a simpler Hessian inverse approximation. Recall that SR1 and LBFGS updates all satisfy the secant equation for each of the `m` previous gradient differences stored in memory. The authors choose to get "close" to satisfying the equations by solving an l_2 penalization of the secant equations. The resulting algorithm is interesting, but it is not clear from the paper what the claimed advantage of doing this is. The LBFGS and SR1 unrolled update rules for H (Hessian inverse approximation) is O(m^2 d) (Sec 7.2 NW 2006), and this seems to be the same for the authors* method, where the main matrix R_k that forms H has the same order. (BTW, did you mean *d* in place of *n* the computational order discussion preceding Sec 4.2?) The experiments show that this method*s performance is impressive compared to an LBFGS implementation provided by Bollapragada 2018 , but as I recall that paper presented a variable/increasing batch method, while the authors* method uses fixed batches (as far as I can tell) so it is not clear that comparison on time alone is sufficient. The advantage over LBFGS and SGD seen in MNIST seems to go away by the CIFAR example, so it is unclear what might happen in larger problems like ImageNet. I am also not able to see the difference between the *stochastic* line search presented here and the standard backtracking method as applied to mini-batch evaluated estimates. What is different, and new that brings in consideration for the noise? I recall that bollapragada 2018 had an additional variance based rule to check. Some more conservative values are chosen for the step length, but I do not see the justification presented in the appendix, esp Eq 47 : p is not independent from g here, being calucated as p=Hg, so E[p^t g] is not equal to the product of the individual expectations. Some key points were left out in the discussion of the experiments. This is a common slip up when writing conf papers these days, but please do consider discussing the settings of parameters like mini-batches sizes , value of .lambda in the H derivation, how one calculates the .sigma^2_g within the algorithm presented in the Appendix. The last must include an extra computational cost, or are you using Adam style online variance estimator? The MSE error alone seems insufficient in the results. Please publish the test mis-classification results too. Also, why is the MSE loss used with the softmax in CIFAR? Shouldn*t cross-entropy, better justified theoretically, be better justified?