Stochastic Quasi-Newton Methods for Nonconvex Stochastic Optimization

Introduction

In this paper, we consider the following stochastic optimization problem:

In this paper, we study stochastic quasi-Newton (SQN) methods for solving the nonconvex stochastic optimization problem (1.1). In the deterministic optimization setting, quasi-Newton methods are more robust and achieve higher accuracy than gradient methods, because they use approximate second-order derivative information. Quasi-Newton methods usually employ the following updates for solving (1.1):

where $B_{k}$ is an approximation to the Hessian matrix $\nabla^{2}f(x_{k})$ at $x_{k}$, or $H_{k}$ is an approximation to $[\nabla^{2}f(x_{k})]^{-1}$. The most widely-used quasi-Newton method, the BFGS method updates $B_{k}$ via

where $s_{k-1}:=x_{k}-x_{k-1}$ and $y_{k-1}:=\nabla f(x_{k})-\nabla f(x_{k-1})$. By using the Sherman-Morrison-Woodbury formula, it is easy to derive that the equivalent update to $H_{k}=B_{k}^{-1}$ is

where $\rho_{k-1}:=1/(s_{k-1}^{\top}y_{k-1})$. For stochastic optimization, there has been some work in designing stochastic quasi-Newton methods that update the iterates via (1.3) using the stochastic gradient $g_{k}$ in place of $\nabla f(x_{k})$. Specific examples include the following. The adaptive subgradient (AdaGrad) method proposed by Duchi, Hazan and Singer , which takes $B_{k}$ to be a diagonal matrix that estimates the diagonal of the square root of the uncentered covariance matrix of the gradients, has been proven to be quite efficient in practice. In , Bordes, Bottou and Gallinari studied SGD with a diagonal rescaling matrix based on the secant condition associated with quasi-Newton methods. Roux and Fitzgibbon discussed the necessity of including both Hessian and covariance matrix information in a stochastic Newton type method. Byrd et al. proposed a quasi-Newton method that uses the sample average approximation (SAA) approach to estimate Hessian-vector multiplications. In , Byrd et al.proposed a stochastic limited-memory BFGS (L-BFGS) method based on SA, and proved its convergence for strongly convex problems. Stochastic BFGS and L-BFGS methods were also studied for online convex optimization by Schraudolph, Yu and Günter in . For strongly convex problems, Mokhtari and Ribeiro proposed a regularized stochastic BFGS method (RES) and analyzed its convergence in and studied an online L-BFGS method in . Recently, Moritz, Nishihara and Jordan proposed a linearly convergent method that integrates the L-BFGS method in with the variance reduction technique (SVRG) proposed by Johnson and Zhang in to alleviate the effect of noisy gradients. A related method that incorporates SVRG into a quasi-Newton method was studied by Lucchi, McWilliams and Hofmann in . In , Gower, Goldfarb and Richtárik proposed a variance reduced block L-BFGS method that converges linearly for convex functions. It should be noted that all of the above stochastic quasi-Newton methods are designed for solving convex or even strongly convex problems.

Challenges. The key challenge in designing stochastic quasi-Newton methods for nonconvex problem lies in the difficulty in preserving the positive-definiteness of $B_{k}$ (and $H_{k}$), due to the non-convexity of the problem and the presence of noise in estimating the gradient. It is known that the BFGS update (1.4) preserves the positive-definiteness of $B_{k}$ as long as the curvature condition

holds, which can be guaranteed for strongly convex problem. For nonconvex problem, the curvature condition (1.6) can be satisfied by performing a line search. However, doing this is no longer feasible for (1.1) in the stochastic setting, because exact function values and gradient information are not available. As a result, an important issue in designing stochastic quasi-Newton methods for nonconvex problems is how to preserve the positive-definiteness of $B_{k}$ (or $H_{k}$) without line search.

Our contributions. Our contributions (and where they appear) in this paper are as follows.

We propose a stochastic damped L-BFGS (SdLBFGS) method that fits into the proposed framework. This method adaptively generates a positive definite matrix $H_{k}$ that approximates the inverse Hessian matrix at the current iterate $x_{k}$. Convergence and complexity results for this method are provided. Moreover, our method does not generate $H_{k}$ explicitly, and only its multiplication with vectors is computed directly. (See Section 3)

Motivated by the recent advance of SVRG for nonconvex minimization , we propose a variance reduced variant of SdLBFGS and analyze its $\mathcal{SFO}$-calls complexity. (See Section 4)

A general framework for stochastic quasi-Newton methods for nonconvex optimization

We now give some assumptions that are required throughout this paper.

where $\sigma>0$ is the noise level of the gradient estimation, and $\xi_{k}$, $k=1,2,\ldots$, are independent samples, and for a given $k$ the random variable $\xi_{k}$ is independent of $\{x_{j}\}_{j=1}^{k}$.

Note that the stochastic BFGS methods studied in require that the noisy gradient is bounded, i.e.,

where $M_{g}>0$ is a constant. Our assumption (2.2) is weaker than (2.3).

Analogous to deterministic quasi-Newton methods, our SQN method takes steps

where $g_{k}$ is defined as a mini-batch estimate of the gradient:

and $\xi_{k,i}$ denotes the random variable generated by the $i$-th sampling in the $k$-th iteration. From AS.2 we can see that $g_{k}$ has the following properties:

There exist two positive constants $C_{l},C_{u}$ such that

We denote by $\xi_{k}=(\xi_{k,1},\ldots,\xi_{k,m_{k}})$, the random samplings in the $k$-th iteration, and denote by $\xi_{[k]}:=(\xi_{1},\ldots,\xi_{k})$, the random samplings in the first $k$ iterations. Since $H_{k}$ is generated iteratively based on historical gradient information by a random process, we make the following assumption on $H_{k}(k\geq 2)$ to control the randomness (note that $H_{1}$ is given in the initialization step).

For any $k\geq 2$, the random variable $H_{k}$ depends only on $\xi_{[k-1]}$.

It then follows directly from AS.4 and (2.6) that

where the expectation is taken with respect to $\xi_{k}$ generated in the computation of $g_{k}$.

We will not specify how to compute $H_{k}$ until Section 3, where a specific updating scheme for $H_{k}$ satisfying both assumptions AS.3 and AS.4 will be proposed.

We now present our SQN method for solving (1.1) as Algorithm 2.1.

In this subsection, we analyze the convergence and complexity of SQN under the condition that the step size $\alpha_{k}$ in (2.4) is diminishing. Specifically, in this subsection we assume $\alpha_{k}$ satisfies the following condition:

which is a standard assumption in stochastic approximation algorithms (see, e.g., ). One very simple choice of $\alpha_{k}$ that satisfies (2.8) is $\alpha_{k}=O(1/k)$.

The following lemma shows that a descent property in terms of the expected objective value holds for SQN. Our analysis is similar to analyses that have been used in .

Suppose that $\{x_{k}\}$ is generated by SQN and assumptions AS.1-4 hold. Further assume that (2.8) holds, and $\alpha_{k}\leq\frac{\underline{\kappa}}{L\bar{\kappa}^{2}}$ for all $k$. (Note that this can be satisfied if $\alpha_{k}$ is non-increasing and the initial step size $\alpha_{1}\leq\frac{\underline{\kappa}}{L\bar{\kappa}^{2}}$). Then the following inequality holds

where the conditional expectation is taken with respect to $\xi_{k}$.

Define $\delta_{k}=g_{k}-\nabla f(x_{k})$. From (2.4), and assumptions AS.1 and AS.3, we have

Taking expectation with respect to $\xi_{k}$ on both sides of (2.10) conditioned on $x_{k}$, we obtain,

which together with (2.11) and AS.3 yields that

Then (2.13) combined with the assumption $\alpha_{k}\leq\frac{\underline{\kappa}}{LC_{u}^{2}}$ implies (2.9). ∎

Before proceeding further, we introduce the definition of a supermartingale (see for more details).

We are now ready to give convergence results for SQN (Algorithm 2.1).

Suppose that assumptions AS.1-4 hold for $\{x_{k}\}$ generated by SQN with batch size $m_{k}={m}$ for all $k$. If the stepsize $\alpha_{k}$ satisfies (2.8) and $\alpha_{k}\leq\frac{\underline{\kappa}}{L\bar{\kappa}^{2}}$ for all $k$, then it holds that

Moreover, there exists a positive constant $M_{f}$ such that

Define $\beta_{k}:=\frac{\alpha_{k}\underline{\kappa}}{2}\|\nabla f(x_{k})\|^{2}$ and $\gamma_{k}:=f(x_{k})+\frac{L\sigma^{2}\bar{\kappa}^{2}}{2m}\sum_{i=k}^{\infty}\alpha_{i}^{2}$. Let ${\cal{F}}_{k}$ be the $\sigma$-algebra measuring $\beta_{k}$, $\gamma_{k}$ and $x_{k}$. From (2.9) we know that for any $k$, it holds that

Since $\sum_{k=1}^{\infty}\alpha_{k}=+\infty$, it follows that (2.14) holds. ∎

Under the assumption (2.3) used in , we now prove a stronger convergence result showing that any limit point of $\{x_{k}\}$ generated by SQN is a stationary point of (1.1) with probability 1.

Assume the same assumptions hold as in Theorem 2.1, and that (2.3) holds. Then

For any given $\epsilon>0$, according to (2.14), there exist infinitely many iterates $x_{k}$ such that $\|\nabla f(x_{k})\|<\epsilon$. Then if (2.18) does not hold, there must exist two infinite sequences of indices $\{m_{i}\}$, $\{n_{i}\}$ with $n_{i}>m_{i}$, such that for $i=1,2,\ldots,$

where the last inequality is due to (2.3) and the convexity of $\|\cdot\|^{2}$. Then it follows from (2.21) that

which together with (2.20) implies that $\|x_{n_{i}}-x_{m_{i}}\|\to 0$ with probability 1, as $i\to+\infty$. Hence, from the Lipschitz continuity of $\nabla f$, it follows that $\|\nabla f(x_{n_{i}})-\nabla f(x_{m_{i}})\|\to 0$ with probability 1 as $i\to+\infty$. However, this contradicts (2.19). Therefore, the assumption that (2.18) does not hold is not true. ∎

Note that our result in Theorem 2.2 is stronger than the ones given in existing works such as and . Moreover, although Bottou also proves that the SA method for nonconvex stochastic optimization with diminishing stepsize is almost surely convergent to stationary point, our analysis requires weaker assumptions. For example, assumes that the objective function is three times continuously differentiable, while our analysis does not require this. Furthermore, we are able to analyze the iteration complexity of SQN, for a specifically chosen step size $\alpha_{k}$ (see Theorem 2.3 below), which is not provided in .

We now analyze the iteration complexity of SQN.

Suppose that assumptions AS.1-4 hold for $\{x_{k}\}$ generated by SQN with batch size $m_{k}={m}$ for all $k$. We also assume that $\alpha_{k}$ is specifically chosen as

with $\beta\in(0.5,1)$. Note that this choice satisfies (2.8) and $\alpha_{k}\leq\frac{\underline{\kappa}}{L\bar{\kappa}^{2}}$ for all $k$. Then

Taking expectation on both sides of (2.9) and summing over $k=1,\ldots,N$ yields

Since $\beta\in(0.5,1)$, it follows that the number of iterations $N$ needed is at most $O(\epsilon^{-\frac{1}{1-\beta}})$. ∎

Note that Theorem 2.3 also provides iteration complexity analysis for the classic SGD method, which can be regarded as a special case of SQN with $H_{k}=I$. To the best of our knowledge, our complexity result in Theorem 2.3 is new for both SGD and stochastic quasi-Newton methods.

2 Complexity of SQN with random output and constant step size

Suppose that assumptions AS.1-4 hold, and that $\alpha_{k}$ in SQN (Algorithm 2.1) is chosen such that $0<\alpha_{k}\leq 2\underline{\kappa}/(L\bar{\kappa}^{2})$ for all $k$ with $\alpha_{k}<2\underline{\kappa}/(L\bar{\kappa}^{2})$ for at least one $k$. Moreover, for a given integer $N$, let $R$ be a random variable with the probability mass function

where $D_{f}:=f(x_{1})-f^{low}$ and the expectation is taken with respect to $R$ and $\xi_{[N]}$. Moreover, if we choose $\alpha_{k}=\underline{\kappa}/(LC_{u}^{2})$ and $m_{k}=m$ for all $k=1,\ldots,N$, then (2.25) reduces to

where $\delta_{k}=g_{k}-\nabla f(x_{k})$. Now summing $k=1,\ldots,N$ and noticing that $\alpha_{k}\leq 2\underline{\kappa}/(L\bar{\kappa}^{2})$, yields

It follows from the definition of $P_{R}$ in (2.24) that

which together with (2.28) implies (2.25). ∎

Note that in Theorem 2.4, $\alpha_{k}$’s are not required to be diminishing, and they can be constant as long as they are upper bounded by $2\underline{\kappa}/(L\bar{\kappa}^{2})$.

We now show that the $\mathcal{SFO}$ complexity of SQN with random output and constant step size is $O(\epsilon^{-2})$.

Assume the conditions in Theorem 2.4 hold, and $\alpha_{k}=\underline{\kappa}/(LC_{u}^{2})$ and $m_{k}=m$ for all $k=1,\ldots,N$. Let $\bar{N}$ be the total number of $\mathcal{SFO}$-calls needed to calculate stochastic gradients $g_{k}$ in SQN (Algorithm 2.1). For a given accuracy tolerance $\epsilon>0$, we assume that

Note that the number of iterations of SQN is at most $N=\lceil\bar{N}/{m}\rceil$. Obviously, $N\geq\bar{N}/(2{m})$. From (2.26) we have that

In Corollary 2.5 we did not consider the $\mathcal{SFO}$-calls that are involved in updating $H_{k}$ in line 3 of SQN. In the next section, we consider a specific updating scheme to generate $H_{k}$, and analyze the total $\mathcal{SFO}$-calls complexity of SQN including the generation of the $H_{k}$.

Stochastic damped L-BFGS method

In this section, we propose a specific way, namely a damped L-BFGS method (SdLBFGS), to generate $H_{k}$ in SQN (Algorithm 2.1) that satisfies assumptions AS.3 and AS.4. We also provide an efficient way to compute $H_{k}g_{k}$ without generating $H_{k}$ explicitly.

Before doing this, we first describe a stochastic damped BFGS method as follows. We generate an auxiliary stochastic gradient at $x_{k}$ using the samplings from the $(k-1)$-st iteration:

Note that we assume that our $\mathcal{SFO}$ can separate two arguments $x_{k}$ and $\xi_{k}$ in the stochastic gradient $g(x_{k},\xi_{k-1})$ and generate an output $g(x_{k};\xi_{k-1,i})$. The stochastic gradient difference is defined as

The iterate difference is still defined as $s_{k-1}=x_{k}-x_{k-1}$. We then define

Note that if $B_{k-1}\succ 0$, then $0<\hat{\theta}_{k-1}\leq 1$. Our stochastic damped BFGS approach updates $B_{k-1}$ as

According to the Sherman-Morrison-Woodbury formula, this corresponds to updating $H_{k}=B_{k}^{-1}$ as

where $\rho_{k-1}=(s_{k-1}^{\top}\bar{y}_{k-1})^{-1}$. The following lemma shows that the damped BFGS updates (3.4) and (3.5) preserve the positive definiteness of $B_{k}$ and $H_{k}$.

For $\bar{y}_{k-1}$ defined in (3.2), $s_{k-1}^{\top}\bar{y}_{k-1}\geq 0.25s_{k-1}^{\top}B_{k-1}s_{k-1}$. Moreover, if $B_{k-1}=H_{k-1}^{-1}\succ 0$, then $B_{k}$ and $H_{k}$ generated by the damped BFGS updates (3.4) and (3.5) are both positive definite.

given that $H_{k-1}\succ 0$. Therefore, both $H_{k}$ and $B_{k}$ defined in (3.5) and (3.4) are positive definite. ∎

where $\rho_{j}=(s_{j}^{\top}y_{j})^{-1}$. The output $H_{k,p}$ is then used as the estimate of the inverse Hessian at $x_{k}$ to compute the search direction at the $k$-th iteration. It can be shown that if the sequence of pairs $\{s_{j},y_{j}\}$ satisfy the curvature condition $s_{j}^{\top}y_{j}>0$, $j=k-1,\ldots,k-p$, then $H_{k,p}$ is positive definite provided that $H_{k,0}$ is positive definite. Recently, stochastic L-BFGS methods have been proposed for solving strongly convex problems in . However, the theoretical convergence analyses in these papers do not apply to nonconvex problems. We now show how to design a stochastic damped L-BFGS formula for nonconvex problems.

Suppose that in the past iterations the algorithm generated $s_{j}$ and $\bar{y}_{j}$ that satisfy

Then at the current iterate, we compute $s_{k-1}=x_{k}-x_{k-1}$ and $y_{k-1}$ by (3.1). Since $s_{k-1}^{\top}y_{k-1}$ may not be positive, motivated by the stochastic damped BFGS update (3.2)-(3.5), we define a new vector $\{\bar{y}_{k-1}\}$ as

Using $s_{j}$ and $\bar{y}_{j}$, $j=k-p,\ldots,k-1$, we define the stochastic damped L-BFGS formula as

where $\rho_{j}=(s_{j}^{\top}\bar{y}_{j})^{-1}$. As in the analysis in Lemma 3.1, by induction we can show that $H_{k,i}\succ 0$, $i=1,\ldots,p$. Note that when $k<p$, we use $s_{j}$ and $\bar{y}_{j}$, $j=1,\ldots,k$ to execute the stochastic damped L-BFGS update.

We next discuss the choice of $H_{k,0}$. A popular choice in the standard L-BFGS method is to set $H_{k,0}=\frac{s_{k-1}^{\top}y_{k-1}}{y_{k-1}^{\top}y_{k-1}}I$. Since $s_{k-1}^{\top}y_{k-1}$ may not be positive for nonconvex problems, we set

To prove that $H_{k}=H_{k,p}$ generated by (3.9)-(3.10) satisfies assumptions AS.3 and AS.4, we need to make the following assumption.

The function $F(x,\xi)$ is twice continuously differentiable with respect to $x$. The stochastic gradient $g(x,\xi)$ is computed as $g(x,\xi)=\nabla_{x}F(x,\xi)$, and there exists a positive constant $\kappa$ such that $\|\nabla_{xx}^{2}F(x,\xi)\|\leq\kappa$, for any $x,\xi$.

Note that AS.5 is equivalent to requiring that $-\kappa I\preceq\nabla_{xx}^{2}F(x,\xi)\preceq\kappa I$, rather than the strong convexity assumption $0\prec\underline{\kappa}I\preceq\nabla_{xx}^{2}F(x,\xi)\preceq\kappa I$ required in . The following lemma shows that the eigenvalues of $H_{k}$ are bounded below away from zero under assumption AS.5.

Suppose that AS.5 holds. Given $H_{k,0}$ defined in (3.10), suppose that $H_{k}=H_{k,p}$ is updated through the stochastic damped L-BFGS formula (3.9). Then all the eigenvalues of $H_{k}$ satisfy

According to Lemma 3.1, $H_{k,i}\succ 0$, $i=1,\ldots,p$. To prove that the eigenvalues of $H_{k}$ are bounded below away from zero, it suffices to prove that the eigenvalues of $B_{k}=H_{k}^{-1}$ are bounded from above. From the damped L-BFGS formula (3.9), $B_{k}=B_{k,p}$ can be computed recursively as

starting from $B_{k,0}=H_{k,0}^{-1}=\gamma_{k}I$. Since $B_{k,0}\succ 0$, Lemma 3.1 indicates that $B_{k,i}\succ 0$ for $i=1,\ldots,p$. Moreover, the following inequalities hold:

From the definition of $\bar{y}_{j}$ in (3.7) and the facts that $s_{j}^{\top}\bar{y}_{j}\geq 0.25s_{j}^{\top}B_{j+1,0}s_{j}$ and $B_{j+1,0}=\gamma_{j+1}I$ from (3.10), we have that for any $j=k-1,\ldots,k-p$

where $\overline{\nabla^{2}_{xx}F}(x_{j},\xi_{j,l},s_{j})=\int_{0}^{1}\nabla_{xx}^{2}F(x_{j}+ts_{j},\xi_{j,l})dt$, because $g(x_{j+1},\xi_{j,l})-g(x_{j},\xi_{j,l})=\int_{0}^{1}\frac{dg}{dt}(x_{j}+ts_{j},\xi_{j,l})dt=\int_{0}^{1}\nabla_{xx}^{2}F(x_{j}+ts_{j},\xi_{j,l})s_{j}dt$. Therefore, for any $j=k-1,\ldots,k-p$, from (3.13), and the facts that $0<\theta_{j}\leq 1$ and $\delta\leq\gamma_{j+1}\leq\kappa+\delta$, and the assumption AS.5 it follows that

We now prove that $H_{k}$ is uniformly bounded above.

Suppose that the assumption AS.5 holds. Given $H_{k,0}$ defined in (3.10), suppose that $H_{k}=H_{k,p}$ is updated through the stochastic damped L-BFGS formula (3.9). Then $H_{k}$ satisfies

where $\alpha=(4\kappa+5\delta)/\delta$, and $\lambda_{\max}(H_{k})$ and $\|H_{k}\|$ denote, respectively, the maximum eigenvalue and operator norm $\|\cdot\|$ of $H_{k}$.

For notational simplicity, let $H=H_{k,i-1}$, $H^{+}=H_{k,i}$, $s=s_{j}$, $\bar{y}=\bar{y}_{j}$, $\rho=(s_{j}^{\top}\bar{y}_{j})^{-1}=(s^{\top}\bar{y})^{-1}$. Now (3.5) can be written as

Using the facts that $\|uv^{\top}\|=\|u\|\cdot\|v\|$ for any vectors $u$ and $v$, $\rho s^{\top}s=\rho\|s\|^{2}=\frac{s^{\top}s}{s^{\top}\bar{y}}\leq\frac{4}{\delta}$, and $\frac{\|\bar{y}\|^{2}}{s^{\top}\bar{y}}\leq 4\left(\frac{\kappa^{2}}{\delta}+\kappa+\delta\right)<\frac{4}{\delta}(\kappa+\delta)^{2}$, which follows from (3.14), we have that

Noting that $\frac{\|\bar{y}\|\|s\|}{s^{\top}\bar{y}}=\left[\frac{\|\bar{y}\|^{2}}{s^{\top}\bar{y}}\cdot\frac{\|s\|^{2}}{s^{\top}\bar{y}}\right]^{1/2}$, it follows that

Lemmas 3.2 and 3.3 indicate that $H_{k}$ generated by (3.7)-(3.9) satisfies assumption AS.3. Moreover, since $y_{k-1}$ defined in (3.1) does not depend on random samplings in the $k$-th iteration, it follows that $H_{k}$ depends only on $\xi_{[k-1]}$ and assumption AS.4 is satisfied.

To analyze the cost of computing the step direction $H_{k}g_{k}$, note that from (3.9), $H_{k}$ can be represented as

which is the same as the classical L-BFGS formula in (3.6), except that $y_{j}$ is replaced by $\bar{y}_{j}$. Hence, we can compute the step direction by the two-loop recursion, implemented in the following procedure.

We now analyze the computational cost of Procedure 3.1. In Step 2, the computation of $\gamma_{k}$ involves $y_{k-1}^{\top}y_{k-1}$ and $s_{k-1}^{\top}y_{k-1}$, which take $2n$ multiplications. In Step 3, from the definition of $\bar{y}_{k}$ in (3.7), since $s_{k-1}^{\top}y_{k-1}$ has been obtained in a previous step, one only needs to compute $s_{k-1}^{\top}s_{k-1}$ and some scalar-vector products, thus the computation of $\bar{y}_{k-1}$ takes $3n$ multiplications. Due to the fact that

all involved computations have been done for $\rho_{k-1}$. Furthermore, the first loop Steps 4-7 involves $p$ scalar-vector multiplications and $p$ vector inner products. So does the second loop Steps 9-12. Including the product $\gamma_{k}^{-1}u_{p}$, the whole procedure takes $(4p+6)n$ multiplications.

Notice that in Step 1 of Procedure 3.1, the computation of $y_{k-1}$ involves the evaluation of $\sum_{i=1}^{m_{k-1}}g(x_{k},\xi_{k-1,i})$, which requires $m_{k-1}$ $\mathcal{SFO}$-calls. As a result, when Procedure 3.1 is plugged into SQN (Algorithm 2.1), the total number of $\mathcal{SFO}$-calls needed in the $k$-th iteration becomes $m_{k}+m_{k-1}$. This leads to the following overall $\mathcal{SFO}$-calls complexity result for our stochastic damped L-BFGS method.

SdLBFGS with a Variance Reduction Technique

Motivated by the recent advance of SVRG for nonconvex minimization proposed in and , we now present a variance reduced SdLBFGS method, which we call SdLBFGS-VR, for solving (1.2). Here, the mini-batch stochastic gradient is defined as $g(x)=\frac{1}{|{\cal{K}}|}\sum_{i\in{\cal{K}}}\nabla f_{i}(x)$, where the subsample set ${\cal{K}}\subseteq[T]$ is randomly chosen from $\{1,\ldots,T\}$. SdLBFGS-VR allows a constant step size, and thus can accelerate the convergence speed of SdLBFGS. SdLBFGS-VR is summarized in Algorithm 4.1.

We now analyze the $\mathcal{SFO}$-calls complexity of Algorithm 4.1. We first analyze the convergence rate of SdLBFGS-VR, essentially following .

Suppose assumptions AS.1, AS.2 and AS.5 hold. Set $c_{t}^{k+1}=c_{t+1}^{k+1}(1+\alpha_{t}^{k+1}\beta_{t}+2L^{2}(\alpha_{t}^{k+1})^{2}\bar{\kappa}^{2}/m)+(\alpha_{t}^{k+1})^{2}L^{3}\bar{\kappa}^{2}/m$. It holds that

It follows from AS.1 and the fact that $\underline{\kappa}I\preceq H_{t}^{k+1}\preceq\bar{\kappa}I$, for any $t=0,\ldots,q-1;k=0,\ldots,N-1$, which follows under assumption AS.5, from Lemmas 3.2 and 3.3, that,

Moreover, for any $\beta_{t}>0$, since $g_{t}^{k+1}$ is an unbiased estimate of $\nabla f(x_{t}^{k+1})$, it holds that,

Furthermore, the following inequality holds:

Combining (4.2), (4.3) and (4.4) yields that,

Suppose assumptions AS.1, AS.2 and AS.5 hold. Set $\beta_{t}=\beta=\frac{L\bar{\kappa}}{T^{1/3}}$, $c_{q}^{k+1}=c_{q}=0$. Suppose that there exist two positive constants $\nu,\mu_{0}\in(0,1)$ such that

holds. Set $\alpha_{t}^{k+1}=\alpha=\frac{\mu_{0}m}{L\bar{\kappa}T^{2/3}}$, and $q=\lfloor\frac{T}{3\mu_{0}m}\rfloor$.

Denote $\theta=\alpha\beta+2\alpha^{2}L^{2}\bar{\kappa}^{2}/m$. It then follows that $\theta=\mu_{0}m/T+2\mu_{0}^{2}m/T^{4/3}\leq 3\mu_{0}m/T$, and $(1+\theta)^{q}\leq e$, where $e$ is the Euler’s number. Because $c_{q}=0$, for any $k\geq 0$, we have

From Theorem 4.1, it follows that to obtain an $\epsilon$-solution, the outer iteration number $N$ of Algorithm 4.1 should be in the order of $O(\frac{T^{2/3}}{qm\epsilon})=O(\frac{T^{-1/3}}{\epsilon}),$ which is due to the fact that $qm=O(T)$. As a result, the total number of component gradient evaluations is $(T+qm)N$, which is $O(T^{2/3}/\epsilon)$. ∎

Numerical Experiments

In this section, we empirically study the performance of the proposed SdLBFGS and SdLBFGS-VR methods. We compare SdLBFGS with SGD with $g_{k}$ given by (2.5) using a diminishing step size $\alpha_{k}=\beta/k$ in both methods, for solving the following nonconvex support vector machine (SVM) problem with a sigmoid loss function, which has been considered in :

In Figure 5.1 we compare the performance of SGD and SdLBFGS with various memory sizes $p$. The batch size was set to $m=100$ and the stepsize to both $\alpha_{k}=10/k$ and $\alpha_{k}=20/k$ for SGD and $10/k$ for SdLBFGS. In the left figure we plot the squared norm of the gradient (SNG) versus the number of iterations, up to a total of $1000$. The SNG was computed using $N=5000$ randomly generated testing points $(u_{i},v_{i}),i=1,\ldots,N$ as:

In Figure 5.2, we report the performance of SdLBFGS with different $\delta$ used in (3.10). From Figure 5.2 we see that SdLBFGS performs best with small $\delta$ such as $\delta=0.01,0.1$ and $1$.

In Figure 5.3 we report the effect of the batch size $m$ on the performance of SGD and SdLBFGS with memory size $p=20$. For SdLBFGS, the left figure shows that $m=500$ gives the best performance among the three choices 50, 100 and 500, tested, with respect to the total number of iterations taken. This is because a larger batch size leads to gradient estimation with lower variance. The right figure shows that if the total number of $\mathcal{SFO}$-calls is fixed, then because of the tradeoff between the number of iterations and the batch size, i.e., because the number of iterations is proportional to the reciprocal of the batch size, the SdLBFGS variant corresponding to $m=100$ slightly outperforms the $m=500$ variant.

In Figure 5.4, we report the percentage of correctly classified data for $5000$ randomly generated testing points. The results are consistent with the one shown in the left figure of Figure 5.3, i.e., the ones with a lower squared norm of the gradient give a higher percentage of correctly classified data.

Moreover, we also counted the number of steps taken by SdLBFGS in which $s_{k-1}^{\top}y_{k-1}<0$. We set the total number of iterations to 1000 and tested the effect of the memory size and batch size of SdLBFGS on the number of such steps. For fixed batch size $m=50$, the average numbers of such steps over 10 runs of SdLBFGS were respectively equal to $(178,50,36,42,15)$ when the memory sizes were $(1,3,5,10,20)$. For fixed memory size $p=20$, the average numbers of such steps over 10 runs of SdLBFGS were respectively equal to $(15,6,1)$ when the batch sizes are $(50,100,500)$. Therefore, the number of such steps roughly decreases as the memory size $p$ and the batch size $m$ increase. This is to be expected because as $p$ increases, there is less negative effect caused by “limited-memory”; and as $m$ increases, the gradient estimation has lower variance.

2 Numerical results for SdLBFGS on the RCV1 dataset

In this subsection, we compare SGD and SdLBFGS for solving (5.1) on a real dataset: RCV1 , which is a collection of newswire articles produced by Reuters in 1996-1997. In our tests, we used a subset downloaded from http://www.cad.zju.edu.cn/home/dengcai/Data/TextData.html of RCV1 used in that contains 9625 articles with 29992 distinct words. The articles are classified into four categories “C15”, “ECAT”, “GCAT” and “MCAT”, each with 2022, 2064, 2901 and 2638 articles respectively. We consider the binary classification problem of predicting whether or not an article is in the second and fourth category, i.e., the entry of each label vector is 1 if a given article appears in category “MCAT” or “ECAT”, and -1 otherwise. We used 60% of the articles (5776) as training data and the remaining 40% (3849) as testing data.

In Figure 5.5, we compare SdLBFGS with various memory sizes and SGD on the RCV1 dataset. For SGD and SdLBFGS, we use the stepsize: $\alpha_{k}=10/k$ and the batch size $m=100$. We also used a second stepsize of $20/k$ for SGD. Note that the SNG computed via (5.3) uses $N=3849$ testing data. The left figure shows that for the RCV1 data set, increasing the memory size improves the performance of SdLBFGS. The performance of SdLBFGS with memory sizes $p=10,20$ was similar, although for $p=10$ it was slightly better. The right figure also shows that larger memory sizes can achieve higher correct classification percentages.

In Figure 5.6, we report the performance of SdLBFGS on RCV1 dataset with different $\delta$ used in (3.10). Similar to Figure 5.2, we see from Figure 5.6 that SdLBFGS works best for small $\delta$ such as $\delta=0.01,0.1$ and $1$.

Figure 5.7 compares SGD and SdLBFGS with different batch sizes. The stepsize of SGD and SdLBFGS was set to $\alpha_{k}=20/k$ and $\alpha_{k}=10/k$, respectively. The memory size of SdLBFGS was chosen as $p=10$. We tested SGD with batch size $m=1,50,100$, and SdLBFGS with batch size $m=1,50,75,100$. From Figure 5.7 we can see that SGD performs worse than SdLBFGS. For SdLBFGS, from Figure 5.7 (a) we observe that larger batch sizes give better results in terms of SNG. If we fix the total number of $\mathcal{SFO}$-calls to $2*10^{4}$, SdLBFGS with $m=1$ performs the worst among the different batch sizes and exhibits dramatic oscillation. The performance gets much better when the batch size becomes larger. In this set of tests, the performance with $m=50,75$ was slightly better than $m=100$. One possible reason is that for the same number of $\mathcal{SFO}$-calls, a smaller batch size leads to larger number of iterations and thus gives better results.

In Figure 5.8, we report the percentage of correctly classified data points for both SGD and SdLBFGS with different batch sizes. These results are consistent with the ones in Figure 5.7. Roughly speaking, the algorithm that gives a lower SNG leads to a higher percentage of correctly classified data points.

We also counted the number of steps taken by SdLBFGS in which $s_{k-1}^{\top}y_{k-1}<0$. We again set the total number of iterations of SdLBFGS to 1000. For fixed batch size $m=50$, the average numbers of such steps over 10 runs of SdLBFGS were respectively equal to $(3,5,8,6,1)$ when the memory sizes were $(1,3,5,10,20)$. For fixed memory size $p=10$, the average numbers of such steps over 10 runs of SdLBFGS were respectively equal to $(308,6,2)$ when the batch sizes were $(1,50,75)$. This is qualitatively similar to our observations in Section 5.1, except that for the fixed batch size $m=50$, a fewer number of such steps were required by the SdLBFGS variants with memory sizes $p=1$ and $p=3$ compared with $p=5$ and $p=10$.

3 Numerical results for SdLBFGS-VR on the RCV1 dataset

Figure 5.9 compares the performance of SdLBFGS-VR with different memory size $p$. It shows that the limited-memory BFGS improves performance, even when $p=1$. Moreover, larger memory size usually provides better performance, but the difference is not very significant.

Figure 5.10 compares the performance of SdLBFGS-VR with different batch sizes $m$ and shows that SdLBFGS-VR is not very sensitive to $m$. In these tests, we always set $q=\lfloor T/m\rfloor$.

The impact of step size on SdLBFGS-VR and SVRG is shown in Figure 5.11 for three step sizes: $\alpha=0.1,0.01$ and $0.001$. Clearly, for the same step size, SdLBFGS-VR gives better result than SVRG. From our numerical tests, we also observed that neither SdLBFGS-VR nor SVRG is stable when $\alpha\geq 1$.

In Figure 5.12, we report the performance of SdLBFGS-VR with different constant step sizes $\alpha$, for $\alpha=0.1,0.01,0.001$ and SdLBFGS with different diminishing step sizes $\beta/k$ for $\beta=10,1,0.1$, since SdLBFGS needs a diminishing step size to guarantee convergence. We see there that SdLBFGS-VR usually performs better than SdLBFGS. The performance of SdLBFGS with $\beta=10$ is in fact already very good, but still inferior to SdLBFGS-VR. This indicates that the variance reduction technique is indeed helpful.

4 Numerical results of SdLBFGS-VR on MNIST dataset

Conclusions

In this paper we proposed a general framework for stochastic quasi-Newton methods for nonconvex stochastic optimization. Global convergence, iteration complexity, and $\mathcal{SFO}$-calls complexity were analyzed under different conditions on the step size and the output of the algorithm. Specifically, a stochastic damped limited memory BFGS method was proposed, which falls under the proposed framework and does not generate $H_{k}$ explicitly. The damping technique was used to preserve the positive definiteness of $H_{k}$, without requiring the original problem to be convex. A variance reduced stochastic L-BFGS method was also proposed for solving the empirical risk minimization problem. Encouraging numerical results were reported for solving nonconvex classification problems using SVM and neural networks.

Acknowledgement

The authors are grateful to two anonymous referees for their insightful comments and constructive suggestions that have improved the presentation of this paper greatly. The authors also thank Conghui Tan for helping conduct the numerical tests in Section 5.4.

References