Global Convergence of Online Limited Memory BFGS

Introduction

SGD is the most popular method used to solve stochastic optimization problems (Bottou (2010); Shalev-Shwartz et al. (2007); Zhang (2004); LeRoux et al. (2012)). However, as we consider problems of ever larger dimension their slow convergence times have limited their practical appeal and fostered the search for alternatives. In this regard, it has to be noted that SGD is slow because of both, the use of gradients as descent directions and their replacement by random estimates. Several alternatives have been proposed to deal with randomness in an effort to render the convergence times of SGD closer to the faster convergence times of gradient descent (Syski (1983); Konecny and Richtarik (2013); Zhang et al. (2013a)). These SGD variants succeed in reducing randomness and end up exhibiting the asymptotic convergence rate of gradient descent. Although they improve asymptotic convergence rates, the latter methods are still often slow in practice. This is not unexpected. Reducing randomness is of no use when the function $F({\mathbf{w}})$ has a challenging curvature profile. In these ill-conditioned functions SGD is limited by the already slow convergence times of deterministic gradient descent. The golden standard to deal with ill-conditioned functions in a deterministic setting is Newton’s method. However, unbiased stochastic estimates of Newton steps can’t be computed in general. This fact limits the application of stochastic Newton methods to problems with specific structure (Birge et al. (1995); Zargham et al. (2013)).

If SGD is slow to converge and stochastic Newton can’t be used in general, the remaining alternative is to modify deterministic quasi-Newton methods that speed up convergence times relative to gradient descent without using Hessian evaluations (J. E. Dennis and More (1974); Powell (1971); Byrd et al. (1987); Nocedal and Wright (1999)). This has resulted in the development of the stochastic quasi-Newton methods known as online (o) Broyden-Fletcher-Goldfarb-Shanno (BFGS) (Schraudolph et al. (2007)), regularized stochastic BFGS (RES) (Mokhtari and Ribeiro (2014a)), and online limited memory (oL)BFGS (Schraudolph et al. (2007)) which occupy the middle ground of broad applicability irrespective of problem structure and conditioning. All three of these algorithms extend BFGS by using stochastic gradients both as descent directions and constituents of Hessian estimates. The oBFGS algorithm is a direct generalization of BFGS that uses stochastic gradients in lieu of deterministic gradients. RES differs in that it further modifies BFGS to yield an algorithm that retains its convergence advantages while improving theoretical convergence guarantees and numerical behavior. The oLBFGS method uses a modification of BFGS to reduce the computational cost of each iteration.

An important observation here is that in trying to adapt to the changing curvature of the objective, stochastic quasi-Newton methods may end up exacerbating the problem. Indeed, since Hessian estimates are stochastic, it is possible to end up with almost singular Hessian estimates. The corresponding small eigenvalues then result in a catastrophic amplification of the noise which nullifies progress made towards convergence. This is not a minor problem. In oBFGS this possibility precludes convergence analyses (Bordes et al. (2009); Schraudolph et al. (2007)) and may result in erratic numerical behavior (Mokhtari and Ribeiro (2014a)). As a matter of fact, the main motivation for the introduction of RES is to avoid this catastrophic noise amplification so as to retain smaller convergence times while ensuring that optimal arguments are found with probability 1 (Mokhtari and Ribeiro (2014a)). However valuable, the convergence guarantees of RES and the convergence time advantages of oBFGS and RES are tainted by an iteration cost of order $O(n^{2})$ and $O(n^{3})$, respectively, which precludes their use in problems where $n$ is very large. In deterministic settings this problem is addressed by limited memory (L)BFGS (Dong C. and Nocedal (1989)) which can be easily generalized to develop the oLBFGS algorithm (Schraudolph et al. (2007)). Numerical tests of oLBFGS are promising but theoretical convergence characterizations are still lacking. The main contribution of this paper is to show that oLBFGS converges with probability 1 to optimal arguments across realizations of the random variables $\theta$. This is the same convergence guarantee provided for RES and is in marked contrast with oBFGS, which fails to converge if not properly regularized. Convergence guarantees for oLBFGS do not require such measures.

We begin the paper with brief discussions of deterministic BFGS (Section 2) and LBFGS (Section 2.1) and the introduction of oLBFGS (Section 2.2). The fundamental idea in BFGS and oLBFGS is to continuously satisfy a secant condition while staying close to previous curvature estimates. They differ in that BFGS uses all past gradients to estimate curvature while oLBFGS uses a fixed moving window of past gradients. The use of this window reduces memory and computational cost (Appendix A). The difference between LBFGS and oLBFGS is the use of stochastic gradients in lieu of their deterministic counterparts.

Convergence properties of oLBFGS are then analyzed (Section 3). Under the assumption that the sample functions $f({\mathbf{w}},{\boldsymbol{\theta}})$ are strongly convex we show that the trace and determinant of the Hessian approximations computed by oLBFGS are upper and lower bounded, respectively (Lemma 3). These bounds are then used to limit the range of variation of the ratio between the Hessian approximations’ largest and smallest eigenvalues (Lemma 4). In turn, this condition number limit is shown to be sufficient to prove convergence to the optimal argument ${\mathbf{w}}^{*}$ with probability 1 over realizations of the sample functions (Theorem 6). This is an important result because it ensures that oLBFGS doesn’t suffer from the numerical problems that hinder oBFGS. We complement this almost sure convergence result with a characterization of the convergence rate which is shown to be at least $O(1/t)$ in expectation (Theorem 7). It is fair to emphasize that, different from the deterministic case, the convergence rate of oLBFGS is not better than the convergence rate of SGD. This is not a limitation of our analysis. The difference between stochastic and regular gradients introduces a noise term that dominates convergence once we are close to the optimum, which is where superlinear convergence rates manifest. In fact, the same convergence rate would be observed if exact Hessians were available. The best that can be proven of oLBFGS is that the convergence rate is not worse than that of SGD. Given that theoretical guarantees only state that the curvature correction does not exacerbate the problem’s condition it is perhaps fairer to describe oLBFGS as an adaptive reconditioning strategy instead of a stochastic quasi-Newton method. The latter description refers to the genesis of the algorithm. The former is a more accurate description of its actual behavior.

To show the advantage of using oLBFGS as an adaptive reconditioning strategy we develop its application to SVM problems (Section 4) and perform a comparative numerical analysis with synthetic data. The conclusions of this numerical analysis are that oLBFGS performs as well as oBFGS and RES while outperforming SGD when convergence is measured with respect to the number of feature vectors processed. In terms of computation time, oLBFGS outperforms all three methods, SGD, oBFGS, and RES. The advantages of oLBFGS grow with the dimension of the feature vector and can be made arbitrarily large (Section 4.1). To further substantiate numerical claims we use oLBFGS to train a logistic regressor to predict the click through rate in a search engine advertising problem (Section 5). The logistic regression uses a heterogeneous feature vector with 174,026 binary entries that describe the user, the search, and the advertisement (Section 5.1). Being a large scale problem with heterogeneous data, the condition number of the logistic log likelihood objective is large and we expect to see significant advantages of oLBFGS relative to SGD. This expectation is fulfilled. The oLBFGS algorithm trains the regressor using less than $1\%$ of the data required by SGD to obtain similar classification accuracy. (Section 5.3). We close the paper with concluding remarks (Section 6).

Algorithm definition

Since the function $F({\mathbf{w}})$ is strongly convex, gradients ${\mathbf{s}}({\mathbf{w}})$ are descent directions that can be used to find the optimal argument ${\mathbf{w}}^{*}$ in (1). Introduce then a time index $t$, a step size $\epsilon_{t}$, and a positive definite matrix ${\mathbf{B}}_{t}^{-1}\succ 0$ to define a generic descent algorithm through the iteration

where we have also defined the descent step ${\mathbf{d}}_{t}={\mathbf{B}}_{t}^{-1}{\mathbf{s}}({\mathbf{w}}_{t})$. When ${\mathbf{B}}_{t}^{-1}={\mathbf{I}}$ is the identity matrix, (3) reduces to gradient descent. When ${\mathbf{B}}_{t}={\mathbf{H}}({\mathbf{w}}_{t}):=\nabla^{2}F({\mathbf{w}}_{t})$ is the Hessian of the objective function, (3) defines Newton’s algorithm. In this paper we focus on quasi-Newton methods whereby we attempt to select matrices ${\mathbf{B}}_{t}$ close to the Hessian ${\mathbf{H}}({\mathbf{w}}_{t})$. Various methods are known to select matrices ${\mathbf{B}}_{t}$, including those by Broyden e.g., Broyden et al. (1973); Davidon, Fletcher, and Powell (DFP) e.g., Fletcher (2013); and Broyden, Fletcher, Goldfarb, and Shanno (BFGS) e.g., Byrd et al. (1987); Powell (1971). We work with the matrices ${\mathbf{B}}_{t}$ used in BFGS since they have been observed to work best in practice (see Byrd et al. (1987)).

In BFGS, the function’s curvature ${\mathbf{B}}_{t}$ is approximated by a finite difference. Let ${\mathbf{v}}_{t}$ denote the variable variation at time $t$ and ${\mathbf{r}}_{t}$ the gradient variation at time $t$ which are respectively defined as

We select the matrix ${\mathbf{B}}_{t+1}$ to be used in the next time step so that it satisfies the secant condition ${\mathbf{B}}_{t+1}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$. The rationale for this selection is that the Hessian ${\mathbf{H}}({\mathbf{w}}_{t})$ satisfies this condition for ${\mathbf{w}}_{t+1}$ tending to ${\mathbf{w}}_{t}$. Notice however that the secant condition ${\mathbf{B}}_{t+1}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$ is not enough to completely specify ${\mathbf{B}}_{t+1}$. To resolve this indeterminacy, matrices ${\mathbf{B}}_{t+1}$ in BFGS are also required to be as close as possible to the previous Hessian approximation ${\mathbf{B}}_{t}$ in terms of differential entropy. These conditions can be resolved in closed form leading to the explicit expression – see, e.g., Nocedal and Wright (1999) –,

While the expression in (5) permits updating the Hessian approximations ${\mathbf{B}}_{t+1}$, implementation of the descent step in (3) requires its inversion. This can be avoided by using the Sherman-Morrison formula in (5) to write

where we defined the scalar $\rho_{t}$ and the matrix ${\mathbf{Z}}_{t}$ as

The updates in (5) and (6) require the inner product of the gradient and variable variations to be positive, i.e., ${\mathbf{v}}_{t}^{T}{\mathbf{r}}_{t}>0$. This is always true if the objective $F({\mathbf{w}})$ is strongly convex and further implies that ${\mathbf{B}}_{t+1}^{-1}$ stays positive definite if ${\mathbf{B}}_{t}^{-1}\succ{\mathbf{0}}$, Nocedal and Wright (1999).

Each BFGS iteration has a cost of $O(n^{2})$ arithmetic operations. This is less than the $O(n^{3})$ of each step in Newton’s method but more than the $O(n)$ cost of each gradient descent iteration. In general, the relative convergence rates are such that the total computational cost of BFGS to achieve a target accuracy is smaller than the corresponding cost of gradient descent. Still, alternatives to reduce the computational cost of each iteration are of interest for large scale problems. Likewise, BFGS requires storage and propagation of the $O(n^{2})$ elements of ${\mathbf{B}}_{t}^{-1}$, whereas gradient descent requires storage of $O(n)$ gradient elements only. This motivates alternatives that have smaller memory footprints. Both of these objectives are accomplished by the limited memory (L)BFGS algorithm that we describe in the following section.

As it follows from (6), the updated Hessian inverse approximation ${\mathbf{B}}_{t}^{-1}$ depends on ${\mathbf{B}}_{t-1}^{-1}$ and the curvature information pairs $\{{\mathbf{v}}_{t-1}$, ${\mathbf{r}}_{t-1}\}$. In turn, to compute ${\mathbf{B}}_{t-1}^{-1}$, the estimate ${\mathbf{B}}_{t-2}^{-1}$ and the curvature pair $\{{\mathbf{v}}_{t-2},{\mathbf{r}}_{t-2}\}$ are used. Proceeding recursively, it follows that ${\mathbf{B}}_{t}^{-1}$ is a function of the initial approximation ${\mathbf{B}}_{0}^{-1}$ and all previous $t$ curvature information pairs $\{{\mathbf{v}}_{u}$, ${\mathbf{r}}_{u}\}_{u=0}^{t-1}$. The idea in LBFGS is to restrict the use of past curvature information to the last $\tau$ pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}_{u=t-\tau}^{t-1}$. Since earlier iterates $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}$ with $u<t-\tau$ are likely to carry little information about the curvature at the current iterate ${\mathbf{w}}_{t}$, this restriction is expected to result in a minimal performance penalty.

For a precise definition, pick a positive definite matrix ${\mathbf{B}}_{t,0}^{-1}$ as the initial Hessian inverse approximation at step $t$. Proceed then to perform $\tau$ updates of the form in (6) using the last $\tau$ curvature information pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}_{u=t-\tau}^{t-1}$. Denoting as ${\mathbf{B}}_{t,u}^{-1}$ the curvature approximation after $u$ updates are performed we have that the refined matrix approximation ${\mathbf{B}}_{t,u+1}^{-1}$ is given by [cf. (6)]

where $u=0,\dots,\tau-1$ and the constants $\rho_{t-\tau+u}$ and rank-one plus identity matrices ${\mathbf{Z}}_{t-\tau+u}$ are as given in (7). The inverse Hessian approximation ${\mathbf{B}}_{t}^{-1}$ to be used in (3) is the one yielded after completing the $\tau$ updates in (8), i.e., ${\mathbf{B}}_{t}^{-1}={\mathbf{B}}_{t,\tau}^{-1}$. Observe that when $t<\tau$ there are not enough pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}$ to perform $\tau$ updates. In such case we just redefine $\tau=t$ and proceed to use the $t=\tau$ available pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}_{u=0}^{t-1}$ .

Implementation of the product ${\mathbf{B}}_{t}^{-1}{\mathbf{s}}({\mathbf{w}}_{t})$ in (3) for matrices ${\mathbf{B}}_{t}^{-1}={\mathbf{B}}_{t,\tau}^{-1}$ obtained from the recursion in (8) does not need explicit computation of the matrix ${\mathbf{B}}_{t,\tau}^{-1}$. Although the details are not straightforward, observe that each iteration in (8) is similar to a rank-one update and that as such it is not unreasonable to expect that the product ${\mathbf{B}}_{t}^{-1}{\mathbf{s}}({\mathbf{w}}_{t})={\mathbf{B}}_{t,\tau}^{-1}{\mathbf{s}}({\mathbf{w}}_{t})$ can be computed using $\tau$ recursive inner products. Assuming that this is possible, the implementation of the recursion in (8) doesn’t need computation and storage of prior matrices ${\mathbf{B}}_{t-1}^{-1}$. Rather, it suffices to keep the $\tau$ most recent curvature information pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}_{u=t-\tau}^{t-1}$, thus reducing storage requirements from $O(n^{2})$ to $O(\tau n)$. Furthermore, each of these inner products can be computed at a cost of $n$ operations yielding a total computational cost of $O(\tau n)$ per LBFGS iteration. Hence, LBFGS decreases both the memory requirements and the computational cost of each iteration from the $O(n^{2})$ required by regular BFGS to $O(\tau n)$. We present the details of this iteration in the context of the online (stochastic) LBFGS that we introduce in the following section.

2 Online (Stochastic) Limited memory BFGS

To define the oLBFGS algorithm we just need to provide stochastic versions of the definitions in (7) and (8). The scalar constants and identity plus rank-one matrices in (7) are redefined to the corresponding stochastic quantities

whereas the LBFGS matrix ${\mathbf{B}}_{t}^{-1}={\mathbf{B}}_{t,\tau}^{-1}$ in (8) is replaced by the oLBFGS Hessian inverse approximation ${\hat{\mathbf{B}}}_{t}^{-1}={\hat{\mathbf{B}}}_{t,\tau}^{-1}$ which we define as the outcome of $\tau$ recursive applications of the update,

where the initial matrix ${\hat{\mathbf{B}}}_{t,0}^{-1}$ is given and the time index is $u=0,\ldots,\tau-1$. The oLBFGS algorithm is defined by the stochastic descent iteration in (10) with matrices ${\hat{\mathbf{B}}}_{t}^{-1}={\hat{\mathbf{B}}}_{t,\tau}^{-1}$ computed by $\tau$ recursive applications of (13). Except for the fact that they use stochastic variables, (10) and (13) are identical to (3) and (8). Thus, as is the case in (3), the Hessian inverse approximation ${\hat{\mathbf{B}}}_{t}^{-1}$ in (13) is a function of the initial Hessian inverse approximation ${\mathbf{B}}_{t,0}^{-1}$ and the $\tau$ most recent curvature information pairs $\{{\mathbf{v}}_{u},{\hat{\mathbf{r}}}_{u}\}_{u=t-\tau}^{t-1}$. Likewise, when $t<\tau$ there are not enough pairs $\{{\mathbf{v}}_{u},{\hat{\mathbf{r}}}_{u}\}$ to perform $\tau$ updates. In such case we just redefine $\tau=t$ and proceed to use the $t=\tau$ available pairs $\{{\mathbf{v}}_{u},{\hat{\mathbf{r}}}_{u}\}_{u=0}^{t-1}$ . We also point out that the update in (13) necessitates ${\hat{\mathbf{r}}}_{u}^{T}{\mathbf{v}}_{u}>0$ for all time indexes $u$. This is true as long as the instantaneous functions $f({\mathbf{w}},{\boldsymbol{\theta}})$ are strongly convex with respect to ${\mathbf{w}}$ as we show in Lemma 2.

Equation (14) shows the relation between the Hessian inverse approximation ${\hat{\mathbf{B}}}_{t}^{-1}$ and the $(\tau-1)$st updated version of the initial Hessian inverse approximation ${\hat{\mathbf{B}}}_{t,\tau-1}^{-1}$ at step $t$. Set now $u=\tau-2$ in (13) to express ${\hat{\mathbf{B}}}_{t,\tau-1}^{-1}$ in terms of ${\hat{\mathbf{B}}}_{t,\tau-2}^{-1}$ and substitute the result in (14) to rewrite ${\hat{\mathbf{B}}}_{t}^{-1}$ as

We can proceed recursively by substituting ${\hat{\mathbf{B}}}_{t,\tau-2}^{-1}$ for its expression in terms of ${\hat{\mathbf{B}}}_{t,\tau-3}^{-1}$ and in the result substitute ${\hat{\mathbf{B}}}_{t,\tau-3}^{-1}$ for its expression in terms of ${\hat{\mathbf{B}}}_{t,\tau-3}^{-1}$ and so on. Observe that a new summand is added in each of these substitutions from which it follows that repeating this process $\tau$ times yields

Consider the oLBFGS Hessian inverse approximation ${\hat{\mathbf{B}}}_{t}^{-1}={\hat{\mathbf{B}}}_{t,\tau}^{-1}$ obtained after $\tau$ recursive applications of the update in (13) with the scalar sequence $\hat{\rho}_{t-\tau+u}$ and identity plus rank-one matrix sequence ${\hat{\mathbf{Z}}}_{t-\tau+u}$ as defined in (12) for given variable and stochastic gradient variation pairs $\{{\mathbf{v}}_{u},{\mathbf{r}}_{u}\}_{u=t-\tau}^{t-1}$. For a given vector ${\mathbf{p}}={\mathbf{p}}_{0}$ define the sequence of vectors ${\mathbf{p}}_{k}$ through the recursion

where we also define the constants $\alpha_{u}:=\hat{\rho}_{t-u-1}{\mathbf{v}}_{t-u-1}^{T}{\mathbf{p}}_{u}$. Further define the sequence of vectors ${\mathbf{q}}_{k}$ with initial value ${\mathbf{q}}_{0}={\hat{\mathbf{B}}}_{t,0}^{-1}{\mathbf{p}}_{\tau}$ and subsequent elements

where we define constants $\beta_{u}:=\hat{\rho}_{t-\tau+u}{\hat{\mathbf{r}}}_{t-\tau+u}^{T}{\mathbf{q}}_{u}$. The product ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}$ equals ${\mathbf{q}}_{\tau}$, i.e., ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}\ =\ {\mathbf{q}}_{\tau}$.

Proof See Appendix A. Proposition 1 asserts that it is possible to reduce the computation of the product ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}$ between the oLBFGS Hessian approximation matrix and arbitrary vector ${\mathbf{p}}$ to the computation of two vector sequences $\{{\mathbf{p}}_{u}\}_{u=0}^{\tau-1}$ and $\{{\mathbf{q}}_{u}\}_{u=0}^{\tau-1}$. The product ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}\ =\ {\mathbf{q}}_{\tau}$ is given by the last element of the latter sequence. Since determination of each of the elements of each sequence requires $O(n)$ operations and the total number of elements in each sequence is $\tau$ the total operation cost to compute both sequences is of order $O(\tau n)$. In computing ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}$ we also need to add the cost of the product ${\mathbf{q}}_{0}={\hat{\mathbf{B}}}_{t,0}^{-1}{\mathbf{p}}_{\tau}$ that links both sequences. To maintain overall computation cost of order $O(\tau n)$ this matrix has to have a sparse or low rank structure. A common choice in LBFGS, that we adopt for oLBFGS, is to make ${\hat{\mathbf{B}}}_{t,0}^{-1}=\hat{\gamma}_{t}{\mathbf{I}}$. The scalar constant $\hat{\gamma}_{t}$ is a function of the variable and stochastic gradient variations ${\mathbf{v}}_{t-1}$ and ${\hat{\mathbf{r}}}_{t-1}$, explicitly given by

with the value at the first iteration being $\hat{\gamma}_{0}=1$. The scaling factor $\hat{\gamma}_{t}$ attempts to estimate one of the eigenvalues of the Hessian matrix at step $t$ and has been observed to work well in practice; see e.g., Dong C. and Nocedal (1989); Nocedal and Wright (1999). Further observe that the cost of computing $\hat{\gamma}_{t}$ is of order $O(n)$ and that since ${\hat{\mathbf{B}}}_{t,0}^{-1}$ is diagonal cost of computing the product ${\mathbf{q}}_{0}={\hat{\mathbf{B}}}_{t,0}^{-1}{\mathbf{p}}_{\tau}$ is also of order $O(n)$. We adopt the initialization in (19) in our subsequent analysis and numerical experiments.

Convergence analysis

Our goal here is to show that as time progresses the sequence of variable iterates ${\mathbf{w}}_{t}$ approaches the optimal argument ${\mathbf{w}}^{*}$. In proving this result we make the following assumptions.

The second moment of the norm of the stochastic gradient is bounded for all ${\mathbf{w}}$. i.e., there exists a constant $S^{2}$ such that for all variables ${\mathbf{w}}$ it holds

The step size sequence is selected as nonsummable but square summable, i.e.,

The bounds in (25) are customary in convergence proofs of descent methods. For the results here the stronger condition spelled in Assumption 1 is needed. This assumption in necessary to guarantee that the inner product ${\hat{\mathbf{r}}}_{t}^{T}{\mathbf{v}}_{t}>0$ is positive as we show in the following lemma.

Furthermore, the ratio of stochastic gradient variation squared norm $\|{\hat{\mathbf{r}}}_{t}\|^{2}={\hat{\mathbf{r}}}_{t}^{T}{\hat{\mathbf{r}}}_{t}$ to inner product of variable and stochastic gradient variations is bounded as

The analysis is easier if we consider the matrix ${\hat{\mathbf{B}}}_{t}$ – as opposed to ${\hat{\mathbf{B}}}_{t}^{-1}$. Consider then the update in (13), and use the Sherman-Morrison formula to rewrite as an update that relates ${\hat{\mathbf{B}}}_{t,u+1}$ to ${\hat{\mathbf{B}}}_{t,u}$,

for $u=0,\dots,\tau-1$ and ${\hat{\mathbf{B}}}_{t,0}=1/\hat{\gamma}_{t}{\mathbf{I}}$ as per (19). As in (13), the Hessian approximation at step $t$ is ${\hat{\mathbf{B}}}_{t}={\hat{\mathbf{B}}}_{t,\tau}$. In the following lemma we use the update formula in (28) to find bounds on the trace and determinant of the Hessian approximation ${\hat{\mathbf{B}}}_{t}$.

Consider the Hessian approximation ${\hat{\mathbf{B}}}_{t}={\hat{\mathbf{B}}}_{t,\tau}$ defined by the recursion in (28) with ${\hat{\mathbf{B}}}_{t,0}=\hat{\gamma}_{t}^{-1}{\mathbf{I}}$ and $\hat{\gamma}_{t}$ as given by (19). If Assumption 1 holds true, the trace $\text{tr}({\hat{\mathbf{B}}}_{t})$ of the Hessian approximation ${\hat{\mathbf{B}}}_{t}$ is uniformly upper bounded for all times $t\geq 1$,

Likewise, if Assumption 1 holds true, the determinant $\det({\hat{\mathbf{B}}}_{t})$ of the Hessian approximation ${\hat{\mathbf{B}}}_{t}$ is uniformly lower bounded for all times $t$

Lemma 3 states that the trace and determinants of the Hessian approximation matrix ${\hat{\mathbf{B}}}_{t}={\hat{\mathbf{B}}}_{t,\tau}$ are bounded for all times $t\geq 1$. For time $t=0$ we can write a similar bound that takes into account the fact that the constant $\gamma_{t}$ that initializes the recursion in (28) is $\gamma_{0}=1$. Given that we are interested in an asymptotic convergence analysis, this bound in inconsequential. The bounds on the trace and determinant of ${\hat{\mathbf{B}}}_{t}$ are respectivey equivalent to bounds in the sum and product of its eigenvalues. Further considering that the matrix ${\hat{\mathbf{B}}}_{t}$ is positive definite, as it follows from Lemma 2, these bounds can be further transformed into bounds on the smalls and largest eigenvalue of ${\hat{\mathbf{B}}}_{t}$. The resulting bounds are formally stated in the following lemma.

The bounds in Lemma 4 imply that their respective inverses are bounds on the range of the eigenvalues of the Hessian inverse approximation matrix ${\hat{\mathbf{B}}}_{t}^{-1}$. Specifically, the minimum eigenvalue of the Hessian inverse approximation ${\hat{\mathbf{B}}}_{t}^{-1}$ is larger than $1/C$ and the maximum eigenvalue of ${\hat{\mathbf{B}}}_{t}^{-1}$ does not exceed $1/c$, or, equivalently,

Consider the online Limited memory BFGS algorithm as defined by the descent iteration in (10) with matrices ${\hat{\mathbf{B}}}_{t}^{-1}={\hat{\mathbf{B}}}_{t,\tau}^{-1}$ obtained after $\tau$ recursive applications of the update in (13) initialized with ${\hat{\mathbf{B}}}_{t,0}^{-1}=\hat{\gamma}_{t}{\mathbf{I}}$ and $\hat{\gamma}_{t}$ as given by (19). If Assumptions 1 and 2 hold true, the sequence of average function values $F({\mathbf{w}}_{t})$ satisfies

Consider the online Limited memory BFGS algorithm as defined by the descent iteration in (10) with matrices ${\hat{\mathbf{B}}}_{t}^{-1}={\hat{\mathbf{B}}}_{t,\tau}^{-1}$ obtained after $\tau$ recursive applications of the update in (13) initialized with ${\hat{\mathbf{B}}}_{t,0}^{-1}=\hat{\gamma}_{t}{\mathbf{I}}$ and $\hat{\gamma}_{t}$ as given by (19). If Assumptions 1-3 hold true the limit infimum of the squared Euclidean distance to optimality $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|^{2}$ converges to zero almost surely, i.e.,

Theorem 7 shows that under specified assumptions the expected error in terms of the objective value after $t$ oLBFGS iterations is of order $O(1/t)$. As is the case of Theorem 6, this result is not better than the convergence rate of conventional SGD. As can be seen in the proof of Theorem 7, the convergence rate is dominated by the noise term introduced by the difference between stochastic and regular gradients. This noise term would be present even if exact Hessians were available and in that sense the best that can be proven of oLBFGS is that the convergence rate is not worse than that of SGD. Given that theorems 6 and 7 parallel the theoretical guarantees of SGD it is perhaps fairer to describe oLBFGS as an adaptive reconditioning strategy instead of a stochastic quasi-Newton method. The latter description refers to the genesis of the algorithm, but the former is more accurate description of its behavior. Do notice that while the convergence rate doesn’t change, improvements in convergence time are significant as we illustrate with the numerical experiments that we present in the next two sections.

Support vector machines

where we have also added the regularization term ${\lambda}\|{\mathbf{w}}\|^{2}/{2}$ for some constant $\lambda>0$. Common selections for the loss function are the hinge loss $l(({\mathbf{x}},y);{\mathbf{w}})=\max(0,1-y({\mathbf{w}}^{T}{\mathbf{x}}))$ and the squared hinge loss $l(({\mathbf{x}},y);{\mathbf{w}})=\max(0,1-y({\mathbf{w}}^{T}{\mathbf{x}}))^{2}$. See, e.g., Bottou (2010). To model (37) as a problem in the form of (1), define $\boldsymbol{\theta}_{i}=({\mathbf{x}}_{i},y_{i})$ as a given training point and the probability distribution of $\theta$ as uniform on the training set ${\mathcal{S}}=\{({\mathbf{x}}_{i},y_{i})\}_{i=1}^{N}=\{\boldsymbol{\theta}_{i}\}_{i=1}^{N}$. It then suffices to define

Figures 1 and 2 show the empirical distributions of the objective function value $F({\mathbf{w}}_{t})$ attained after processing $Lt=4\times 10^{4}$ feature vectors using $J=10^{3}$ realizations for the cases that $n=10^{2}$ and $n=10^{3}$, respectively. According to Figure 1 the averages of objective value function for oLBFGS, oBFGS and RES are $1.7\times 10^{-5}$, $1.4\times 10^{-5}$ and $1.9\times 10^{-5}$, respectively. These numbers show that the performance of oLBFGS is very close to the performances of oBFGS and RES. This similarity holds despite the fact that oLBFGS uses only the last $\tau=10$ stochastic gradients to estimate curvature whereas oBFGS and RES utilize all past stochastic gradients to do so. The advantage of oLBFGS is in the smaller computational cost of processing feature vectors as we discuss in Section 4.2. The corresponding average objective values achieved by SGD and SAG after processing $Lt=4\times 10^{4}$ feature vectors are $1.6\times 10^{-3}$ and $5.7\times 10^{-4}$, respectively. Both of these are at least an order of magnitude larger than the average objective value achieved by oLBFGS – or RES and oBFGS for that matter.

Figure 2 repeats the study in Figure 1 for the case in which the feature vector dimension is increased to $n=10^{3}$. The performance of oLBGS is still about the same as the performances of oBFGS and RES. The average objective function values achieved after processing $Lt=4\times 10^{4}$ feature vectors are $9.9\times 10^{-6}$, $9.8\times 10^{-6}$ and $9.5\times 10^{-6}$ for oLBFGS, oBFGS and RES, respectively. The relative performance with respect to SGD and SAG, however, is now larger. The averages of objective function values for SAG and SGD in this case are $2.1\times 10^{-2}$ and $4.5\times 10^{-2}$, respectively. These values are more than 3 orders of magnitude larger than the corresponding values achieved by oLBFGS. This relative improvement can be further increased if we consider problems of even larger dimension. Further observe that oBFGS and RES start to become impractical if we further increase the feature vector dimension since the respective iterations have computational costs of order $O(n^{2})$ and $O(n^{3})$. We analyze this in detail in the following section.

2 Convergence versus processing time

The analysis in Section 4.1 is relevant for online implementations in which the goal is to make the best possible use of the information provided by each new acquired feature vector. In implementations where computational cost is of dominant interest we have to account for the fact that the respective iteration costs are of order $O(n)$ for SGD and SAG, of order $O(\tau n)$ for oLBFGS, and of orders $O(n^{2})$ and $O(n^{3})$ for oBFGS and RES. As we increase the problem dimension we expect the convergence time advantages of oBFGS and RES in terms of number of feature vectors processed to be overwhelmed by the increased computational cost of each iteration. For oLBFGS, on the contrary, we expect the convergence time advantages in terms of number of feature vectors processed to persist in terms of processing time. To demonstrate that this is the case we repeat the experiments in Section 4.1 but record the processing time required to achieve a target objective value. The parameters used here are the same parameters of Section 4.1.

In Figure 3 we consider $n=10^{2}$ and record the processing time required to achieve the objective function value $F({\mathbf{w}}_{t})=10^{-4}$. Histograms representing empirical distributions of execution times measured in seconds (s) are shown for oLBFGS, oBFGS, RES, SGD, and SAG. We also summarize the average minimum and maximum times observed for each algorithm. The average run times for oBFGS and RES are $0.14\,\text{s}$ and $0.26\,\text{s}$ which are better than the average run times of SGD and SAG that stand at $0.63\,\text{s}$ and $0.50\,\text{s}$. The advantage, however, is less marked than when measured with respect to the number of feature vector processed. For oLBGS the advantage with respect to SGD and SAG is still close to one order of magnitude since the average convergence time stands at $0.073\,\text{s}$. When measured in computation time oLBGS is also better than RES and oBFGS, as expected.

Figure 4 presents the analogous histograms and summary statistics when the feature vector dimension is $n=10^{3}$ and the algorithm is run until achieving the objective value $F({\mathbf{w}}_{t})=10^{-5}$. For this problem and metric the performances of RES and oBFGS are worse than the corresponding performances of SGD and SAG. The respective average convergence times are $7.7\,\text{s}$ and $4.1\,\text{s}$ for RES and oBFGS and $1.4\,\text{s}$ and $2.0\,\text{s}$ for SAG and SGD. The oLBFGS algorithm, however, has an average convergence time of $0.11\,\text{s}$. This is still an order of magnitude faster than the first order methods SAG and SGD – and has an even larger advantage with respect to oBFGS and RES, by extension. The relative reduction of execution times of oLBGS relative to all other 4 methods becomes more marked for problems of larger dimension. We investigate these advantages on the search engine advertising problem that we introduce in the following section.

Search engine advertising

We apply oLBFGS to the problem of predicting the click-through rate (CTR) of an advertisement displayed in response to a specific search engine query by a specific visitor. In these problems we are given meta information about an advertisement, the words that appear in the query, as well as some information about the visitor and are asked to predict the likelihood that this particular ad is clicked by this particular user when performing this particular query. The information specific to the ad includes descriptors of different characteristics such as the words that appear in the title, the name of the advertiser, keywords that identify the product, and the position on the page where the ad is to be displayed. The information specific to the user is also heterogeneous and includes gender, age, and propensity to click on ads. To train a classifier we are given information about past queries along with the corresponding click success of the ads displayed in response to the query. The ad metadata along with user data and search words define a feature vector that we use to train a logistic regressor that predicts the CTR of future ads. Given the heterogeneity of the components of the feature vector we expect a logistic cost function with skewed level sets and consequent large benefits from the use of oLBFGS.

For the CTR problem considered here we use the Tencent search engine data set Sun (2012). This data set contains the outcomes of 236 million ($236\times 10^{6}$) searches along with information about the ad, the query, and the user. The information contained in each sample point is the following:

User profile: If known, age and gender of visitor performing query.

Depth: Total number of advertisements displayed in the search results page.

Position: Position of the advertisement in the search page.

Impression: Number of times the ad was displayed to the user who issued the query.

Query: The words that appear in the user’s query.

Title: The words that appear in the title of ad.

Keywords: Selected keywords that specify the type of product.

Ad ID: Unique identifier assigned to each specific advertisement.

Advertiser ID: Unique identifier assigned to each specific advertiser.

Clicks: Number of times the user clicked on the ad.

From this information we create a set of feature vectors $\{{\mathbf{x}}_{i}\}_{i=1}^{N}$, with corresponding labels $y_{i}\in\{-1,1\}$. The label associated with feature vector ${\mathbf{x}}_{i}$ is $y_{i}=1$ if the number of clicks in the ad is more than . Otherwise the label is $y_{i}=-1$. We use a binary encoding for all the features in the vector ${\mathbf{x}}_{i}$. For the age of the user we use the six age intervals $(0,12]$, $(12,18]$, $(18,24]$, $(24,30]$, $(30,40]$, and $(40,\infty)$ to construct six indicator entries in ${\mathbf{x}}_{i}$ that take the value 1 if the age of the user is known to be in the corresponding interval. E.g., a 21 year old user has an age that falls in the third interval which implies that we make $[{\mathbf{x}}_{i}]_{3}=1$ and $[{\mathbf{x}}_{i}]_{k}=0$ for all other $k$ between 1 and 6. If the age of the user is unknown we make $[{\mathbf{x}}_{i}]_{k}=0$ for all $k$ between 1 and 6. For the gender of the visitors we use the next three components of ${\mathbf{x}}_{i}$ to indicate male, female, or unknown gender. For a male user we make $[{\mathbf{x}}_{i}]_{7}=1$, for a female user $[{\mathbf{x}}_{i}]_{8}=1$, and for visitors of unknown gender we make $[{\mathbf{x}}_{i}]_{9}=1$. The next three components of ${\mathbf{x}}_{i}$ are used for the depth feature. If the the number of advertisements displayed in the search page is $1$ we make $[{\mathbf{x}}_{i}]_{10}=1$, if $2$ different ads are shown we make $[{\mathbf{x}}_{i}]_{11}=1$, and for depths of $3$ or more we make $[{\mathbf{x}}_{i}]_{12}=1$. To indicate the position of the ad in the search page we also use three components of ${\mathbf{x}}_{i}$. We use $[{\mathbf{x}}_{i}]_{13}=1$, $[{\mathbf{x}}_{i}]_{14}=1$, and $[{\mathbf{x}}_{i}]_{15}=1$ to indicate that the ad is displayed in the first, second, and third position, respectively. Likewise we use $[{\mathbf{x}}_{i}]_{16}$, $[{\mathbf{x}}_{i}]_{17}$ and $[{\mathbf{x}}_{i}]_{18}$ to indicate that the impression of the ad is $1$, $2$ or more than $3$.

For the words that appear in the query we have in the order of $10^{5}$ distinct words. To reduce the number of elements necessary for this encoding we create 20,000 bags of words through random hashing with each bag containing 5 or 6 distinct words. Each of these bags is assigned an index $k$. For each of the words in the query we find the bag in which this word appears. If the word appears in the $k$th bag we indicate this occurrence by setting the $k+18$th component of the feature vector to $[{\mathbf{x}}_{i}]_{k+18}=1$. Observe that since we use 20,000 bags, components 19 through 20,018 of ${\mathbf{x}}_{i}$ indicate the presence of specific words in the query. Further note that we may have more than one ${\mathbf{x}}_{i}$ component different from zero because there may be many words in the query, but that the total number of nonzero elements is much smaller than 20,000. On average, $3.0$ of these elements of the feature vector are nonzero. The same bags of words are used to encode the words that appear in the title of the ad and the product keywords. We encode the words that appear in the title of the ad by using the next $20,000$ components of vector ${\mathbf{x}}_{i}$, i.e. components $20,019$ through $40,018$. Components $40,019$ through $60,018$ are used to encode product keywords. As in the case of the words in the search just a few of these components are nonzero. On average, the number of non-zero components of feature vectors that describe the title features is $8.8$. For product keywords the average is $2.1$. Since the number of distinct advertisers in the training set is $5,184$ we use feature components $60,019$ through $65202$ to encode this information. For the $k$th advertiser ID we set the $k+60,018{th}$ component of the feature vector to $[{\mathbf{x}}_{i}]_{k+60,018}=1$. Since the number of distinct advertisements is $108,824$ we allocate the last $108,824$ components of the feature vector to encode the ad ID. Observe that only one out of $5,184$ advertiser ID components and one of the $108,824$ advertisement ID components are nonzero.

In total, the length of the feature vector is 174,026 where each of the components are either or $1$. The vector is very sparse. We observe a maximum of 148 nonzero elements and an average of $20.9$ nonzero elements in the training set – see Table 1. This is important because the cost of implementing inner products in the oLBFGS training of the logistic regressor that we introduce in the following section is proportional to the number of nonzero elements in ${\mathbf{x}}_{i}$.

2 Logistic regression of click-through rate

We read (39) as stating that for a feature vector ${\mathbf{x}}$ the CTR is determined by the inner product ${\mathbf{x}}^{T}{\mathbf{w}}$ through the given logistic transformation.

Consider now the training set ${\mathcal{S}}=\{({\mathbf{x}}_{i},y_{i})\}_{i=1}^{N}$ which contains $N$ realizations of features ${\mathbf{x}}_{i}$ and respective click outcomes $y_{i}$ and further define the sets ${\mathcal{S}}_{1}:=\{({\mathbf{x}}_{i},y_{i})\in{\mathcal{S}}:y_{i}=1\}$ and ${\mathcal{S}}_{-1}:=\{({\mathbf{x}}_{i},y_{i})\in{\mathcal{S}}:y_{i}=-1\}$ containing clicked and unclicked advertisements, respectively. With the data given in ${\mathcal{S}}$ we define the optimal classifier ${\mathbf{w}}^{*}$ as a maximum likelihood estimate (MLE) of ${\mathbf{w}}$ given the model in (39) and the training set ${\mathcal{S}}$. This MLE can be found as the minimizer of the log-likelihood loss

where we have added the regularization term $\lambda\|{\mathbf{w}}\|^{2}/2$ to disincentivize large values in the weight vector ${\mathbf{w}}^{*}$; see e.g., Ng (2004).

The practical use of (39) and (5.2) is as follows. We use the data collected in the training set ${\mathcal{S}}$ to determine the vector ${\mathbf{w}}^{*}$ in (5.2). When a user issues a query we concatenate the user and query specific elements of the feature vector with the ad specific elements of several candidate ads. We then proceed to display the advertisement with, say, the largest CTR. We can interpret the set ${\mathcal{S}}$ as having been acquired offline or online. In the former case we want to use a stochastic optimization algorithm because computing gradients is infeasible – recall that we are considering training samples with a number of elements $N$ in the order of $10^{6}$. The performance metric of interest in this case is the logistic cost as a function of computational time. If elements of ${\mathcal{S}}$ are acquired online we update ${\mathbf{w}}$ whenever a new vector becomes available so as to adapt to changes in preferences. In this case we want to exploit the information in new samples as much as possible. The correct metric in this case is the logistic cost as a function of the number of feature vectors processed. We use the latter metric for the numerical experiments in the following section.

3 Numerical Results

Out of the $236\times 10^{6}$ in the Tencent dataset we select $10^{6}$ sample points to use as the training set ${\mathcal{S}}$ and $10^{5}$ sample points to use as a test set ${\mathcal{T}}$. To select elements of the training and test set we divide the first $1.1\times 10^{6}$ sample points of the complete dataset in $10^{5}$ consecutive blocks with 11 elements. The first 10 elements of the block are assigned to the training set and the 11th element to the test set. To solve for the optimal classifier we implement SGD and oLBFGS by selecting feature vectors ${\mathbf{x}}_{i}$ at random from the training set ${\mathcal{S}}$. In all of our numerical experiments the regularization parameter in (5.2) is $\lambda=10^{-6}$. The stepsizes for both algorithms are of the form $\epsilon_{t}=\epsilon_{0}T_{0}/(T_{0}+t)$. We set $\epsilon_{0}=10^{-2}$ and $T_{0}=10^{4}$ for oLBFGS and $\epsilon_{0}=10^{-1}$ and $T_{0}=10^{6}$ for SGD. For SGD the sample size in (9) is set to $L=20$ whereas for oLBFGS it is set to $L=100$. The values of parameters $\epsilon_{0}$, $T_{0}$, and $L$ are chosen to yield best convergence times in a rough parameter optimization search. Observe the relatively large values of $L$ that are used to compute stochastic gradients. This is necessary due to the extreme sparsity of the feature vectors ${\mathbf{x}}_{i}$ that contain an average of only $20.9$ nonzero out 174,026 elements. Even when considering $L=100$ vectors they are close to orthogonal. The size of memory for oLBFGS is set to $\tau=10$. With $L=100$ features with an average sparsity of $20.9$ nonzero elements and memory $\tau=10$ the cost of each LBGS iteration is in the order of $2.1\times 10^{4}$ operations.

Figure 5 illustrates the convergence path of SGD and oLBFGS on the advertising training set. We depict the value of the log likelihood objective in (5.2) evaluated at ${\mathbf{w}}={\mathbf{w}}_{t}$ where ${\mathbf{w}}_{t}$ is the classifier iterate determined by SGD or oLBFGS. The horizontal axis is scaled by the number of feature vectors $L$ that are used in the evaluation of stochastic gradients. This results in a plot of log likelihood cost versus the number $Lt$ of feature vectors processed. To read iteration indexes from Figure 5 divide the horizontal axis values by $L=100$ for oLBGS and $L=20$ for SGD. Consistent with the synthetic data results in Section 4, the curvature correction of oLBFGS results in significant reductions in convergence time. For way of illustration observe that after processing $Lt=3\times 10^{4}$ feature vectors the objective value achieved by oLBFGS is $F({\mathbf{w}}_{t})=0.65$, while for SGD it still stands at $F({\mathbf{w}}_{t})=16$ which is a meager reduction from the random initialization point at which $F({\mathbf{w}}_{0})=30$. In fact, oLBFGS converges to the minimum possible log likelihood cost $F({\mathbf{w}}_{t})=0.65$ after processing $1.7\times 10^{4}$ feature vectors. This illustration hints that oLBGS makes better use of the information available in feature vectors.

To corroborate that the advantage of oLBGS is not just an artifact of the structure of the log likelihood cost in (5.2) we process $2\times 10^{4}$ feature vectors with SGD and oLBFGS and evaluate the predictive accuracy of the respective classifiers on the test set. As measures of predictive accuracy we adopt the frequency histogram of the predicted click through rate $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for all clicked ads and the frequency histogram of the complementary predicted click through rate $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for all the ads that were not clicked. To do so we separate the test set by defining the set ${\mathcal{T}}_{1}:=\{({\mathbf{x}}_{i},y_{i})\in{\mathcal{T}}:y_{i}=1\}$ of clicked ads and the set ${\mathcal{T}}_{-1}:=\{({\mathbf{x}}_{i},y_{i})\in{\mathcal{T}}:y_{i}=-1\}$ of ads in the test set that were not clicked. For a given classifier ${\mathbf{w}}$ we compute the predicted probability $\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ for each of the ads in the clicked set ${\mathcal{T}}_{1}$. We then consider a given interval $[a,b]$ and define the frequency histogram of the predicted click through rate as the fraction of clicked ads for which the prediction $\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ falls in $[a,b]$,

where $\#({\mathcal{T}}_{1})$ denotes the cardinality of the set ${\mathcal{T}}_{1}$. Likewise, we consider the ads in the set ${\mathcal{T}}_{-1}$ that were not clicked and compute the prediction $1-\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ on the probability of the ad not being clicked. We then consider a given interval $[a,b]$ and define the frequency histogram ${\mathcal{H}}_{-1}({\mathbf{w}};a,b)$ as the fraction of unclicked ads for which the prediction $1-\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ falls in $[a,b]$,

The histogram ${\mathcal{H}}_{1}({\mathbf{w}};a,b)$ in (41) allows us to study how large the predicted probability $\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ is for the clicked ads. Conversely, the histogram ${\mathcal{H}}_{-1}({\mathbf{w}};a,b)$ in (42) gives an indication of how large the predicted probability $1-\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})$ is for the unclicked ads. An ideal classifier is one for which the frequency counts in ${\mathcal{H}}_{1}({\mathbf{w}};a,b)$ accumulate at $\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})=1$ and for which ${\mathcal{H}}_{-1}({\mathbf{w}};a,b)$ accumulates observations at $1-\text{CTR}({\mathbf{x}}_{i};{\mathbf{w}})=1$. This corresponds to a classifier that predicts a click probability of 1 for all ads that were clicked and a click probability of 0 for all ads that were not clicked.

Fig. 6(a) shows the histograms of predicted click through rate $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for all clicked ads by oLBFGS and SGD classifiers after processing $2\times 10^{4}$ training sample points. oLBFGS classifier for $88\%$ of test points in ${\mathcal{T}}_{1}$ predicts $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ in the interval $[0,0.1]$ and the classifier computed by SGD estimates the click through rate $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ in the same interval for $37\%$ of clicked ads in the test set. These numbers shows the inaccurate click through rate predictions of both classifiers for the test points with label $y=1$. Although, SGD and oLBFGS classifiers have catastrophic performances in predicting click through rate $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for the clicked ads in the test set, they perform well in estimating complementary predicted click through rate $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for the test points with label $y=-1$. This observation implied by Fig. 6(b) which shows the histograms of complementary predicted click through rate $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for all not clicked ads by oLBFGS and SGD classifiers after processing $2\times 10^{4}$ training sample points. As it shows after processing $2\times 10^{4}$ sample points of the training set the predicted probability $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ by the SGD classifier for $38.8\%$ of the test points are in the interval $[0.9,1]$, while for the classifier computed by oLBFGS $97.3\%$ of predicted probability $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ are in the interval $[0.9,1]$ which is a significant performance.

The reason for the inaccurate predictions of both classifiers is that most elements in the training set ${\mathcal{S}}$ are unclicked ads. Thus, the minimizer ${\mathbf{w}}^{*}$ of the log likelihood cost in (5.2) is close to a classifier that predicts $\text{CTR}({\mathbf{x}};{\mathbf{w}}^{*})\approx 0$ for most ads. Indeed, out of the $10^{6}$ elements in the training set, $94.8\%$ of them have labels $y_{i}=-1$ and only the remaining $5.2\times 10^{4}$ feature vectors correspond to clicked ads. To overcome this problem we replicate observations with labels $y_{i}=1$ to balance the representation of both labels in the training set. Equivalently, we introduce a constant $\gamma$ and redefine the log likelihood objective in (5.2) to give a larger weight to feature vectors that correspond to clicked ads,

where we defined $M:=\gamma\#({\mathcal{S}}_{1})+\#({\mathcal{S}}_{-1})$ to account for the replication of clicked featured vectors that is implicit in (43). To implement SGD and oLBFGS in the weighted log function in (43) we need to bias the random choice of feature vector so that vectors in ${\mathcal{S}}_{1}$ are $\gamma$ times more likely to be selected than vectors in ${\mathcal{S}}_{2}$. Although our justification to introduce $\gamma$ is to balance the types of feature vectors, $\gamma$ is just a tradeoff constant to increase the percentage of correct predictions for clicked ads – which is close to zero in Figure 6 – at the cost of reducing the accuracy of correct predictions of unclicked ads – which is close to one in Figure 6.

We repeat the experiment of processing $2\times 10^{4}$ feature vectors that we sumamrized in Figure 6 but now we use the objective cost in (43) instead of the cost in (5.2). We set $\gamma=18.2$ which makes replicated clicked ads as numerous as unclicked ads. The resulting SGD and oLBFGS histograms of the predicted click through rates for all clicked ads and complementary predicted click through rates for all unclicked ads are shown in Figure 7. In particular, Figure 7(a) shows the histograms of predicted click through rate $\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for all clicked ads after processing $2\times 10^{4}$ training sample points. The modification of the log likelihood cost increases the accuracy of the oLBFGS classifier which is now predicting a click probability $\text{CTR}({\mathbf{x}};{\mathbf{w}})\in[0.9,1]$ for $54.7\%$ of the ads that were indeed clicked. There is also improvement for the SGD classifier but the prediction is much less impressive. Only $15.5\%$ of the clicked ads are associated with a click probability prediction in the interval $[0.9,1]$. This improvement is at the cost of reducing the complementary predicted click through rate $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})$ for the ads that were indeed not clicked. However, the classifier computed by oLBFGS after processing $2\times 10^{4}$ feature vectors still predicts a probability $1-\text{CTR}({\mathbf{x}};{\mathbf{w}})\in[0.9,1]$ for $46.3\%$ of the unclicked ads. The corresponding frequency for the SGD classifier is $10.8\%$.

Do note that the relatively high prediction accuracies in Figure 7 are a reflection of sample bias to some extent. Since ads were chosen for display because they were deemed likely to be clicked they are not a completely random test set. Still, the point to be made here is that oLBFGS succeeds in finding an optimal classifier when SGD fails. It would take the processing of about $10^{6}$ feature vectors for SGD to achieve the same accuracy of oLBFGs.

Conclusions

An online limited memory version of the (oL)BFGS algorithm was studied for solving strongly convex optimization problems with stochastic objectives. Almost sure convergence was established by bounding the traces and determinants of curvature estimation matrices under the assumption that sample functions have well behaved Hessians. The convergence rate of oLBFGS was further determined to be at least of order $O(1/t)$ in expectation. This rate is customary of stochastic optimization algorithms which are limited by their ability to smooth out the noise in stochastic gradient estimates. The application of oLBFGS to support vector machines was also developed and numerical tests on synthetic data were provided. The numerical results show that oLBFGS affords important reductions with respect to stochastic gradient descent (SGD) in terms of the number of feature vectors that need to be processed to achieve a target accuracy as well as in the associated execution time. Moreover, oLBFGS also exhibits a significant execution time reduction when compared to other stochastic quasi-Newton methods. These reductions increase with the problem dimension and can become arbitrarily large. A detailed comparison between oLBFGS and SGD for training a logistic regressor in a large scale search engine advertising problem was also presented. The numerical tests show that oLBFGS trains the regressor using less than $1\%$ of the data required by SGD to obtain similar classification accuracy.

We acknowledge the support of the National Science Foundation (NSF CAREER CCF-0952867) and the Office of Naval Research (ONR N00014-12-1-0997).

A Proof of Proposition 1

We begin by observing that the ${\mathbf{p}}_{u}$ sequence in (17) is defined so that we can write ${\mathbf{p}}_{u+1}={\hat{\mathbf{Z}}}_{t-u-1}{\mathbf{p}}_{u}$ with ${\mathbf{p}}_{0}={\mathbf{p}}$. Indeed, use the explicit expression for ${\hat{\mathbf{Z}}}_{t-u-1}$ in (12) to write the product ${\hat{\mathbf{Z}}}_{t-u-1}{\mathbf{p}}_{u}$ as

where the second equality follows from the definition $\alpha_{u}:=\hat{\rho}_{t-u-1}{\mathbf{v}}_{t-u-1}^{T}{\mathbf{p}}_{u}$ and the third equality from the definition of the ${\mathbf{p}}_{u}$ sequence in (17).

Recall now the oLBFGS Hessian inverse approximation expression in (2.2). It follows that for computing the product ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}$ we can multiply each of the $\tau+1$ summands in the right hand side of (2.2) by ${\mathbf{p}}={\mathbf{p}}_{0}$. Implementing this procedure yields

The fundamental observation in (A) is that all summands except the last contain the product ${\hat{\mathbf{Z}}}_{t-1}{\mathbf{p}}_{0}$. This product cannot only be computed efficiently but, as shown in (44), is given by ${\mathbf{p}}_{1}={\hat{\mathbf{Z}}}_{t-1}{\mathbf{p}}_{0}$. A not so fundamental, yet still important observation, is that the last term can be simplified to $\hat{\rho}_{t-1}{\mathbf{v}}_{t-1}{\mathbf{v}}_{t-1}^{T}{\mathbf{p}}_{0}=\alpha_{0}{\mathbf{v}}_{t-1}$ given the definition of $\alpha_{0}:=\hat{\rho}_{t-1}{\mathbf{v}}_{t-1}^{T}{\mathbf{p}}_{0}$. Implementing both of these substitutions in (A) yields

The structure of (A) is analogous to the structure of (A). In all terms except the last two we require determination of the product ${\hat{\mathbf{Z}}}_{t-2}{\mathbf{p}}_{1}$, which, as per (44) can be computed with $2n$ multiplications and is given by ${\mathbf{p}}_{2}={\hat{\mathbf{Z}}}_{t-2}{\mathbf{p}}_{1}$. Likewise, in the second to last term we can simplify the product $\hat{\rho}_{t-2}{\mathbf{v}}_{t-2}{\mathbf{v}}_{t-2}^{T}{\mathbf{p}}_{1}=\alpha_{1}{\mathbf{v}}_{t-2}$ using the definition $\alpha_{1}=\hat{\rho}_{t-2}{\mathbf{v}}_{t-2}^{T}{\mathbf{p}}_{1}$. Implementing these substitutions in (A) yields an expression that is, again, analogous. In all of the resulting summands except the last three we need to compute the product ${\hat{\mathbf{Z}}}_{t-3}{\mathbf{p}}_{2}$, which is given by ${\mathbf{p}}_{3}={\hat{\mathbf{Z}}}_{t-3}{\mathbf{p}}_{2}$ and in the third to last term we can simplify the product $\hat{\rho}_{t-3}{\mathbf{v}}_{t-3}{\mathbf{v}}_{t-3}^{T}{\mathbf{p}}_{2}=\alpha_{2}{\mathbf{v}}_{t-3}$. Repeating this process keeps yielding terms with analogous structure and, after $\tau-1$ repetitions we simplify (A) to

In the first summand in (47) we can substitute the definition of the first element of the ${\mathbf{q}}_{u}$ sequence ${\mathbf{q}}_{0}:={{\hat{\mathbf{B}}}_{t,0}}^{-1}{\mathbf{p}}_{\tau}$. More important, observe that the matrix ${\hat{\mathbf{Z}}}_{t-1}^{T}$ is the first factor in all but the last summand. Likewise, the matrix ${\hat{\mathbf{Z}}}_{t-2}^{T}$ is the second factor in all but the last two summands and, in general, the matrix ${\hat{\mathbf{Z}}}_{t-u}^{T}$ is the $u$th factor in all but the last $u$ summands. Pulling these common factors recursively through (47) it follows that ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}_{t}$ can be equivalently written as

To conclude the proof we just need to note that the recursive definition of ${\mathbf{q}}_{u}$ in (18) is a computation of the nested elements of (48). To see this consider the innermost element of (48) and use the definition of $\beta_{0}:=\hat{\rho}_{t-\tau}{\hat{\mathbf{r}}}_{t-\tau}^{T}{\mathbf{q}}_{0}$ to conclude that $\alpha_{\tau-1}{\mathbf{v}}_{t-\tau}+{\hat{\mathbf{Z}}}_{t-\tau}^{T}{\mathbf{q}}_{0}$ is given by

where in the last equality we use the definition of ${\mathbf{q}}_{1}$ [cf. (18). Substituting this simplification into (48) eliminates the innermost nested term and leads to

Mimicking the computations in (49) we can see that the innermost term in (50) is $\alpha_{\tau-2}{\mathbf{v}}_{t-\tau+1}+{\hat{\mathbf{Z}}}_{t-\tau+1}^{T}{\mathbf{q}}_{1}={\mathbf{q}}_{2}$ and obtain an analogous expression that we can substitute for ${\mathbf{q}}_{3}$ and so on. Repeating this process $\tau-2$ times leads to the last term being ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}=\alpha_{0}{\mathbf{v}}_{t-1}+{\hat{\mathbf{Z}}}_{t-1}^{T}{\mathbf{q}}_{\tau-1}$ which we can write as $\alpha_{0}{\mathbf{v}}_{t-1}+{\hat{\mathbf{Z}}}_{t-1}^{T}{\mathbf{q}}_{\tau-1}={\mathbf{q}}_{\tau}$ by repeating the operations in (49). This final observation yields ${\hat{\mathbf{B}}}_{t}^{-1}{\mathbf{p}}={\mathbf{q}}_{\tau}$.

B Proof of Lemma 2

For given ${\mathbf{w}}_{t}$ and ${\mathbf{w}}_{t+1}$ define the mean instantaneous Hessian ${\hat{\mathbf{G}}}_{t}$ as the average Hessian value along the segment $[{\mathbf{w}}_{t},{\mathbf{w}}_{t+1}]$

Using the definitions of the mean instantaneous Hessian ${\hat{\mathbf{G}}}_{t}$ in (52) as well as the definitions of the stochastic gradient variations ${\hat{\mathbf{r}}}_{t}$ and variable variations ${\mathbf{v}}_{t}$ in (11) and (4) we can rewrite (53) as

The claim in (26) follows from (54) and (55). Indeed, consider the ratio of inner products ${\hat{\mathbf{r}}}_{t}^{T}{\mathbf{v}}_{t}/{\mathbf{v}}_{t}^{T}{\mathbf{v}}_{t}$ and use (54) and the first inequality in (55) to write

It follows that (26) is true for all times $t$.

To prove (27) we operate (54) and (55). Considering the ratio of inner products ${{\hat{\mathbf{r}}}_{t}^{T}{\hat{\mathbf{r}}}_{t}}/{{\hat{\mathbf{r}}}_{t}^{T}{\mathbf{v}}_{t}}$ and observing that (54) states ${\hat{\mathbf{G}}}_{t}{\mathbf{v}}_{t}={\hat{\mathbf{r}}}_{t}$, we can write

Since the mean instantaneous Hessian ${\hat{\mathbf{G}}}_{t}$ is positive definite according to (55), we can define ${\mathbf{z}}_{t}={\hat{\mathbf{G}}}_{t}^{1/2}{\mathbf{v}}_{t}$. Substituting this observation into (57) we can conclude

Observing (58) and the inequalities in (55), it follows that (27) is true.

C Proof of Lemma 3

We begin with the trace upper bound in (29). Consider the recursive update formula for the Hessian approximation ${\hat{\mathbf{B}}}_{t}$ as defined in (28). To simplify notation we define $s$ as a new index such that $s=t-\tau+u$. Introduce this simplified notation in (28) and compute the trace of both sides. Since traces are linear function of their arguments we obtain

Recall that the trace of a matrix product is independent of the order of the factors to conclude that the second summand of (59) can be simplified to

where the second equality follows because ${\mathbf{v}}_{s}^{T}{\hat{\mathbf{B}}}_{t,u}{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}$ is a scalar and the second equality by observing that the term ${\mathbf{v}}_{s}^{T}{\hat{\mathbf{B}}}_{t,u}{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}$ is the inner product of the vector ${\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}$ with itself. Use the same procedure for the last summand of (59) so as to write $\text{tr}({\hat{\mathbf{r}}}_{s}{\hat{\mathbf{r}}}_{s}^{T})={\hat{\mathbf{r}}}_{s}^{T}{\hat{\mathbf{r}}}_{s}=\|{\hat{\mathbf{r}}}_{s}\|^{2}$. Substituting this latter observation as well as (60) into (59) we can simplify the trace of ${\hat{\mathbf{B}}}_{t,u+1}$ to

The second term in the right hand side of (61) is negative because, as we have already shown, the matrix ${\hat{\mathbf{B}}}_{t,u}$ is positive definite. The third term is the one for which we have derived the bound that appears in (27) of Lemma 2. Using this two observations we can conclude that the trace of ${\hat{\mathbf{B}}}_{t,u+1}$ can be bounded as

By considering (62) as a recursive expression for $u=0,\ldots\tau-1$, we can conclude that

To finalize the proof of (29) we need to find a bound for the initial trace $\text{tr}({\hat{\mathbf{B}}}_{t,0})$. To do so we consider the definition ${\hat{\mathbf{B}}}_{t,0}={\mathbf{I}}/\hat{\gamma}_{t}$ with $\hat{\gamma}_{t}$ as given by (19). Using this definition of ${\hat{\mathbf{B}}}_{t,0}$ as a scaled identity it follows that we can write the trace of ${\hat{\mathbf{B}}}_{t,0}$ as

Substituting the definition of $\hat{\gamma}_{t}$ into the rightmost side of (19) it follows that for all times $t\geq 1$,

The term $\|{\hat{\mathbf{r}}}_{t-1}\|^{2}/{\mathbf{v}}_{t-1}^{T}{\hat{\mathbf{r}}}_{t-1}$ in (77) is of the same form of the rightmost term in (61). We can then, as we did in going from (61) to (62) apply the bound that we provide in (27) of Lemma 2 to conclude that for all times $t\geq 1$

Substituting (66) into (63) and pulling common factors leads to the conclusion that for all times $t\geq 1$ and indices $0\leq u\leq\tau$ it holds

We consider now the determinant lower bound in (30). As we did in (59) begin by considering the recursive update in (28) and define $s$ as a new index such that $s=t-\tau+u$ to simplify notation. Compute the determinant of both sides of (28), factorize ${\hat{\mathbf{B}}}_{t,u}$ on the right hand side, and use the fact that the determinant of a product is the product of the determinants to conclude that

To simplify the right hand side of (68) we should first know that for any vectors ${\mathbf{u}}_{1}$, ${\mathbf{u}}_{2}$, ${\mathbf{u}}_{3}$ and ${\mathbf{u}}_{4}$, we can write $\det({\mathbf{I}}+{\mathbf{u}}_{1}{\mathbf{u}}_{2}^{T}+{\mathbf{u}}_{3}{\mathbf{u}}_{4}^{T})=(1+{\mathbf{u}}_{1}^{T}{\mathbf{u}}_{2})(1+{\mathbf{u}}_{3}^{T}{\mathbf{u}}_{4})-({\mathbf{u}}_{1}^{T}{\mathbf{u}}_{4})({\mathbf{u}}_{2}^{T}{\mathbf{u}}_{3})$ – see, e.g., Li and Fukushima (2001), Lemma $3.3$). Setting ${\mathbf{u}}_{1}={\mathbf{v}}_{s}$, ${\mathbf{u}}_{2}={{{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}}/{{\mathbf{v}}_{s}^{T}{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}}}$, ${\mathbf{u}}_{3}={\hat{\mathbf{B}}}_{t,u}^{-1}{\hat{\mathbf{r}}}_{s}$ and ${\mathbf{u}}_{4}={{{\hat{\mathbf{r}}}_{s}}/{{\hat{\mathbf{r}}}_{s}^{T}{\mathbf{v}}_{s}}}$, implies that $\det({\mathbf{I}}+{\mathbf{u}}_{1}{\mathbf{u}}_{2}^{T}+{\mathbf{u}}_{3}{\mathbf{u}}_{4}^{T})$ is equivalent to the last term in the right hand side of (68). Applying these substitutions implies that $(1+{\mathbf{u}}_{1}^{T}{\mathbf{u}}_{2})=1-{\mathbf{v}}_{s}^{T}{{{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}}/{{\mathbf{v}}_{s}{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}}}=0$ and ${\mathbf{u}}_{1}^{T}{\mathbf{u}}_{4}=-{\mathbf{v}}_{s}^{T}{{{\hat{\mathbf{r}}}_{s}}/{{\hat{\mathbf{r}}}_{s}^{T}{\mathbf{v}}_{s}}}=-1$. Hence, the term $\det({\mathbf{I}}+{\mathbf{u}}_{1}{\mathbf{u}}_{2}^{T}+{\mathbf{u}}_{3}{\mathbf{u}}_{4}^{T})$ can be simplified as ${\mathbf{u}}_{2}^{T}{\mathbf{u}}_{3}$. By this simplification we can write the right hand side of (68) as

To further simplify (69) write $({\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s})^{T}={\mathbf{v}}_{s}^{T}{\hat{\mathbf{B}}}_{t,u}^{T}$ and observer that since ${\hat{\mathbf{B}}}_{t,u}$ is symmetric we have ${\hat{\mathbf{B}}}_{t,u}^{T}{\hat{\mathbf{B}}}_{t,u}^{-1}={\hat{\mathbf{B}}}_{t,u}{\hat{\mathbf{B}}}_{t,u}^{-1}={\mathbf{I}}$. Therefore,

Substitute the simplification in (70) for the corresponding factor in (68). Further multiply and divide the right hand side by the nonzero norm $\|{\mathbf{v}}_{s}\|$ and regroup terms to obtain

To bound the third factor in (71) observe that the largest possible value for the normalized quadratic form ${{\mathbf{v}}_{s}^{T}{\hat{\mathbf{B}}}_{t,u}{\mathbf{v}}_{s}}/\|{\mathbf{v}}_{s}\|^{2}$ occurs when ${\mathbf{v}}_{s}$ is an eigenvector of ${\hat{\mathbf{B}}}_{t,u}$ associated with its largest eigenvalue. In such case the value attained is precisely the largest eigenvalue of ${\hat{\mathbf{B}}}_{t,u}$ implying that we can write