RES: Regularized Stochastic BFGS Algorithm

I Introduction

Recourse to quasi-Newton methods then arises as a natural alternative. Indeed, quasi-Newton methods achieve superlinear convergence rates in deterministic settings while relying on gradients to compute curvature estimates . Since unbiased gradient estimates are computable at manageable cost, stochastic generalizations of quasi-Newton methods are not difficult to devise . Numerical tests of these methods on simple quadratic objectives suggest that stochastic quasi-Newton methods retain the convergence rate advantages of their deterministic counterparts . The success of these preliminary experiments notwithstanding, stochastic quasi-Newton methods are prone to yield near singular curvature estimates that may result in erratic behavior (see Section V-A).

In this paper we introduce a stochastic regularized version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method to solve problems with the generic structure in (1). The proposed regularization avoids the near-singularity problems of more straightforward extensions and yields an algorithm with provable convergence guarantees when the functions $f({\mathbf{w}},{\boldsymbol{\theta}})$ are strongly convex.

We begin the paper with a brief discussion of SGD (Section II) and deterministic BFGS (Section II-A). The fundamental idea of BFGS is to continuously satisfy a secant condition that captures information on the curvature of the function being minimized while staying close to previous curvature estimates. To regularize deterministic BFGS we retain the secant condition but modify the proximity condition so that eigenvalues of the Hessian approximation matrix stay above a given threshold (Section II-A). This regularized version is leveraged to introduce the regularized stochastic BFGS algorithm (Section II-B). Regularized stochastic BFGS differs from standard BFGS in the use of a regularization to make a bound on the largest eigenvalue of the Hessian inverse approximation matrix and on the use of stochastic gradients in lieu of deterministic gradients for both, the determination of descent directions and the approximation of the objective function’s curvature. We abbreviate regularized stochastic BFGS as RESThe letters “R and “E” appear in “regularized” as well as in the names of Broyden, Fletcher, and Daniel Goldfarb; “S” is for “stochastic” and Shanno..

Convergence properties of RES are then analyzed (Section III). We prove that lower and upper bounds on the Hessians of the sample functions $f({\mathbf{w}},{\boldsymbol{\theta}})$ are sufficient to guarantee convergence to the optimal argument ${\mathbf{w}}^{*}$ with probability 1 over realizations of the sample functions (Theorem 1). We complement this result with a characterization of the convergence rate which is shown to be at least linear in expectation (Theorem 2). Linear expected convergence rates are typical of stochastic optimization algorithms and, in that sense, no better than SGD. Advantages of RES relative to SGD are nevertheless significant, as we establish in numerical results for the minimization of a family of quadratic objective functions of varying dimensionality and condition number (Section IV). As we vary the condition number we observe that for well conditioned objectives RES and SGD exhibit comparable performance, whereas for ill conditioned functions RES outperforms SGD by an order of magnitude (Section IV-A). As we vary problem dimension we observe that SGD becomes unworkable for large dimensional problems. RES however, exhibits manageable degradation as the number of iterations required for convergence doubles when the problem dimension increases by a factor of ten (Section IV-C).

An important example of a class of problems having the form in (1) are support vector machines (SVMs) that reduce binary classification to the determination of a hyperplane that separates points in a given training set; see, e.g., . We adapt RES for SVM problems (Section V) and show the improvement relative to SGD in convergence time, stability, and classification accuracy through numerical analysis (SectionV-A). We also compare RES to standard (non-regularized) stochastic BFGS. The regularization in RES is fundamental in guaranteeing convergence as standard (non-regularized) stochastic BFGS is observed to routinely fail in the computation of a separating hyperplane.

II Algorithm definition

Introducing now a time index $t$, an initial iterate ${\mathbf{w}}_{0}$, and a step size sequence $\epsilon_{t}$, a stochastic gradient descent algorithm is defined by the iteration

A customary step size choice for which (5) holds is to make $\epsilon_{t}=\epsilon_{0}T_{0}/(T_{0}+t)$, for given parameters $\epsilon_{0}$ and $T_{0}$ that control the initial step size and its speed of decrease, respectively. Convergence notwithstanding, the number of iterations required to approximate ${\mathbf{w}}^{*}$ is very large in problems that don’t have small condition numbers. This motivates the alternative methods we discuss in subsequent sections.

To speed up convergence of (4) resort to second order methods is of little use because evaluating Hessians of the objective function is computationally intensive. A better suited methodology is the use of quasi-Newton methods whereby gradient descent directions are premultiplied by a matrix ${\mathbf{B}}_{t}^{-1}$,

The idea is to select positive definite matrices ${\mathbf{B}}_{t}\succ 0$ close to the Hessian of the objective function ${\mathbf{H}}({\mathbf{w}}_{t}):=\nabla^{2}F({\mathbf{w}}_{t})$. Various methods are known to select matrices ${\mathbf{B}}_{t}$, including those by Broyden e.g., ; Davidon, Feletcher, and Powell (DFP) ; and Broyden, Fletcher, Goldfarb, and Shanno (BFGS) e.g., . We work here with the matrices ${\mathbf{B}}_{t}$ used in BFGS since they have been observed to work best in practice .

In BFGS – and all other quasi-Newton methods for that matter – the function’s curvature is approximated by a finite difference. Specifically, define the variable and gradient variations at time $t$ as

respectively, and 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 ${\mathbf{B}}_{t}$ in terms of the Gaussian differential entropy,

The constraint ${\mathbf{Z}}\succeq{\mathbf{0}}$ in (II-A) restricts the feasible space to positive semidefinite matrices whereas the constraint ${\mathbf{Z}}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$ requires ${\mathbf{Z}}$ to satisfy the secant condition. The objective $\text{tr}({\mathbf{B}}_{t}^{-1}{\mathbf{Z}})-\log\det({\mathbf{B}}_{t}^{-1}{\mathbf{Z}})-n$ represents the differential entropy between random variables with zero-mean Gaussian distributions ${\mathcal{N}}({\mathbf{0}},{\mathbf{B}}_{t})$ and ${\mathcal{N}}({\mathbf{0}},{\mathbf{Z}})$ having covariance matrices ${\mathbf{B}}_{t}$ and ${\mathbf{Z}}$. The differential entropy is nonnegative and equal to zero if and only if ${\mathbf{Z}}={\mathbf{B}}_{t}$. The solution ${\mathbf{B}}_{t+1}$ of the semidefinite program in (II-A) is therefore closest to ${\mathbf{B}}_{t}$ in the sense of minimizing the Gaussian differential entropy among all positive semidefinite matrices that satisfy the secant condition ${\mathbf{Z}}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$.

Strongly convex functions are such that the inner product of the gradient and variable variations is positive, i.e., ${\mathbf{v}}_{t}^{T}{\mathbf{r}}_{t}>0$. In that case the matrix ${\mathbf{B}}_{t+1}$ in (II-A) is explicitly given by the update – see, e.g., and the proof of Lemma 1 –,

In principle, the solution to (II-A) could be positive semidefinite but not positive definite, i.e., we can have ${\mathbf{B}}_{t+1}\succeq{\mathbf{0}}$ but ${\mathbf{B}}_{t+1}\not\succ{\mathbf{0}}$. However, through direct operation in (9) it is not difficult to conclude that ${\mathbf{B}}_{t+1}$ stays positive definite if the matrix ${\mathbf{B}}_{t}$ is positive definite. Thus, initializing the curvature estimate with a positive definite matrix ${\mathbf{B}}_{0}\succ{\mathbf{0}}$ guarantees ${\mathbf{B}}_{t}\succ{\mathbf{0}}$ for all subsequent times $t$. Still, it is possible for the smallest eigenvalue of ${\mathbf{B}}_{t}$ to become arbitrarily close to zero which means that the largest eigenvalue of ${\mathbf{B}}_{t}^{-1}$ can become arbitrarily large. This has been proven not to be an issue in BFGS implementations but is a more significant challenge in the stochastic version proposed here.

To avoid this problem we introduce a regularization of (II-A) to enforce the eigenvalues of ${\mathbf{B}}_{t+1}$ to exceed a positive constant $\delta$. Specifically, we redefine ${\mathbf{B}}_{t+1}$ as the solution of the semidefinite program,

The curvature approximation matrix ${\mathbf{B}}_{t+1}$ defined in (II-A) still satisfies the secant condition ${\mathbf{B}}_{t+1}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$ but has a different proximity requirement since instead of comparing ${\mathbf{B}}_{t}$ and ${\mathbf{Z}}$ we compare ${\mathbf{B}}_{t}$ and ${\mathbf{Z}}-\delta{\mathbf{I}}$. While (II-A) does not ensure that all eigenvalues of ${\mathbf{B}}_{t+1}$ exceed $\delta$ we can show that this will be the case under two minimally restrictive assumptions. We do so in the following proposition where we also give an explicit solution for (II-A) analogous to the expression in (9) that solves the non regularized problem in (II-A).

Consider the semidefinite program in (II-A) where the matrix ${\mathbf{B}}_{t}\succ{\mathbf{0}}$ is positive definite and define the corrected gradient variation

Furthermore, ${\mathbf{B}}_{t+1}$ is explicitly given by the expression

II-B RES: Regularized Stochastic BFGS

by using ${\hat{\mathbf{r}}}_{t}$ instead of ${\mathbf{r}}_{t}$. The Hessian approximation ${\hat{\mathbf{B}}}_{t+1}$ for the next iteration is defined as the matrix that satisfies the stochastic secant condition ${\mathbf{Z}}{\mathbf{v}}_{t}={\hat{\mathbf{r}}}_{t}$ and is closest to ${\hat{\mathbf{B}}}_{t}$ in the sense of (II-A). As per Proposition 1 we can compute ${\hat{\mathbf{B}}}_{t+1}$ explicitly as

III Convergence

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

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 (25) as well as the definitions of the stochastic gradient variations ${\hat{\mathbf{r}}}_{t}$ and variable variations ${\mathbf{v}}_{t}$ in (15) and (7) we can rewrite (III) as

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

Consider the RES algorithm as defined by (14)-(17). If assumptions 1, 2 and 3 hold true, the sequence of average function $F({\mathbf{w}}_{t})$ satisfies

where the constant $K:={MS^{2}}({1/\delta}+\Gamma)^{2}/2$.

We proceed to bound the third term in the right hand side of (34). Start by observing that the 2-norm of a product is not larger than the product of the 2-norms and that, as noted above, with ${\mathbf{w}}_{t}$ given the matrix ${\hat{\mathbf{B}}}_{t}^{-1}$ is also given to write

We now find a lower bound for the second term in the right hand side of (III). Since the Hessian approximation matrices ${\hat{\mathbf{B}}}_{t}$ are positive definite their inverses ${\hat{\mathbf{B}}}_{t}^{-1}$ are positive semidefinite. In turn, this implies that all the eigenvalues of ${\hat{\mathbf{B}}}_{t}^{-1}+\Gamma{\mathbf{I}}$ are not smaller than $\Gamma$ since $\Gamma{\mathbf{I}}$ increases all the eigenvalues of ${\hat{\mathbf{B}}}_{t}^{-1}$ by $\Gamma$. This lower bound for the eigenvalues of ${\hat{\mathbf{B}}}_{t}^{-1}+\Gamma{\mathbf{I}}$ implies that

Substituting the lower bound in (38) for the corresponding summand in (III) and further noting the definition of $K:={MS^{2}}({1/\delta}+\Gamma)^{2}/2$ in the statement of the lemma, the result in (33) follows.

Consider the RES algorithm as defined by (14)-(17). If assumptions 1, 2 and 3 hold true and the sequence of stepsizes satisfies (5), the limit infimum of the squared Euclidean distance to optimality $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|^{2}$ satisfies

Proof : The proof uses the relationship in the statement (32) of Lemma 2 to build a supermartingale sequence. For that purpose define the stochastic process $\gamma_{t}$ with values

Observe that $\gamma_{t}$ is well defined because the $\sum_{u=t}^{\infty}{{\epsilon_{u}^{2}}}<\sum_{u=0}^{\infty}{{\epsilon_{u}^{2}}}<\infty$ is summable. Further define the sequence $\beta_{t}$ with values

Let now ${\mathcal{F}}_{t}$ be a sigma-algebra measuring $\gamma_{t}$, $\beta_{t}$, and ${\mathbf{w}}_{t}$. The conditional expectation of $\gamma_{t+1}$ given ${\mathcal{F}}_{t}$ can be written as

because the term $K\sum_{u=t}^{\infty}{{\epsilon_{u}^{2}}}$ is just a deterministic constant. Substituting (32) of Lemma 2 into (42) and using the definitions of $\gamma_{t}$ in (40) and $\beta_{t}$ in (41) yields

Since the sequences $\gamma_{t}$ and $\beta_{t}$ are nonnegative it follows from (43) that they satisfy the conditions of the supermartingale convergence theorem – see e.g. theorem E$7.4$ . Therefore, we conclude that: (i) The sequence $\gamma_{t}$ converges almost surely. (ii) The sum $\sum_{t=0}^{\infty}\beta_{t}<\infty$ is almost surely finite. Using the explicit form of $\beta_{t}$ in (41) we have that $\sum_{t=0}^{\infty}\beta_{t}<\infty$ is equivalent to

Since the sequence of stepsizes is nonsummable for (44) to be true we need to have a vanishing subsequence embedded in $\|\nabla F({\mathbf{w}}_{t})\|^{2}$. By definition, this miles that the limit infimum of the sequence $\|\nabla F({\mathbf{w}}_{t})\|^{2}$ is null,

To transform the gradient bound in (45) into a bound pertaining to the squared distance to optimality $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|^{2}$ simply observe that the lower bound $m$ on the eigenvalues of $\textbf{H}({\mathbf{w}}_{t})$ applied to a Taylor’s expansion around the optimal argument ${\mathbf{w}}^{*}$ implies that

Observe now that since ${\mathbf{w}}^{*}$ is the minimizing argument of $F({\mathbf{w}})$ we must have $F({\mathbf{w}}^{*})-\ F({\mathbf{w}}_{t})\leq 0$ for all ${\mathbf{w}}$. Using this fact and reordering terms we simplify (46) to

Further observe that the Cauchy-Schwarz inequality implies that $\nabla F({\mathbf{w}}_{t})^{T}({\mathbf{w}}_{t}-{\mathbf{w}}^{*})\leq\|\nabla F({\mathbf{w}}_{t})\|\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|$. Substitution of this bound in (47) and simplification of a $\|{\mathbf{w}}^{*}-{\mathbf{w}}_{t}\|$ factor yields

Since the limit infimum of $\|\nabla F({\mathbf{w}}_{t})\|$ is null as stated in (45) the result in (39) follows from considering the bound in (48) in the limit as the iteration index $t\to\infty$.

We complement the convergence result in Theorem 1 with a characterization of the expected convergence rate that we introduce in the following theorem.

Consider the RES algorithm as defined by (14)-(17) and let the sequence of step sizes be given by $\epsilon_{t}=\epsilon_{0}T_{0}/(T_{0}+t)$ with the parameter $\epsilon_{0}$ sufficiently small and the parameter $T_{0}$ sufficiently large so as to satisfy the inequality

Theorem 2 shows that under specified assumptions, the expected error in terms of the objective value after $t$ RES iterations is of order $O(1/t)$. This implies that the rate of convergence for RES is at least linear in expectation. Linear expected convergence rates are typical of stochastic optimization algorithms and, in that sense, no better than conventional SGD. While the convergence rate doesn’t change, improvements in convergence time are marked as we illustrate with the numerical experiments of sections IV and V-A.

IV Numerical analysis

For a given $\rho$ we study the convergence metric

which represents the time needed to achieve a given relative distance to optimality $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|/\|{\mathbf{w}}^{*}\|\leq\rho$ as measured in terms of the number $Lt$ of stochastic functions that are processed to achieve such accuracy.

To study the effect of the problem’s condition number we generate instances of (IV) by choosing ${\mathbf{b}}$ uniformly at random from the box $^{n}$ and the matrix ${\mathbf{A}}$ as diagonal with elements $a_{ii}$ uniformly drawn from the discrete set $\{1,10^{-1},\ldots,10^{-\xi}\}$. This choice of ${\mathbf{A}}$ yields problems with condition number $10^{\xi}$.

As expected for a problem with a large condition number RES is much faster than SGD. After $t=1,200$ the distance to optimality for the SGD iterate is $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|/\|{\mathbf{w}}^{*}\|=3.8\times 10^{-2}$. Comparable accuracy $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|/\|{\mathbf{w}}^{*}\|=3.8\times 10^{-2}$ for RES is achieved after $t=38$ iterations. Since we are using $L=5$ for RES this corresponds to $Lt=190$ random function evaluations. Conversely, upon processing $Lt=1,200$ random functions – which corresponds to $t=240$ iterations – RES achieves accuracy $\|{\mathbf{w}}_{t}-{\mathbf{w}}^{*}\|/\|{\mathbf{w}}^{*}\|=6.6\times 10^{-3}$. This relative performance difference can be made arbitrarily large by modifying the condition number of ${\mathbf{A}}$.

A more comprehensive analysis of the relative advantages of RES appears in figs. 2 and 3. We keep the same parameters used to generate Fig. 1 except that we use $\xi=0$ for Fig. 2 and $\xi=2$ for Fig. 3. This yields a family of well-condition functions with condition number $10^{\xi}=1$ and a family of ill-conditioned functions with condition number $10^{\xi}=10^{2}$. In both figures we consider $\rho=10^{-2}$ and study the convergence times $\tau$ and $\tau^{\prime}$ of RES and SGD, respectively [cf. (53)]. Resulting empirical distributions of $\tau$ and $\tau^{\prime}$ across $J=1,000$ instances of the functions $F({\mathbf{w}})$ in (IV) are reported in figs. 2 and 3 for the well conditioned and ill conditioned families, respectively. For the well conditioned family RES reduces the number of functions processed from an average of $\bar{\tau}^{\prime}=601$ in the case of SGD to an average of $\bar{\tau}=144$. This nondramatic improvement becomes more significant for the ill conditioned family where the reduction is from an average of $\bar{\tau}^{\prime}=7.2\times 10^{3}$ for SGD to an average of $\bar{\tau}=3.2\times 10^{2}$ for RES. The spread in convergence times is also smaller for RES.

IV-B Choice of stochastic gradient average

The trends in convergence times $\tau$ apparent in Fig. 4 are: (i) As we increase $L$ the variance of convergence times decreases. (ii) The average convergence time decreases as we go from small to moderate values of $L$ and starts increasing as we go from moderate to large values of $L$. Indeed, the empirical standard deviations of convergence times decrease monotonically from $\sigma_{\tau_{1}}=2.8\times 10^{3}$ to $\sigma_{\tau_{2}}=2.6\times 10^{2}$, $\sigma_{\tau_{5}}=31.7$, $\sigma_{\tau_{10}}=28.8$, and $\sigma_{\tau_{20}}=22.7$, when $L$ increases from $L=1$ to $L=2$, $L=5$, $L=10$, and $L=20$. The empirical mean decreases from $\bar{\tau}_{1}=3.5\times 10^{3}$ to $\bar{\tau}_{2}=6.3\times 10^{2}$ as we move from $L=1$ to $L=2$, stays at about the same value $\bar{\tau}_{5}=3.3\times 10^{2}$ for $L=5$ and then increases to $\bar{\tau}_{10}=5.8\times 10^{2}$ and $\bar{\tau}_{20}=1.2\times 10^{3}$ for $L=10$ and $L=20$. This behavior is expected since increasing $L$ results in curvature estimates ${\hat{\mathbf{B}}}_{t}$ closer to the Hessian ${\mathbf{H}}({\mathbf{w}}_{t})$ thereby yielding better convergence times. As we keep increasing $L$, there is no payoff in terms of better curvature estimates and we just pay a penalty in terms of more function evaluations for an equally good ${\hat{\mathbf{B}}}_{t}$ matrix. This can be corroborated by observing that the convergence times $\tau_{5}$ are about half those of $\tau_{10}$ which in turn are about half those of $\tau_{20}$. This means that the actual convergence times $\tau/L$ have similar distributions for $L=5$, $L=10$, and $L=20$. The empirical distributions in Fig. 4 show that moderate values of $L$ suffice to provide workable curvature approximations. This justifies the use $L=5$ in sections IV-A and IV-C

IV-C Effect of problem’s dimension

To evaluate performance for problems of different dimensions we consider functions of the form in (IV) with ${\mathbf{b}}$ uniformly chosen from the box $^{n}$ and diagonal matrix ${\mathbf{A}}$ as in Section IV-A. However, we select the elements $a_{ii}$ as uniformly drawn from the interval $$. This results in problems with more moderate condition numbers and allows for a comparative study of performance degradations of RES and SGD as the problem dimension$n$ grows.

The variability parameter for the random vector $\boldsymbol{\theta}$ is set to $\theta_{0}=0.5$. The RES parameters are $L=5$, $\delta=10^{-3}$, and $\Gamma=10^{-4}$. For SGD we use $L=1$. In both methods the step size sequence is $\epsilon_{t}=\epsilon_{0}T_{0}/(T_{0}+t)$ with $\epsilon_{0}=10^{-1}$ and $T_{0}=10^{3}$. For a problem of dimension $n$ we study convergence times $\tau_{n}$ and $\tau^{\prime}_{n}$ of RES and SGD as defined in (53) with $\rho=1$. For each value of $n$ considered we determine empirical distributions of $\tau_{n}$ and $\tau^{\prime}_{n}$ across $J=1,000$ problem instances. If $\tau>5\times 10^{5}$ we report $\tau=5\times 10^{5}$ and interpret this outcome as a convergence failure. The resulting histograms are shown in Fig. 5 for $n=5$, $n=10$, $n=20$, and $n=50$.

For problems of small dimension having $n=5$ the average performances of RES and SGD are comparable, with SGD performing slightly better. E.g., the medians of these times are $\text{median}(\tau_{5})=400$ and $\text{median}(\tau^{\prime}_{5})=265$, respectively. A more significant difference is that times $\tau_{5}$ of RES are more concentrated than times $\tau^{\prime}_{5}$ of SGD. The latter exhibits large convergence times $\tau^{\prime}_{5}>10^{3}$ with probability $0.06$ and fails to converge altogether in a few rare instances – we have $\tau^{\prime}_{5}=5\times 10^{5}$ in 1 out of 1,000 realizations. In the case of RES all realizations of $\tau_{5}$ are in the interval $70\leq\tau_{5}\leq 1095$.

As we increase $n$ we see that RES retains the smaller spread advantage while eventually exhibiting better average performance as well. Medians for $n=10$ are still comparable at $\text{median}(\tau_{10})=575$ and $\text{median}(\tau^{\prime}_{10})=582$, as well as for $n=20$ at $\text{median}(\tau_{20})=745$ and $\text{median}(\tau^{\prime}_{20})=1427$. For $n=50$ the RES median is decidedly better since $\text{median}(\tau_{50})=950$ and $\text{median}(\tau^{\prime}_{50})=7942$.

For large dimensional problems having $n=50$ SGD becomes unworkable. It fails to achieve convergence in $5\times 10^{5}$ iterations with probability $0.07$ and exceeds $10^{4}$ iterations with probability $0.45$. For RES we fail to achieve convergence in $5\times 10^{5}$ iterations with probability $3\times 10^{-3}$ and achieve convergence in less than $10^{4}$ iterations in all other cases. Further observe that RES degrades smoothly as $n$ increases. The median number of gradient evaluations needed to achieve convergence increases by a factor of $\text{median}(\tau^{\prime}_{50})/\text{median}(\tau^{\prime}_{5})=29.9$ as we increase $n$ by a factor of $10$. The spread in convergence times remains stable as $n$ grows.

V Support vector machines

where we also added the regularization term ${\lambda}\|{\mathbf{w}}\|^{2}/{2}$ for some constant $\lambda>0$. The vector ${\mathbf{w}}^{*}$ in (54) balances the minimization of the sum of distances to the separating hyperplane, as measured by the loss function $l(({\mathbf{x}},y);{\mathbf{w}})$, with the minimization of the $L_{2}$ norm $\|{\mathbf{w}}\|_{2}$ to enforce desirable properties in ${\mathbf{w}}^{*}$. 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}}))$, the squared hinge loss $l(({\mathbf{x}},y);{\mathbf{w}})=\max(0,1-y({\mathbf{w}}^{T}{\mathbf{x}}))^{2}$ and the log loss $l(({\mathbf{x}},y);{\mathbf{w}})=\log(1+\exp(-y({\mathbf{w}}^{T}{\mathbf{x}})))$. See, e.g., .

In order to model (54) as a stochastic optimization problem in the form of problem (1), we define $\boldsymbol{\theta}_{i}=({\mathbf{x}}_{i},y_{i})$ as a given training point and $m_{\boldsymbol{\theta}}(\boldsymbol{\theta})$ as a uniform probability distribution on the training set ${\mathcal{S}}=\{({\mathbf{x}}_{i},y_{i})\}_{i=1}^{N}=\{\boldsymbol{\theta}_{i}\}_{i=1}^{N}$. Upon defining the sample functions

it follows that we can rewrite the objective function in (54) as

since each of the functions $f({\mathbf{w}},\boldsymbol{\theta})$ is drawn with probability $1/N$ according to the definition of $m_{\boldsymbol{\theta}}(\boldsymbol{\theta})$. Substituting (56) into (54) yields a problem with the general form of (1) with random functions $f({\mathbf{w}},\boldsymbol{\theta})$ explicitly given by (55).

The specific form of Step 5 is obtained by replacing ${\mathbf{w}}_{t+1}$ for ${\mathbf{w}}_{t}$ in (V). We analyze the behavior of Algorithm (1) in the implementation of a SVM in the following section.

An illustration of the relative performances of SGD and RES for $n\!=4$ is presented in Fig. 6. The value of the objective function $F({\mathbf{w}}_{t})$ is represented with respect to the number of feature vectors processed, which is given by the product $Lt$ between the iteration index and the sample size used to compute stochastic gradients. This is done because the sample sizes in RES ($L=5$) and SGD ($L=1$) are different. The curvature correction of RES results in significant reductions in convergence time. E.g., RES achieves an objective value of $F({\mathbf{w}}_{t})=6.5\times 10^{-2}$ upon processing of $Lt=315$ feature vectors. To achieve the same objective value $F({\mathbf{w}}_{t})=6.5\times 10^{-2}$ SGD processes $1.74\times 10^{3}$ feature vectors. Conversely, after processing $Lt=2.5\times 10^{3}$ feature vectors the objective values achieved by RES and SGD are $F({\mathbf{w}}_{t})=4.14\times 10^{-2}$ and $F({\mathbf{w}}_{t})=6.31\times 10^{-2}$, respectively.

The performance difference between the two methods is larger for feature vectors of larger dimension $n$. The plot of the value of the objective function $F({\mathbf{w}}_{t})$ with respect to the number of feature vectors processed $Lt$ is shown in Fig. 7 for $n=40$. The convergence time of RES increases but is still acceptable. For SGD the algorithm becomes unworkable. After processing $3.5\times 10^{3}$ RES reduces the objective value to $F({\mathbf{w}}_{t})=5.55\times 10^{-4}$ while SGD has barely made progress at $F({\mathbf{w}}_{t})=1.80\times 10^{-2}$.

Differences in convergence times translate into differences in classification accuracy when we process all $N$ vectors in the training set. This is shown for dimension $n=4$ and training set size $N=2.5\times 10^{3}$ in Fig. 8. To build Fig. 8 we process $N=2.5\times 10^{3}$ feature vectors with RES and SGD with the same parameters used in Fig. 6. We then use these vectors to classify $10^{4}$ observations in the test set and record the percentage of samples that are correctly classified. The process is repeated $10^{3}$ times to estimate the probability distribution of the correct classification percentage represented by the histograms shown. The dominance of RES with respect to SGD is almost uniform. The vector ${\mathbf{w}}_{t}$ computed by SGD classifies correctly at most $65\%$ of the of the feature vectors in the test set. The vector ${\mathbf{w}}_{t}$ computed by RES exceeds this accuracy with probability $0.98$. Perhaps more relevant, the classifier computed by RES achieves a mean classification accuracy of $82.2\%$ which is not far from the clairvoyant classification accuracy of $98\%$. Although performance is markedly better in general, RES fails to compute a working classifier with probability $0.02$. We omit comparison of classification accuracy for $n=40$ due to space considerations. As suggested by Fig. 7 the differences are more significant than for the case $n=4$.

We also investigate the difference between regularized and non-regularized versions of stochastic BFGS for feature vectors of dimension $n=10$. Observe that non-regularized stochastic BFGS corresponds to making $\delta=0$ and $\Gamma=0$ in Algorithm 1. To illustrate the advantage of the regularization induced by the proximity requirement in (II-A), as opposed to the non regularized proximity requirement in (II-A), we keep a constant stepsize $\epsilon_{t}=10^{-1}$. The corresponding evolutions of the objective function values $F({\mathbf{w}}_{t})$ with respect to the number of feature vectors processed $Lt$ are shown in Fig. 9 along with the values associated with stochastic gradient descent. As we reach convergence the likelihood of having small eigenvalues appearing in ${\hat{\mathbf{B}}}_{t}$ becomes significant. In regularized stochastic BFGS (RES) this results in recurrent jumps away from the optimal classifier ${\mathbf{w}}^{*}$. However, the regularization term limits the size of the jumps and further permits the algorithm to consistently recover a reasonable curvature estimate. In Fig. 9 we process $10^{4}$ feature vectors and observe many occurrences of small eigenvalues. However, the algorithm always recovers and heads back to a good approximation of ${\mathbf{w}}^{*}$. In the absence of regularization small eigenvalues in ${\hat{\mathbf{B}}}_{t}$ result in larger jumps away from ${\mathbf{w}}^{*}$. This not only sets back the algorithm by a much larger amount than in the regularized case but also results in a catastrophic deterioration of the curvature approximation matrix ${\hat{\mathbf{B}}}_{t}$. In Fig. 9 we observe recovery after the first two occurrences of small eigenvalues but eventually there is a catastrophic deviation after which non-regularized stochastic BFSG behaves not better than SGD.

VI Conclusions

Convex optimization problems with stochastic objectives were considered. RES, a stochastic implementation of a regularized version of the Broyden-Fletcher-Goldfarb-Shanno quasi-Newton method was introduced to find corresponding optimal arguments. Almost sure convergence was established under the assumption that sample functions have well behaved Hessians. A linear convergence rate in expectation was further proven. Numerical results showed that RES affords important reductions in terms of convergence time relative to stochastic gradient descent. These reductions are of particular significance for problems with large condition numbers or large dimensionality since RES exhibits remarkable stability in terms of the total number of iterations required to achieve target accuracies. An application of RES to support vector machines was also developed. In this particular case the advantages of RES manifest in improvements of classification accuracies for training sets of fixed cardinality. Future research directions include the development of limited memory versions as well as distributed versions where the function to be minimized is spread over agents of a network.

Appendix A: Proof of Proposition 1

We first show that (13) is true. Since the optimization problem in (II-A) is convex in ${\mathbf{Z}}$ we can determine the optimal variable ${\mathbf{B}}_{t+1}={\mathbf{Z}}^{*}$ using Lagrangian duality. Introduce then the multiplier variable $\boldsymbol{\mu}$ associated with the secant constraint ${\mathbf{Z}}{\mathbf{v}}_{t}={\mathbf{r}}_{t}$ in (II-A) and define the Lagrangian

The dual function is defined as $d(\boldsymbol{\mu}):=\min_{{\mathbf{Z}}\succeq{\mathbf{0}}}{\mathcal{L}}({\mathbf{Z}},\boldsymbol{\mu})$ and the optimal dual variable is $\boldsymbol{\mu}^{*}:=\operatornamewithlimits{argmin}_{\boldsymbol{\mu}}d(\boldsymbol{\mu})$. We further define the primal Lagrangian minimizer associated with dual variable $\boldsymbol{\mu}$ as

Observe that combining the definitions in (59) and (Appendix A: Proof of Proposition 1) we can write the dual function $d(\boldsymbol{\mu})$ as

We will determine the optimal Hessian approximation ${\mathbf{Z}}^{*}={\mathbf{Z}}(\boldsymbol{\mu}^{*})$ as the Lagrangian minimizer associated with the optimal dual variable $\boldsymbol{\mu}^{*}$. To do so we first find the Lagrangian minimizer (59) by nulling the gradient of ${\mathcal{L}}({\mathbf{Z}},\boldsymbol{\mu})$ with respect to ${\mathbf{Z}}$ in order to show that ${\mathbf{Z}}(\boldsymbol{\mu})$ must satisfy

Multiplying the equality in (61) by ${\mathbf{B}}_{t}$ from the right and rearranging terms it follows that the inverse of the argument of the log-determinant function in (Appendix A: Proof of Proposition 1) can be written as

If, instead, we multiply (61) by $({\mathbf{Z}}(\boldsymbol{\mu})-\delta{\mathbf{I}})$ from the right it follows after rearranging terms that

Further considering the trace of both sides of (63) and noting that $\text{tr}({\mathbf{I}})=n$ we can write the trace in (Appendix A: Proof of Proposition 1) as

Observe now that since the trace of a product is invariant under cyclic permutations of its arguments and the matrix ${\mathbf{Z}}$ is symmetric we have $\text{tr}[\boldsymbol{\mu}{\mathbf{v}}_{t}^{T}({\mathbf{Z}}(\boldsymbol{\mu})-\delta{\mathbf{I}})]=\text{tr}[{\mathbf{v}}\boldsymbol{\mu}_{t}^{T}({\mathbf{Z}}(\boldsymbol{\mu})-\delta{\mathbf{I}})]=\text{tr}[\boldsymbol{\mu}^{T}({\mathbf{Z}}(\boldsymbol{\mu})-\delta{\mathbf{I}}){\mathbf{v}}_{t}]$. Since the argument in the latter is a scalar the trace operation is inconsequential from where it follows that we can rewrite (64) as

Observing that the log-determinant of a matrix is the opposite of the log-determinant of its inverse we can substitute (62) for the argument of the log-determinant in (Appendix A: Proof of Proposition 1). Further substituting (65) for the trace in (Appendix A: Proof of Proposition 1) and rearranging terms yields the explicit expression for the dual function

In order to compute the optimal dual variable $\boldsymbol{\mu}^{*}$ we set the gradient of (66) to zero and manipulate terms to obtain

Applying the Sherman-Morrison formula to compute the inverse of the right hand side of (68) leads to

which can be verified by direct multiplication. The result in (13) follows after solving (69) for ${\mathbf{Z}}(\boldsymbol{\mu}^{*})$ and noting that for the convex optimization problem in (II-A) we must have ${\mathbf{Z}}(\boldsymbol{\mu}^{*})={\mathbf{Z}}^{*}={\mathbf{B}}_{t+1}$ as we already argued.

Consider now the term ${\mathbf{B}}_{t}-{{{\mathbf{B}}_{t}{\mathbf{v}}_{t}{\mathbf{v}}_{t}^{T}{{\mathbf{B}}_{t}}}/{{\mathbf{v}}_{t}^{T}{\mathbf{B}}_{t}{\mathbf{v}}_{t}}}$ and factorize ${\mathbf{B}}_{t}^{{1}/{2}}$ from the left and right side so as to write

Define the vector ${\mathbf{u}}_{t}:={\mathbf{B}}_{t}^{1/2}{\mathbf{v}}_{t}$ and write ${\mathbf{v}}_{t}^{T}{\mathbf{B}}_{t}{\mathbf{v}}_{t}=({\mathbf{B}}_{t}^{1/2}{\mathbf{v}}_{t})^{T}({\mathbf{B}}_{t}^{1/2}{\mathbf{v}}_{t})={\mathbf{u}}_{t}^{T}{\mathbf{u}}_{t}$ as well as ${\mathbf{B}}_{t}^{1/2}{\mathbf{v}}_{t}{\mathbf{v}}_{t}^{T}{\mathbf{B}}_{t}^{1/2}={\mathbf{u}}_{t}{\mathbf{u}}_{t}^{T}$. Substituting these observation into (71) we can conclude that

because the eigenvalues of the matrix ${\mathbf{u}}_{t}{\mathbf{u}}_{t}^{T}/{\mathbf{u}}_{t}^{T}{\mathbf{u}}_{t}$ belong to the interval $$. The only term in (13) which has not been considered is$\delta{\mathbf{I}}$. Since the rest add up to a positive semidefinite matrix it then must be that (12) is true.

Appendix B: Proof of Theorem 2

Let $c>1$, $b>0$ and $t_{0}>0$ be given constants and $u_{t}\geq 0$ be a nonnegative sequence that satisfies the inequality

for all times $t\geq 0$. The sequence $u_{t}$ is then bounded as

for all times $t\geq 0$, where the constant $Q$ is defined as

Proof : We prove (74) using induction. To prove the claim for $t=0$ simply observe that the definition of $Q$ in (75) implies that

because the maximum of two numbers is at least equal to both of them. By rearranging the terms in (76) we can conclude that

Comparing (77) and (74) it follows that the latter inequality is true for $t=0$.

Introduce now the induction hypothesis that (74) is true for $t=s$. To show that this implies that (74) is also true for $t=s+1$ substitute the induction hypothesis $u_{s}\leq Q/(s+t_{0})$ into the recursive relationship in (73). This substitution shows that $u_{s+1}$ is bounded as

Observe now that according to the definition of $Q$ in (75), we know that $b/(c-1)\leq Q$ because $Q$ is the maximum of $b/(c-1)$ and $t_{0}u_{0}$. Reorder this bound to show that $b\leq Q(c-1)$ and substitute into (78) to write

Pulling out $Q/(s+t_{0})^{2}$ as a common factor and simplifying and reordering terms it follows that (79) is equivalent to

To complete the induction step use the difference of squares formula for $(s+t_{0})^{2}-1$ to conclude that

Reordering terms in (81) it follows that \big{[}(s+t_{0})-1\big{]}/(s+t_{0})^{2}\leq 1/\big{[}(s+t_{0})+1\big{]}, which upon substitution into (80) leads to the conclusion that

Eq. (82) implies that the assumed validity of (74) for $t=s$ implies the validity of (74) for $t=s+1$. Combined with the validity of (74) for $t=0$, which was already proved, it follows that (74) is true for all times $t\geq 0$. ∎

Proof of Theorem 2: Consider the result in (32) of Lemma 2 and subtract the average function optimal value $F({\mathbf{w}}^{*})$ from both sides of the inequality to conclude that the sequence of optimality gaps in the RES algorithm satisfies

where, we recall, $K:={MS^{2}}(({1/\delta})+\Gamma)^{2}/2$ by definition.

We proceed to find a lower bound for the gradient norm $\|\nabla F({\mathbf{w}}_{t})\|$ in terms of the error of the objective value $F({\mathbf{w}}_{t})-\ F({\mathbf{w}}^{*})$ – this is a standard derivation which we include for completeness, see, e.g., . As it follows from Assumption 1 the eigenvalues of the Hessian ${\mathbf{H}}({\mathbf{w}}_{t})$ are bounded between $0<m$ and $M<\infty$ as stated in (22). Taking a Taylor’s expansion of the objective function $F({\mathbf{y}})$ around ${\mathbf{w}}$ and using the lower bound in the Hessian eigenvalues we can write

For fixed ${\mathbf{w}}$, the right hand side of (84) is a quadratic function of ${\mathbf{y}}$ whose minimum argument we can find by setting its gradient to zero. Doing this yields the minimizing argument ${\hat{\mathbf{y}}}={\mathbf{w}}-(1/m)\nabla F({\mathbf{w}})$ implying that for all ${\mathbf{y}}$ we must have

The bound in (Appendix B: Proof of Theorem 2) is true for all ${\mathbf{w}}$ and ${\mathbf{y}}$. In particular, for ${\mathbf{y}}={\mathbf{w}}^{*}$ and ${\mathbf{w}}={\mathbf{w}}_{t}$ (Appendix B: Proof of Theorem 2) yields

Rearrange terms in (86) to obtain a bound on the gradient norm squared $\|\nabla F({\mathbf{w}}_{t})\|^{2}$. Further substitute the result in (83) and regroup terms to obtain the bound

Furhter substituting $\epsilon_{t}\!=\!\epsilon_{0}T_{0}/(T_{0}+t)$, which is the assumed form of the step size sequence by hypothesis, we can rewrite (88) as

References