Year: 2018

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.

heading: 1 INTRODUCTION
text: 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 x∈Rd f(x), (1)
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.

heading: 2 RELATED WORK
text: 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; Koneč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 quasiNewton 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.

heading: 3 STOCHASTIC LINE SEARCH
text: In deterministic line search algorithms we first compute a search direction pk and then decide how far to move along that direction according to
xk+1 = xk + αkpk, (2)
where αk > 0 is referred to as the step length. The search direction is of the form
pk = −Hkgk, (3)
where Hk 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
min α f(xk + αpk), α > 0. (4)
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
f(xk + αpk) ≤ f(xk) + cαkgTk pk, (5)
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
f̂k = f(xk) + ek, ĝk = gk + vk, (6)
where gk , ∇f(x)|x=xk , and ek ∈ R and vk ∈ Rd×1 denote noise on the function and gradient evaluations, respectively. Furthermore we assume that
E [ek] = b, Cov[ek] = σ 2 f , E [vk]= 0, Cov[vk] = σ 2 gI. (7a)
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,
E [ f̂(xk + αp̂k)− f̂(xk)− cαkĝTk p̂k ] ≤ 0, (8)
where p̂k = −Hkĝ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̂Tk ĝk ] < 0. Then (for small α:s)
E [ f̂(xk + αp̂k) ] ≤ E [ f̂(xk) + cαkĝ T k p̂k ] , where 0 < c < c̄ =
pTkgk pTkgk − σ2g Tr (H) . (9)
Proof 1 See Appendix A.1.
lemma 2 There exists a step length αk > 0 such that E [ f̂(xk + αkp̂k)− f̂(xk) ] < 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 xk + α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 xk, 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̂(xk + αp̂k) > f̂(xk) + cαkĝTk 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 xk+1 = xk + αkp̂k

heading: 4 COMPUTING THE SEARCH DIRECTION
text: 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 LBFGS 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.

heading: 4.1 QUASI-NEWTON INVERSE HESSIAN APPROXIMATIONS
text: According to the secant condition (see e.g. Fletcher (1987)), the inverse Hessian matrix Hk should satisfy
Hkyk = sk, (10)
where yk = gk − gk−1 and sk = xk − xk−1. Since there are generally more unknown values in Hk than can be determined from yk and sk alone, quasi-Newton methods update Hk from a previous estimate Hk−1 by solving regularised problems of the type
Hk = arg min H
‖H −Hk−1‖2F,W s.t. H = HT, Hyk = sk, (11)
where ‖X‖2F,W = ‖XW‖2F = trace(WTXTXW ) 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 Hk as the solution to the following regularised leastsquares problem
Hk = arg min H ‖HYk − Sk‖2F + λ‖H − H̄k‖2F , (12)
where Yk and Sk hold a limited number of past ŷk’s and sk’s according to
Yk , [ŷk−m+1, . . . , ŷk] , Sk , [sk−m+1, . . . , sk] . (13)
Here, m << d is the memory limit and ŷk = ĝk − ĝk−1 is an estimate of yk. The regulator matrix H̄k 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
Hk = ( λH̄k + SkY T k ) ( λI + YkY T k )−1 , (14)
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
p̂k = −Hkĝk. (15) 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 Hk. To describe this, note that the dimensions of the matrices involved are
Hk ∈ Rd×d, Yk ∈ Rd×m, Sk ∈ Rd×m. (16) We can employ the Sherman–Morrison–Woodbury formula to arrive at the following equivalent expression for Hk
Hk = ( H̄k + λ −1SkY T k ) [ I − Yk ( λI + Y Tk Yk )−1 Y Tk ] . (17)
Importantly, the matrix inverse ( λI + Y Tk Yk )−1 is now by construction a positive definite matrix of size m×m. Therefore, we construct and maintain a Cholesky factor of λI + Y Tk Yk since this leads to efficient solutions. In particular, if we express this matrix via a Cholesky decomposition
RTkRk = λI + Y T k Yk, (18)
where Rk ∈ Rm×m is an upper triangular matrix, then the search direction p̂k = −Hkĝk can be computed via
p̂k = −H̄kzk − λ−1Sk(Y Tk zk), zk = ĝk − Ykwk, wk = R−1k ( R−Tk ( Y Tk ĝk )) . (19)
Note that the above expressions for p̂k involve first computing wk, which itself involves a computationally efficient forward-backward substitution (recalling that Rk is an m ×m upper triangular matrix and the memory length m is typically 10–50). Furthermore, H̄k is typically diagonal so that H̄kzk 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 Rk can be achieved in several ways. The so-called normal-equation method constructs the (upper triangular) part of λI + Y Tk Yk and then employs a Cholesky routine, which produces Rk inO(nm(m+1)2 +m 3/3) operations. Alternatively, we can computeRk by applying Givens rotations
or Householder reflections to the matrixMk = [√ λI Y Tk ]T . This costsO(2m2((n+m)−m/3)) operations, and is therefore more expensive, but typically offers better numerical accuracy (Golub & Van Loan, 2012).

heading: 4.2 FAST AND ROBUST INCLUSION OF NEW MEASUREMENTS
text: In order to maximise the speed, we have developed a method for updating a Cholesky factor given the new measurement pair (sk+1, ŷk+1). Suppose we start with a Cholesky factor Rk at iteration k such that
RTkRk = λI + Y T k Yk. (20)
Assume, without loss of generality, that Yk and Sk are ordered in the following manner
Yk , [Y1, ŷk−m+1,Y2] , Sk , [S1, sk−m+1,S2] , (21) where Y1, Y2, S1 and S2 are defined as
Y1 , [ŷk−m+`+1, . . . , ŷk] , Y2 , [ŷk−m+2, . . . , ŷk−m+`] , (22a) S1 , [sk−m+`+1, . . . , sk] , S2 , [sk−m+2, . . . , sk−m+`] , (22b)
and ` is an appropriate integer so that Yk and Sk have m columns. The above ordering arises from “wrapping-around” the index when storing the measurements. We create the new Yk+1 and Sk+1 by replacing the oldest column entries, ŷk−m+1 and sk−m+1, with the latest measurements ŷk+1 and sk+1, respectively, so that
Yk+1 , [Y1, ŷk+1,Y2] , Sk+1 , [S1, sk+1,S2] . (23)
The aim is to generate a new Cholesky factor Rk+1 such that
RTk+1Rk+1 = λI + Y T k+1Yk+1. (24)
To this end, let the upper triangular matrix Rk be written conformally with the columns of Yk as
Rk = [R1 r1 R2 · r2 r3 · · R4 ] , (25)
so thatR1 andR2 have the same number of columns as Y1 and Y2, respectively. Furthermore, r1 is a column vector, r2 is a scalar and r3 is a row vector. Therefore,
RTkRk = RT1R1 RT1 r1 RT1R2· r22 + rT1 r1 rT1R2 + r2r3 · · RT4R4 +RT2R2 + rT3 r3  =
λI + YT1 Y1 YT1 ŷk−m+1 YT1 Y2· λ+ ŷTk−m+1ŷk−m+1 ŷTk−m+1Y2 · · λI + YT2 Y2  . (26) By observing a common structure for the update λI + Y Tk+1Yk+1 it is possible to write
λI + Y Tk+1Yk+1 = λI + YT1 Y1 YT1 ŷk+1 YT1 Y2· λ+ ŷTk+1ŷk−m+1 ŷTk+1Y2 · · λI + YT2 Y2  =
RT1R1 RT1 r4 RT1R2· r25 + rT4 r4 rT4R2 + r5r6 · · RT6R6 +RT2R2 + rT6 r6  , (27) where r4, r5 and r6 are determined by
r4 = R−T1 (YT1 ŷk+1), r5 = ( λ+ ŷTk+1ŷk+1 − rT4 r4 )1/2 , r6 = 1
r5
( ŷTk+1Y2 − rT4R2 ) . (28)
The final termR6 can be obtained by noticing that
RT6R6 +RT2R2 + rT6 r6 = RT4R4 +RT2R2 + rT3 r3, (29) which implies
RT6R6 = RT4R4 − rT6 r6 + rT3 r3. (30) ThereforeR6 can be obtained in a computationally very efficient manner by down-dating and updating the Cholesky factor R4 with the rank-1 matrices rT6 r6 and rT3 r3, respectively (see e.g. Section 12.5.3 in Golub & Van Loan (2012)).
4.3 SELECTING H̄k
There is no magic way of selecting the prior matrix H̄k. However, in practise it has proved very useful to employ a simple strategy of H̄k , γkI, where the positive scalar γk > 0 is adaptively chosen in each iteration. As a crude measure of progress we adopt the following rule
γk = { κγk−1, if αk−1 = 1, γk−1/κ, if αk−1 < 1/ρq, γk−1, otherwise,
(31)
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 Hk. Since p̂k = −Hĝk, this limits the magnitude of p̂k and the change ‖xk+1− xk‖ 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 H̄k comes with an efficiency benefit, since the product H̄kzk in equation 19 then is obtained as the element-wise product between two vectors.

heading: 4.4 ENSURING A DESCENT DIRECTION
text: In deterministic quasi-Newton methods, the search direction pk must be chosen to ensure a descent direction such that pTkgk < 0, since this guarantees reduction in the cost function for sufficiently small step lengths αk. Since pk = −Hgk, we have that pTkgk = −gTkHkgk which is always negative if the approximation Hk of the inverse Hessian is positive definite. Otherwise, we can modify the search direction by subtracting a multiple of the gradient pk ← pk − βkgk. This is motivated by noticing that
(pk − βgk)Tgk = pTkgk − βkgTk gk, (32)
which always can be made negative by selecting βk large enough, i.e. if
βk > pTkgk gTk gk . (33)
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.
lemma 3 If
βk > pTkgk − σ2g Tr (H) gTk gk + dσ 2 g , (34)
then
E [ (p̂k − βkĝk)Tĝk ] < 0. (35)
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 gk and pk, nor the noise variance σ2g . Instead we suggest a pragmatic approach in which these quantities
are replaced by their estimates ĝk, p̂k and σ̂2g . The noise estimate σ̂2g 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.

heading: 5 NUMERICAL EXPERIMENTS
text: 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.

heading: 5.1 LOGISTIC LOSS AND A 2-NORM REGULARISER
text: 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).
1yann.lecun.com/exdb/mnist/
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) (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.

heading: 5.2 DEEP LEARNING
text: 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.

heading: 6 CONCLUSION AND FUTURE WORK
text: 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.

heading: A PROOFS OF THE LEMMAS
text: 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. For 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 xk. This means that the noise contaminating the measurement corresponds to the realisation of z.
Consider the stochastic version of the Wolfe condition Ez,z′ [ f̂z′(x+ αp̂z)− f̂z(x)− cαp̂Tz ĝz ] ≤ 0, (36)
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
f̂z′(x+ αp̂z) ≈ f̂z′ + αp̂Tz ĝz′ . (37)
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
Ez,z′ [ f̂z′ − f̂z − αp̂Tz ĝz′ + cαp̂Tz ĝz ] ≤ 0. (38)
Let b denote a possible bias in the cost function estimator, then we have that Ez′ [ f̂z′ ] = f + b and
Ez [ f̂z ] = f + b. Hence, the first two terms in equation 38 cancel and we have
c ≤ Ez,z′
[ p̂Tz ĝz′ ] Ez,z′ [p̂Tz ĝz] = Ez [ p̂Tz ] Ez′ [ĝz′ ] Ez [p̂Tz ĝz] =
pTg
Ez [p̂Tz ĝz] . (39)
The denominator can be written as Ez [ p̂Tz ĝz ] = Ez [ −ĝTzHĝz ] = −Ez [ Tr ( Hĝz ĝ T z )] = −Tr ( H Ez [ ĝz ĝ T z ]) . (40)
Using the definition of covariance we have that Ez [ ĝz ĝ T z ] = Ez [ĝz] Ez [ ĝTz ] + Covĝ = ggT + σ2gI, (41)
where we in the last equality made use of the assumption that the gradients are unbiased and that Covĝ = σgI . Summarising we have that
Ez [ p̂Tz ĝz ] = −Tr ( HggT ) − σ2g Tr (H) = pTg − σ2g Tr (H) , (42)
and hence we get
c ≤ p Tg
pTg − σ2g Tr (H) . (43)
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 Ez′ [ f̂z′(x+ αp) ] < Ez [ f̂z(x) ] . (44)
We do this by studying the difference Ez′ [ f̂z′(x+ αp) ] − Ez [ f̂z(x) ] = Ez,z′ [ f̂z′(x+ αp)− f̂z(x) ] , (45)
in which the first term can be expressed using a Taylor series according to
f̂z′(x+ αp̂z) = f̂z′ + αp̂ T z ĝz′ +O(α‖p̂z‖2). (46)
For small values of α we can discard the termO(α‖p̂z‖2) from this expression. Hence we have that Ez,z′ [ f̂z′ + αp̂ T z ĝz′ − f̂z ] = αEz,z′ [ p̂Tz ĝz′ ] = αEz [ p̂Tz ] Ez′ [ĝz′ ] = αp Tg. (47)
Following from Section 4 we can ensure pTg 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 Ez [ (p̂z − βĝz)Tĝz ] = Ez [ p̂Tz ĝz ] − β Ez [ ĝTz ĝz ] = Ez [ −ĝTzHĝz ] − β Ez [ ĝTz ĝz ] . (48)
Using equation 42 we have Ez [ −ĝTzHĝz ] = −Tr ( HggT ) − σ2g Tr (H) = pTg − σ2g Tr (H) , (49)
and it directly follows from this that Ez [ ĝTz ĝz ] = Ez [ −ĝTz Iĝz ] = gTg + dσ2g . (50)
Hence
Ez [ (p̂z − βĝz)Tĝz ] = pTg − σ2g Tr (H)− β ( gTg + dσ2g ) , (51)
and we see that
Ez [ (p̂z − βĝz)Tĝz ] < 0⇔ β >
pTg − σ2g Tr (H) gTg + dσ2g . (52)

heading: B QUASI-NEWTON METHODS
text: The historically most popular quasi-Newton method is given by the BFGS algorithm (Broyden, 1970; Fletcher, 1970; Goldfarb, 1970; Shanno, 1970)
Hk+1 =
( I − sky T k
yTk sk
) Hk ( I − yks T k
yTk sk
) + sks T k
yTk sk . (53)
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
Hk+1 = Hk − Hkyky
T kHk
yTkHkyk + sks
T k
yTk sk . (54)
Both of these are rank 2 updates that ensure Hk 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
Hk = Hk−1 + (sk −Hk−1yk)(sk −Hk−1yk)T
(sk −Hk−1yk)Tyk . (55)
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 Hk 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 pk directly without the need of storing the full matrix Hk. 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).

heading: C RESULTING ALGORITHM
text: 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.

heading: D EXPERIMENT DETAILS
text: 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 H̄k did not have much impact in most of the logistic regression examples, and thus it was kept constant H̄k = γ0I except for in the covtype and URL problems. In the CIFAR and MNIST problems, however, the procedure was found to be of significant importance.
Algorithm 2 Stochastic quasi-Newton with line search Require: An initial estimate x1, a maximum number of iterations kmax, memory limit m, regular-
isation parameter λ, scale factor ρ ∈ (0, 1), reduction limit ξ ≥ 1, backtracking limit τ > 0
1: Set k = 1 2: while k < kmax do 3: Obtain a measurement of the cost function and its gradient
f̂k = f(xk) + ek,
ĝk = gk + vk.
4: Set H̄k = γkI where
γk = { κγk−1, if αk−1 = 1, γk−1/κ, if αk−1 < 1/ρq, γk−1, otherwise.
5: if ŷTk sk > ‖sk‖22 then 6: if k > m then 7: Form Yk and Sk by replacing the oldest vector-pairs in Yk−1 and Sk−1 with ŷk and sk 8: else 9: Form Yk and Sk by adding ŷk and sk to Yk−1 and Sk−1
10: end if 11: else 12: Set Yk = Yk−1 and Sk = Sk−1 13: end if 14: Select p̂k as
p̂k = −H̄kzk − λ−1Sk(Y Tk zk), zk = ĝk − Ykwk, wk = R −1 k ( R−Tk ( Y Tk ĝk )) ,
with details provided in Section 4.2. 15: Set p̂k ← p̂k − βkĝk with βk where
βk > p̂Tk ĝk − σ̂2g Tr (H)
ĝTk ĝk + dσ̂ 2 g
.
16: Set αk = min{1, ξ/k} 17: Set i = 1 18: while f̂(xk + αp̂k) > f̂(xk) + cαkĝTk p̂k and i ≤ max{0, τ − k} do 19: Reduce the step length αk ← ραk 20: Set i← i+ 1 21: end while 22: Update xk+1 = xk + αkp̂k 23: Set k ← k + 1 24: end while

Review summaries:
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.