Exploiting Negative Curvature in Deterministic and Stochastic Optimization

TL;DR

This paper introduces new algorithms that leverage negative curvature directions in optimization, demonstrating improved performance in both deterministic and stochastic settings through decision-making based on objective function models.

Contribution

The paper proposes novel frameworks for combining descent and negative curvature steps, enabling fixed stepsizes and applicability to stochastic optimization.

Findings

Deterministic instances show performance gains over descent-only methods.

In stochastic settings, the methods improve progress when standard stochastic-gradient methods are slow.

Frameworks are based on upper-bounding models, allowing theoretical and practical advantages.

Abstract

This paper addresses the question of whether it can be beneficial for an optimization algorithm to follow directions of negative curvature. Although prior work has established convergence results for algorithms that integrate both descent and negative curvature steps, there has not yet been extensive numerical evidence showing that such methods offer consistent performance improvements. In this paper, we present new frameworks for combining descent and negative curvature directions: alternating two-step approaches and dynamic step approaches. The aspect that distinguishes our approaches from ones previously proposed is that they make algorithmic decisions based on (estimated) upper-bounding models of the objective function. A consequence of this aspect is that our frameworks can, in theory, employ fixed stepsizes, which makes the methods readily translatable from deterministic to stochastic settings. For deterministic problems, we show that instances of our dynamic framework yield gains in performance compared to related methods that only follow descent steps. We also show that gains can be made in a stochastic setting in cases when a standard stochastic-gradient-type method might make slow progress.

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }f(x),$$

$$d_k^{T}H(x_k)d_k$$

$$g(x_k{)}^T{d}_k$$

$$\text{and}\Vert d_k{\Vert }_2$$

$$\frac{-g({\hat{x}}_k{)}^T{\hat{s}}_k}{\Vert {\hat{s}}_k{\Vert }_2\Vert g({\hat{x}}_k){\Vert }_2}\ge \delta \text{and}\zeta \le \frac{\Vert {\hat{s}}_k{\Vert }_2}{\Vert g({\hat{x}}_k){\Vert }_2}\le \eta$$

$$\mathcal{D}:=\{k\in {\mathbb{N}}_{+}:d_k\ne 0\}\equiv \{k\in {\mathbb{N}}_{+}:{\lambda }_k<0\}$$

$$\underset{k\to \infty }{\lim }\Vert g(x_k){\Vert }_2=0\text{and}\underset{k\to \infty }{\mathrm{liminf}}{\lambda }_k\ge 0.$$

$$f({\hat{x}}_k)$$

$$\le f(x_k)+\frac{1}{2}{\beta }^2\gamma {\lambda }_k\Vert d_k{\Vert }_2^2+\frac{1}{6}\sigma {\beta }^3{\theta }^3|{\lambda }_k{|}^3$$

$$\displaystyle\leq f(x_{k})-\tfrac{1}{2}\beta^{2}\theta^{2}\big{(}\gamma-\tfrac{1}{3}\sigma\beta\theta\big{)}|\lambda_{k}|^{3}$$

$$=f(x_k)-c_1(\beta )|{\lambda }_k{|}^3,$$

$$f(x_{k+1})$$

$$\le f({\hat{x}}_k)-\alpha \delta \zeta \Vert g({\hat{x}}_k){\Vert }_2^2+\frac{1}{2}L{\alpha }^2{\eta }^2\Vert g({\hat{x}}_k){\Vert }_2^2$$

$$=f({\hat{x}}_k)-\alpha (\delta \zeta -\frac{1}{2}L\alpha {\eta }^2)\Vert g({\hat{x}}_k){\Vert }_2^2$$

$$=f({\hat{x}}_k)-c_2(\alpha )\Vert g({\hat{x}}_k){\Vert }_2^2,$$

$$f(x_{k+1})\le f(x_k)-{\mathcal{I}}_{\mathcal{D}}(k)c_1(\beta )|{\lambda }_k{|}^3-c_2(\alpha )\Vert g({\hat{x}}_k){\Vert }_2^2.$$

$$f(x_1)-f(x_{\ell +1})$$

$$\displaystyle\geq\sum_{k=1}^{\ell}\Big{(}{\cal I}_{\cal D}(k)c_{1}(\beta)|\lambda_{k}|^{3}+c_{2}(\alpha)\|g(\hat{x}_{k})\|_{2}^{2}\Big{)}.$$

$$\sum _{k=1,k\in \mathcal{D}}^{\infty }|{\lambda }_k{|}^3\equiv \sum _{k=1}^{\infty }|({\lambda }_k{)}_{-}{|}^3$$

$$\text{and}\sum _{k=1}^{\infty }\Vert g({\hat{x}}_k){\Vert }_2^2$$

$$\underset{k\to \infty }{\lim }\Vert g({\hat{x}}_k){\Vert }_2=0,$$

$$\sum _{k=1}^{\infty }\Vert {\hat{x}}_k-x_k{\Vert }_2^3={\beta }^3\sum _{k=1}^{\infty }\Vert d_k{\Vert }_2^3\le {\beta }^3{\theta }^3\sum _{k=1}^{\infty }|({\lambda }_k{)}_{-}{|}^3<\infty ,$$

$$\underset{k\to \infty }{\lim }\Vert {\hat{x}}_k-x_k{\Vert }_2=0.$$

$$0$$

$$=\underset{k\to \infty }{\mathrm{limsup}}\Vert g(x_k)-g({\hat{x}}_k)+g({\hat{x}}_k){\Vert }_2$$

$$\le \underset{k\to \infty }{\mathrm{limsup}}\Vert g(x_k)-g({\hat{x}}_k){\Vert }_2+\underset{k\to \infty }{\mathrm{limsup}}\Vert g({\hat{x}}_k){\Vert }_2$$

$$\le L\underset{k\to \infty }{\mathrm{limsup}}\Vert x_k-{\hat{x}}_k{\Vert }_2+\underset{k\to \infty }{\mathrm{limsup}}\Vert g({\hat{x}}_k){\Vert }_2=0,$$

$$\sum _{k=1}^{\infty }|({\lambda }_k{)}_{-}{|}^3\ge \sum _{k\in \mathcal{K}}|({\lambda }_k{)}_{-}{|}^3\ge \sum _{k\in \mathcal{K}}{ϵ}^3=\infty ,$$

$$-g(x_k{)}^T{s}_k\ge \delta \Vert s_k{\Vert }_2\Vert g(x_k){\Vert }_2.$$

$$m_{s,k}(\alpha )$$

$$\text{and}{m}_{d,k}(\beta )$$

$$f(x_k+\alpha s_k)$$

$$\text{and}f(x_k+\beta d_k)$$

$${\alpha }_k:=\frac{-g(x_k{)}^T{s}_k}{L_k\Vert s_k{\Vert }_2^2}\text{and}{\beta }_k:=\frac{\left(-c_k+\sqrt{c_k^2-2{\sigma }_k\Vert d_k{\Vert }_2^{3}g(x_k{)}^T{d}_k}\right)}{{\sigma }_k\Vert d_k{\Vert }_2^3},$$

$$L_k\le L_{\max }:=\max \{L_1,\rho L\}\text{and}{\sigma }_k\le {\sigma }_{\max }:=\max \{{\sigma }_1,\rho \sigma \}.$$

$$\underset{k\to \infty }{\lim }\Vert g(x_k){\Vert }_2=0\text{and}\underset{k\to \infty }{\mathrm{liminf}}{\lambda }_k\ge 0.$$

$$m_{s,k}({\alpha }_k)=\frac{1}{2L_k}{\left(\frac{g(x_k{)}^T{s}_k}{\Vert s_k{\Vert }_2}\right)}^2\ge \frac{{\delta }^2}{2L_k}\Vert g(x_k){\Vert }_2^2.$$

$$m_{d,k}({\beta }_k)$$

Full text

1 Introduction

There has been a recent surge of interest in solving nonconvex optimization problems. A prime example is the dramatic increase in interest in the training of deep neural networks. Another example is the task of clustering data that arise from the union of low dimensional subspaces. In this setting, the nonconvexity typically results from sophisticated modeling approaches that attempt to accurately capture corruptions in the data elhamifar2009sparse ; JiaRV17 . It is now widely accepted that the design of new methods for solving nonconvex problems (at least locally) is sorely needed.

First consider deterministic optimization problems. For solving such problems, most algorithms designed for minimizing smooth objective functions only ensure convergence to first-order stationarity, i.e., that the gradient of the objective asymptotically vanishes. This characterization is certainly accurate for line search methods, which seek to reduce the objective function by searching along descent directions. Relatively few researchers have designed line search (or other, such as trust region or regularization) algorithms that generate iterates that provably converge to second-order stationarity. The reason for this is three-fold: $(i)$ such methods are more complicated and expensive, necessarily involving the computation of directions of negative curvature when they exist; $(ii)$ methods designed only to achieve first-order stationarity rarely get stuck at saddle points that are first-order, but not second-order stationary pmlr-v49-lee16 ; and $(iii)$ there has not been sufficient evidence showing benefits of integrating directions of negative curvature.

For solving stochastic optimization problems, the methods most commonly invoked are variants of the stochastic gradient (SG) method. During each iteration of SG, a stochastic gradient is computed and a step opposite that direction is taken to obtain the next iterate. Even for nonconvex problems, convergence guarantees (e.g., in expectation or almost surely) for SG methods to first-order stationarity have been established under reasonable assumptions; e.g., see BottCurtNoce16 . In fact, SG and its variants represent the current state-of-the-art for training deep neural networks. As for methods that compute and follow negative curvature directions, it is no surprise that such methods (e.g., liu2017noisy ) have not been used or studied extensively since practical benefits have not even been shown in deterministic optimization.

The main purpose of this paper is to revisit and provide new perspectives on the use of negative curvature directions in deterministic and stochastic optimization. Whereas previous work in deterministic settings has focused on line search and other methods, we focus on methods that attempt to construct upper-bounding models of the objective function. This allows them, at least in theory, to employ fixed stepsizes while ensuring convergence guarantees. For a few instances of such methods of interest, we provide theoretical convergence guarantees and empirical evidence showing that an optimization process can benefit by following negative curvature directions. The fact that our methods might employ fixed stepsizes is also important as it allows us to offer new strategies for stochastic optimization where, e.g., line search strategies are not often viable.

1.1 Contributions

The contributions of our work are the following.

• •

For deterministic optimization, we first provide conditions on descent and negative curvature directions that allow us to guarantee convergence to second-order stationarity with a simple two-step iteration with fixed stepsizes; see §2.1. Using the two-step method as motivation, we then propose a dynamic choice for the direction and stepsize; see §2.2. In particular, the dynamic algorithm makes decisions based on which available step appears to offer a more significant objective reduction. The details of our dynamic strategy represent the main novelty of this paper with respect to deterministic optimization.

• •

We prove convergence rate guarantees for our deterministic optimization methods that provide upper bounds for the numbers of iterations required to achieve (approximate) first- and second-order stationarity; see §2.3. We follow this with a discussion, in §2.4, on the issue of how methods might behave in the neighborhood of so-called strict saddle points, a topic of much interest in the literature.

• •

In §2.5 and §2.6, we discuss different techniques for computing the search directions in our deterministic algorithms and provide the results of numerical experiments showing the benefits of following negative curvature directions, as opposed to only descent directions.

• •

For solving stochastic optimization problems, we propose two methods. Our first method shows that one can maintain the convergence guarantees of a stochastic gradient method by adding an appropriately scaled negative curvature direction for a stochastic Hessian estimate; see §3.1. This approach can be seen as a refinement of that in NeelVilnLeSutsKaisKuraMart15 , which adds noise to each SG step.

• •

Our second method for stochastic optimization is an adaptation of our dynamic (deterministic) method when stochastic gradient and Hessian estimates are involved; see §3.2. Although we are unable to establish a convergence theory for this approach, we do illustrate some gain in performance in neural network training; see §3.3. We view this as a first step in the design of a practical algorithm for stochastic optimization that efficiently exploits negative curvature.

Computing directions of negative curvature carries an added cost that should be taken into consideration when comparing algorithm performance. We remark along with the results of our numerical experiments how these added costs can be worthwhile and/or marginal relative to the other per-iteration costs.

1.2 Prior Related Work

For solving deterministic optimization problems, relatively little research has been directed towards the use of negative curvature directions. Perhaps the first exception is the work in more1979use in which convergence to second-order stationary points is proved using a curvilinear search formed by descent and negative curvature directions. In a similar vein, the work in forsgren1995computing offers similar convergence properties based on using a curvilinear search, although the primary focus in that work is describing how a partial Cholesky factorization of the Hessian matrix could be used to compute descent and negative curvature directions. More recently, a linear combination of descent and negative curvature directions was used in an optimization framework to establish convergence to second-order stationary solutions under loose assumptions GilKR16 ; GilKR16c , and a strategy for how to use descent and negative curvature directions was combined with a backtracking linesearch to provide worst case iteration bounds in royer2017complexity . For further instances of work employing negative curvature directions, see Birgin2001 ; goldfarb1980curvilinear ; doi:10.1080/10556780008805794 ; doi:10.1137/15M1022100 . Importantly, none of the papers above (or any others to the best of our knowledge) have established consistent gains in computational performance as a result of using negative curvature directions.

Another recent trend in the design of deterministic methods for nonconvex optimization is to focus on the ability of an algorithm to escape regions around a saddle point, i.e., a first-order stationary point that is not a minimizer nor a maximizer. A prime example of this trend is the work in PateMokhRibe17 in which a standard type of regularized Newton method is considered. (For a more general presentation of regularized Newton methods of which the method in PateMokhRibe17 is a special case, see, e.g., NoceWrig06 .) While the authors do show a probabilistic convergence result for their method to (approximate) second-order stationarity, their main emphasis is on the number of iterations required to escape neighborhoods of saddle points. This result, similar to those in some (but not all) recent papers discussing the behavior of descent methods when solving nonconvex problems (e.g., see jin2017escape ; pmlr-v49-lee16 ), requires that all saddle points are non-degenerate in the sense that the Hessian of the objective at any such point does not have any zero eigenvalues. In particular, in jin2017escape they show that a perturbed form of gradient descent converges to a second-order stationary point in a number iterations that depends poly-logarithmically on the dimension of the problem. Our convergence theory does not require such a non-degeneracy assumption; see §2.4 for additional discussion comparing our results to others. (See also LeePanPilSimJorRec17 , where it is shown, without a non-degeneracy assumption and for almost all starting points, that certain first-order methods do not converge to saddle points at which the Hessian of the objective has a negative eigenvalue. That said, this work does not prove bounds on the number of iterations such a method might spend near saddle points.)

For solving stochastic optimization problems, there has been little work that focuses explicitly on the use of directions of negative curvature. For two examples, see NIPS2015_5870 ; liu2017noisy . Meanwhile, it has recently become clear that the nonconvex optimization problems used to train deep neural networks have a rich and interesting landscape dauphin2014identifying ; keskar2016large . Instead of using negative curvature to their advantage, modern methods ignore it, introduce random perturbations (e.g., see GeHuanJinYuan15 ), or employ regularized/modified Newton methods that attempt to avoid its potential ill-effects (e.g., dauphin2014identifying ; martens2010deep ). Although such approaches are reasonable, one might intuitively expect better performance if directions of negative curvature are exploited instead of avoided. In this critical respect, the methods that we propose are different from these prior algorithms, with the exception of liu2017noisy .

2 Deterministic Optimization

Consider the unconstrained optimization problem

[TABLE]

where the objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ is twice continuously differentiable and bounded below by a scalar $f_{\inf}\in\mathbb{R}$. We define the gradient function $g:=\nabla f$ and Hessian function $H:=\nabla^{2}f$. We assume that both of these functions are Lipschitz continuous on the path defined by the iterates computed in an algorithm, the gradient function $g$ with Lipschitz constant $L\in(0,\infty)$ and the Hessian function $H$ with Lipschitz constant $\sigma\in(0,\infty)$. Given an invertible matrix $M\in\mathbb{R}^{n\times n}$, its Euclidean norm condition number is written as $\kappa(M)=\|M\|_{2}\|M^{-1}\|_{2}$. Given a scalar $\lambda\in\mathbb{R}$, we define $(\lambda)_{-}:=\min\{0,\lambda\}$.

In the remainder of this section, we present a two-step method that is guaranteed to converge toward second-order stationarity, as well as a dynamic approach that chooses between two types of steps at each point based on lower bounds on objective function decrease. Both algorithms are presented in a generic manner that offers flexibility in the ways the steps are computed.

2.1 Two-Step Method

Our first method alternates between negative curvature and descent steps using fixed stepsizes for each. (Either can be taken first; arbitrarily, we state our algorithm and analysis assuming that one starts with a negative curvature step.) At a given iterate $x_{k}\in\mathbb{R}^{n}$, let $\lambda_{k}$ denote the left-most eigenvalue of $H(x_{k})$. If $\lambda_{k}\geq 0$ (i.e., $H(x_{k})\succeq 0$), the algorithm sets $d_{k}\leftarrow 0$; otherwise, $d_{k}$ is computed so that

[TABLE]

for some $\gamma\in(0,1]$ and $\theta\in(0,\infty)$ that are independent of $k$. A step in this direction is then taken to obtain $\hat{x}_{k}\leftarrow x_{k}+\beta d_{k}$ for some $\beta\in(0,(3\gamma)/(\sigma\theta))$. At this point, if $g(\hat{x}_{k})=0$, then the algorithm sets $\hat{s}_{k}\leftarrow 0$; otherwise, $\hat{s}_{k}$ is computed satisfying

[TABLE]

for some $(\delta,\zeta)\in(0,1]\times(0,1]$ and $\eta\in[1,\infty)$ that are independent of $k$. The iteration ends by taking a step along this direction to obtain $x_{k+1}\leftarrow\hat{x}_{k}+\alpha\hat{s}_{k}\equiv x_{k}+\alpha\hat{s}_{k}+\beta d_{k}$ for some $\alpha\in(0,(2\delta\zeta)/(L\eta^{2}))$. We remark that (2) and (3) are satisfiable; e.g., if $\hat{s}_{k}=-g(\hat{x}_{k})$ and $d_{k}$ is an eigenvector corresponding to the left-most eigenvalue $\lambda_{k}$ scaled so that $g(x_{k})^{T}d_{k}\leq 0$ and $\|d_{k}\|_{2}=|\lambda_{k}|$, then (2) and (3) are satisfied with $\gamma=\theta=\delta=\zeta=\eta=1$. For further remarks on step computation techniques, see §2.5.

We now show that Algorithm 1 converges toward second-order stationarity. Critical for this analysis are bounds on $(\alpha,\beta)$ that we have stated (and are also stated in Algorithm 1). Also, for the analysis, it is convenient to define

[TABLE]

along with the indicator ${\cal I}_{\cal D}(k)$, which evaluates as 1 if $k\in{\cal D}$ and as 0 otherwise.

Theorem 2.1

If Algorithm 1 terminates finitely in iteration $k\in\mathbb{N}_{+}$, then $g(x_{k})=0$ and $\lambda_{k}\geq 0$, i.e., $x_{k}$ is second-order stationary. Otherwise, the computed iterates satisfy

[TABLE]

Proof

Algorithm 1 terminates finitely if and only if, for some $k\in\mathbb{N}_{+}$, $d_{k}=\hat{s}_{k}=0$. We can then observe that $d_{k}=0$ if and only if $\lambda_{k}\geq 0$, while $\hat{s}_{k}=0$ if and only if $0=g(\hat{x}_{k})=g(x_{k})$, where the last equality holds because $\hat{x}_{k}=x_{k}$ whenever $d_{k}=0$. These represent the desired conclusions for this case. Otherwise, if Algorithm 1 does not terminate finitely, consider arbitrary $k\in\mathbb{N}_{+}$. If $k\notin{\cal D}$, then $d_{k}=0$ and $\hat{x}_{k}=x_{k}$, meaning that $f(\hat{x}_{k})=f(x_{k})$. Otherwise, if $k\in{\cal D}$, then $d_{k}\neq 0$ and $\lambda_{k}<0$, and, by the step computation conditions in (2), it follows that

[TABLE]

where c_{1}(\beta):=\tfrac{1}{2}\beta^{2}\theta^{2}\big{(}\gamma-\tfrac{1}{3}\sigma\beta\theta\big{)}\in(0,\infty). Similarly, it follows from (3) that

[TABLE]

where $c_{2}(\alpha):=\alpha(\delta\zeta-\tfrac{1}{2}L\alpha\eta^{2})\in(0,\infty)$. Overall,

[TABLE]

Now observe that, for any $\ell\in\mathbb{N}_{+}$, it follows that

[TABLE]

This inequality and $f$ being bounded below show that

[TABLE]

The latter bound yields

[TABLE]

while the former bound and (2c) yield

[TABLE]

from which it follows that

[TABLE]

It follows from Lipschitz continuity of $g$ along with (7) and (8) that

[TABLE]

which implies the first limit in (4). Finally, in order to derive a contradiction, suppose that $\liminf_{k\to\infty}\lambda_{k}<0$, meaning that there exists some $\epsilon>0$ and infinite index set ${\cal K}\subseteq{\cal D}$ such that $\lambda_{k}\leq-\epsilon$ for all $k\in{\cal K}$. This implies that

[TABLE]

contradicting (6a). This yields the second limit in (4). ∎

There are two potential weaknesses of this two-step approach. First, it simply alternates back-and-forth between descent and negative curvature directions, which might not always lead to the most productive step from each point. Second, even though our analysis holds for all stepsizes $\alpha$ and $\beta$ in the intervals provided in Algorithm 1, the algorithm might suffer from poor performance if these values are chosen poorly. We next present a method that addresses these weaknesses.

2.2 Dynamic Method

Suppose that, in any iteration $k\in\mathbb{N}_{+}$ when $\lambda_{k}<0$, one computes a nonzero direction of negative curvature satisfying (2a)–(2b) for some $\gamma\in(0,1]$. Suppose also that, if $g(x_{k})\neq 0$, then one computes a nonzero direction $s_{k}$ satisfying the equivalent of the first condition in (3), namely, for some $\delta\in(0,1]$,

[TABLE]

Defining, for $(L_{k},\sigma_{k})\in(0,\infty)\times(0,\infty)$, the model reductions

[TABLE]

ones finds that, if $L_{k}\geq L$ and $\sigma_{k}\geq\sigma$, then

[TABLE]

These two inequalities suggest that, during iteration $k\in\mathbb{N}_{+}$, one could choose which of the two steps ($s_{k}$ or $d_{k}$) to take based on which model reduction predicts the larger decrease in the objective. One can verify that the reductions $m_{s,k}$ and $m_{d,k}$ are maximized (over the positive real numbers) by

[TABLE]

where $c_{k}:=d_{k}^{T}H(x_{k})d_{k}$ is a measure of $H(x_{k})$-curvature for $d_{k}$.

Algorithm 2, stated below, follows this dynamic strategy of choosing between $s_{k}$ and $d_{k}$ for all $k\in\mathbb{N}_{+}$. It also involves dynamic updates for Lipschitz constant estimates, represented in iteration $k\in\mathbb{N}_{+}$ by $L_{k}$ and $\sigma_{k}$. In this deterministic setting, a step is only taken if it yields an objective function decrease. Otherwise, a null step is effectively taken and a Lipschitz constant estimate is increased.

In the next two results, we establish that Algorithm 2 is well-defined and that it has convergence guarantees on par with Algorithm 1 (recall Theorem 2.1).

Lemma 1

Algorithm 2 is well defined in the sense that it either terminates finitely or generates infinitely many iterates. In addition, at the end of each iteration $k\in\mathbb{N}_{+}$,

[TABLE]

Proof

We begin by showing that the loop in Step 6 terminates finitely anytime it is entered. For a proof by contradiction, if that loop were never to terminate, then the updates to $L_{k}$ and/or $\sigma_{k}$ would cause at least one of them to become arbitrarily large. Since (10a) holds whenever $L_{k}\geq L$ and (10b) holds whenever $\sigma_{k}\geq\sigma$, it follows that the loop would eventually terminate, thus reaching a contradiction. Therefore, we have proved that the loop in Step 6 terminates finitely anytime it is entered, and moreover that (12) holds. This completes the proof once we observe that during iteration $k\in\mathbb{N}_{+}$, Algorithm 2 might finitely terminate in Step 5. ∎

Theorem 2.2

If Algorithm 2 terminates finitely in iteration $k\in\mathbb{N}_{+}$, then $g(x_{k})=0$ and $\lambda_{k}\geq 0$, i.e., $x_{k}$ is second-order stationary. Otherwise, the computed iterates satisfy

[TABLE]

Proof

Algorithm 2 terminates finitely only if, for some $k\in\mathbb{N}_{+}$, $d_{k}=s_{k}=0$. This can only occur if $\lambda_{k}\geq 0$ and $g(x_{k})=0$, which are the desired conclusions. Otherwise, Algorithm 2 does not terminate finitely and during each iteration at least one of $s_{k}$ and $d_{k}$ is nonzero. We use this fact in the following arguments.

Consider arbitrary $k\in\mathbb{N}_{+}$. If $s_{k}\neq 0$, then the definition of $\alpha_{k}$ in (11) and the step computation condition on $s_{k}$ in (9) ensure that

[TABLE]

Similarly, if $d_{k}\neq 0$, then by the fact that $\beta_{k}$ maximizes $m_{d,k}(\beta)$ over $\beta>0$, (2a)–(2b), and defining $\hat{\beta}_{k}:=-2d_{k}^{T}H(x_{k})d_{k}/(\sigma_{k}\|d_{k}\|_{2}^{3})>0$, one finds

[TABLE]

Overall, for all $k\in\mathbb{N}_{+}$, since at least one of $s_{k}$ and $d_{k}$ is nonzero, $\|g_{k}\|_{2}=0$ if and only if $\|s_{k}\|_{2}=0$, and $|(\lambda_{k})_{-}|=0$ if and only if $\|d_{k}\|_{2}=0$, it follows that

[TABLE]

Indeed, to show that (16) holds, let us consider two cases. First, suppose that the update $x_{k+1}\leftarrow x_{k}+\alpha_{k}s_{k}$ is completed, meaning that (10a) holds with $\alpha=\alpha_{k}$ and $m_{s,k}(\alpha_{k})\geq m_{d,k}(\beta_{k})$. Combining these facts with (14) and (15) establishes (16) in this case. Second, suppose that $x_{k+1}\leftarrow x_{k}+\beta_{k}d_{k}$ is completed, meaning that (10b) holds with $\beta=\beta_{k}$ and $m_{s,k}(\alpha_{k})<m_{d,k}(\beta_{k})$. Combining these facts with (14) and (15) establishes (16) in this case. Thus, (16) holds for all $k\in\mathbb{N}_{+}$.

It now follows from (16), the bounds in (12), and a proof similar to that used in Theorem 2.1 (in particular, to establish (6a) and (6b)) that

[TABLE]

One may now establish the desired results in (13) using the same arguments as used in the proof of Theorem 2.1. ∎

Let us add a few remarks about Algorithm 2.

• •

Our convergence theory allows $L_{k+1}\leftarrow L_{k}$ and $\sigma_{k+1}\leftarrow\sigma_{k}$ in Step 24 for all $k\in\mathbb{N}_{+}$, in which case the Lipschitz constant estimates are monotonically increasing. However, in Algorithm 2, we allow these estimates to decrease since this might yield better results in practice.

• •

If the Lipschitz constants $L$ and $\sigma$ for $g:=\nabla f$ and $H:=\nabla^{2}f$, respectively, are known, then one could simply set $L_{k}=L$ and $\sigma_{k}=\sigma$ during each iteration; in this case, the loop would not actually be needed. Although this would simplify the presentation, it would generally result in more iterations being required to obtain (approximate) first- and/or second-order stationarity. However, if the cost of evaluating $f$ is substantial, then static parameters might work well.

• •

Each time through the loop, condition (10a) or (10b) is tested, but not both, since this would require an extra evaluation of $f$. If the cost of evaluating the objective function is not a concern, then one could choose between the two steps based on actual objective function decrease rather than model decrease.

2.3 Complexity Analysis

We have the following complexity result for Algorithm 2. We claim that a similar result could be stated and proved for Algorithm 1 as well, where one employs the inequality (5) in place of (16) in the proof below. However, for brevity and since we believe it might often be the better method in practice, we focus on the following result for our dynamic method, Algorithm 2.

Theorem 2.3

Consider any scalars $\epsilon_{g}\in(0,\bar{\epsilon}_{g}]$ and $\epsilon_{H}\in(0,\bar{\epsilon}_{H}]$ for some constants $(\bar{\epsilon}_{g},\bar{\epsilon}_{H})\in(0,\infty)\times(0,\infty)$. With respect to Algorithm 2, the cardinality of the index set

[TABLE]

is at most ${\cal O}(\epsilon_{g}^{-2})$. In addition, the cardinality of the index set

[TABLE]

is at most ${\cal O}(\epsilon_{H}^{-3})$. Hence, the number of iterations and derivative (i.e., gradient and Hessian) evaluations required until iteration $k\in\mathbb{N}_{+}$ is reached with

[TABLE]

is at most ${\cal O}(\max\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\})$. Moreover, if Step 24 is modified such that

[TABLE]

for some $(L_{\min},\sigma_{\min})\in(0,\infty)\times(0,\infty)$ with $L_{1}\geq L_{\min}$ and $\sigma_{1}\geq\sigma_{\min}$, then the number of iterations in the loop for all $k\in\mathbb{N}_{+}$ is uniformly bounded, meaning that the complexity bound above also holds for the number of function evaluations.

Proof

As in the proof of Theorem 2.2, one has that, for all $k\in\mathbb{N}_{+}$, the inequality (16) holds, which we restate here as

[TABLE]

From this inequality and the bounds in (12), it follows that

[TABLE]

Since $f$ is bounded below by $f_{\inf}$ and (18) ensures that $\{f(x_{k})\}$ is monotonically decreasing, the inequalities above imply that ${\cal G}(\epsilon_{g})$ and ${\cal H}(\epsilon_{H})$ are both finite. In addition, by summing the reductions achieved in $f$ over the iterations in the above index sets, we obtain

[TABLE]

Rearranging, one obtains that

[TABLE]

as desired. Finally, the last conclusion follows by the fact that, with (17), the number of iterations within the loop for any $k\in\mathbb{N}_{+}$ is bounded above by

[TABLE]

which is the maximum number of updates to $L_{k}$ and/or $\sigma_{k}$ to bring values within $[L_{\min},L]$ and $[\sigma_{\min},\sigma]$ up to satisfy $L_{k}\geq L$ and $\sigma_{k}\geq\sigma$. ∎

2.4 Behavior Near Strict Saddle Points

As previously mentioned (recall §1.2), much recent attention has been directed toward the behavior of nonconvex optimization algorithms in the neighborhood of saddle points. In particular, in order to prove guarantees about avoiding saddle points, an assumption is often made about all saddle points being strict (or ridable) or even non-degenerate. Strict saddle points are those at which the Hessian has at least one negative eigenvalue; intuitively, these are saddle points from which an algorithm should usually be expected to avoid. Nondegenerate saddle points are ones at which the Hessian has no zero eigenvalues, a much stronger assumption.

One can show that in certain problems of interest, all saddle points are strict. This is interesting. However, much of the work that has discussed the behavior of nonconvex algorithms in the neighborhood of strict saddle points have focused on standard types of descent methods, about which one can only prove high probability results LeePanPilSimJorRec17 . By contrast, when one considers explicitly computing directions of negative curvature, this is all less of an issue. After all, notice that our convergence and complexity analyses for our framework did not require careful consideration of the nature of any potential saddle points.

In any case, for ease of comparison to other recent analyses, let us discuss the convergence/complexity properties for our method in the context of a strict saddle point assumption. Suppose that the set of maximizers and saddle points of $f$, say $\{\bar{x}_{i}\}_{i\in\cal{I}}$ for some index set $\cal{I}$, which must necessarily have $g(\mkern 1.5mu\overline{\mkern-1.5mux}_{i})=0$ for all $i\in{\cal I}$, also has $\bar{\lambda}_{i}<0$ for all $i\in{\cal I}$, where $\{\bar{\lambda}_{i}\}_{i\in{\cal I}}$ are the leftmost eigenvalues of $\{H(\bar{x}_{i})\}_{i\in{\cal I}}$ and are uniformly bounded away from zero. In this setting, we may draw the following conclusions from Theorems 2.2 and 2.3.

• •

Any limit point of the sequence $\{x_{k}\}$ computed by Algorithm 2 is a minimizer of $f$. This follows since Theorem 2.2 ensures that for any limit point, say $\mkern 1.5mu\overline{\mkern-1.5mux}$, the gradient must be zero and the leftmost eigenvalue must be nonnegative. It follows from these facts and the assumptions on the maximizers and saddle points $\{\mkern 1.5mu\overline{\mkern-1.5mux}_{i}\}$ that $\mkern 1.5mu\overline{\mkern-1.5mux}$ must be a minimizer, as claimed.

• •

From the discussion in the previous bullet and Theorem 2.3 we know, in fact, that the iterates of Algorithm 2 must eventually enter a region consisting only of minimizers (i.e., one that does not contain any maximizers or saddle points) in a number of iterations that is polynomial in $\epsilon\in(0,\infty)$, where the negative value $-\epsilon$ is greater than the largest of the leftmost eigenvalues of $\{H(\mkern 1.5mu\overline{\mkern-1.5mux}_{i})\}_{i\in{\cal I}}$. This complexity result holds without assuming that the minimizers, maximizers, and saddle points are all nondegenerate (i.e., the Hessian matrix at these points are nonsingular), which, e.g., should be contrasted with the analysis in (PateMokhRibe17, , see Assumption 3). The primary reason for our stronger convergence/complexity properties is that our method incorporates negative curvature directions when they exist, as opposed to only descent directions.

2.5 Step Computation Techniques

There is flexibility in the ways in which the steps $d_{k}$ and $s_{k}$ (or $\hat{s}_{k}$) are computed in order to satisfy the desired conditions (in (2) and (3)/(9)). For example, $s_{k}$ might be the steepest descent direction $-g(x_{k})$, for which (3) holds with $\delta=\zeta=\eta=1$. Another option, with a symmetric positive definite $B_{k}^{-1}\in\mathbb{R}^{n\times n}$ that has $\kappa(B_{k}^{-1})\leq\delta^{-1}$ with a spectrum falling in $[\zeta,\eta]$, is to compute $s_{k}$ as the modified Newton direction $-B_{k}^{-1}g(x_{k})$. A particularly attractive option for certain applications (when Hessian-vector products are readily computed) is to compute $s_{k}$ via a Newton-CG routine NoceWrig06 with safeguards that terminate the iteration before a CG iterate is computed that violates (3) for prescribed $(\delta,\zeta,\eta)$.

There are multiple ways to compute the negative curvature direction. In theory, the most straightforward approach is to set $d_{k}=\pm v_{k}$, where $v_{k}$ is a leftmost eigenvector of $H(x_{k})$. With this choice, it follows that $d_{k}^{T}H(x_{k})d_{k}=\lambda_{k}\|d_{k}\|^{2}$, meaning that (2a) is satisfied, and one can scale $d_{k}$ to ensure that (2b) and (2c) hold. A second approach for large-scale settings (i.e., when $n$ is large) is to compute $d_{k}$ via matrix-free Lanczos iterations lanczos1950iteration . Such an approach can produce a direction $d_{k}$ satisfying (2a), which can then be scaled to yield (2b) and (2c).

2.6 Numerical Results

In this section, we demonstrate that there can be practical benefits of following directions of negative curvature if one follows the dynamic approach of Algorithm 2. To do this, we implemented software in Matlab that, for all $k\in\mathbb{N}_{+}$, has the option to compute $s_{k}$ via several options (see §2.6.1 and §2.6.2) and $d_{k}=\pm v_{k}$ (recall §2.5). Our test problems include a subset of the CUTEst collection gould2015cutest . Specifically, we selected all of the unconstrained problems with $n\leq 500$ and second derivatives explicitly available. This left us with a test set of $97$ problems.

Using this test set, we considered two variants of Algorithm 2: $(i)$ a version in which the if condition in Step 7 is always presumed to test true, i.e., the descent step $s_{k}$ is chosen for all $k\in\mathbb{N}_{+}$ (which we refer to as Algorithm 2($s_{k}$)), and $(ii)$ a version that, as in our formal statement of the algorithm, chooses between descent and negative curvature steps by comparing model reduction values for each $k\in\mathbb{N}_{+}$ (which we refer to as Algorithm 2($s_{k},d_{k}$)). In our experiments, we used $L_{1}\leftarrow 1$ and $\sigma_{1}\leftarrow 1$, updating them and setting subsequent values in their respective sequences using the following strategy. For increasing one of these values in Step 13 or Step 20, we respectively set the quantities

[TABLE]

then, with $\rho\leftarrow 2$, use the update

[TABLE]

The quantity $\hat{L}_{k}$ is defined so that (10a) holds at equality with $\alpha=\alpha_{k}$ and $L_{k}$ replaced by $\hat{L}_{k}$ in the definition of $m_{s,k}(\alpha_{k})$, i.e., $\hat{L}_{k}$ is the value that makes the model decrease agree with the exact function decrease at $x_{k}+\alpha_{k}s_{k}$. The procedure we use to set $\hat{\sigma}_{k}$ is analogous. We remark that the updates in (19) ensure that $L_{k}\in[\rho L_{k},10^{3}L_{k}]$ and $\sigma_{k}\in[\rho\sigma_{k},10^{3}\sigma_{k}]$, which we claim maintains the convergence and complexity guarantees of Algorithm 2. Moreover, when Step 13 (respectively, Step 20) is reached, then it must be the case that $\hat{L}_{k}>L_{k}$ (respectively, $\hat{\sigma}_{k}>\sigma_{k}$). On the other hand, in Step 24 we use the following updates:

[TABLE]

Over the next two sections, we discuss numerical results when two different options for computing the descent direction $s_{k}$ are used.

2.6.1 Choosing $s_{k}$ as the direction of steepest descent

The tests in this section use the steepest descent direction, i.e., $s_{k}=-g(x_{k})$ for all $k\in\mathbb{N}_{+}$. Although this is a simple choice, it gives a starting point for understanding the potential benefits of using directions of negative curvature.

We ran Algorithm 2($s_{k}$) and Algorithm 2($s_{k},d_{k}$) on the previously described CUTEst problems with commonly used stopping conditions. Specifically, an algorithm terminates with an approximate second-order stationary solution if

[TABLE]

We also terminate an algorithm if an iteration limit of $10,\!000$ is reached or if a trial step smaller than $10^{-16}$ is computed. The results can be found in Figure 1.

In Figure 1(a), letting $f_{\text{final}}(s_{k})$ and $f_{\text{final}}(s_{k},d_{k})$ be the final computed objective values for Algorithm 2($s_{k}$) and Algorithm 2($s_{k},d_{k}$), respectively, we plot

[TABLE]

for each problem that Algorithm 2($s_{k},d_{k}$) used at least one negative curvature direction. In this manner, an upward pointing bar implies that Algorithm 2($s_{k},d_{k}$) terminated with a lower value of the objective function, while its magnitude represents how much better was this value. We can observe from Figure 1(a) that among the $39$ problems, Algorithm 2($s_{k}$) terminated with a significantly lower objective value compared to Algorithm 2($s_{k},d_{k}$) only $2$ times. (While the value in (21) falls in the interval $[-2,2]$, the bars in the figure all fall within $[-1,1]$ since, for each problem, the final function value reached by each algorithm had the same sign.)

We are also interested in the number of iterations and objective function evaluations needed to obtain these final objective function values. The relative performances of the algorithms with respect to these measures are shown in Figure 1(b) and Figure 1(c). In particular, letting $\#its({s_{k}})$ be the total number of iterations required by Algorithm 2($s_{k}$) and $\#its({s_{k},d_{k}})$ be the total number of iterations required by Algorithm 2($s_{k},d_{k}$), we plot in Figure 1(b) the values

[TABLE]

and, letting $\#fevals({s_{k}})$ be the total number of objective function evaluations for Algorithm 2($s_{k}$) and $\#fevals({s_{k},d_{k}})$ be the total number of objective function evaluations for Algorithm 2($s_{k},d_{k}$), we plot in Figure 1(c) the values

[TABLE]

These two plots show that Algorithm 2($s_{k},d_{k}$) tends to perform better than Algorithm 2($s_{k}$) in terms of both the number of required iterations and function evaluations. Overall, we find these results interesting since Algorithm 2($s_{k},d_{k}$) tends to find lower values of the objective function (see Figure 1(a)) while typically requiring fewer iterations (see Figure 1(b)) and function evaluations (see Figure 1(c)).

2.6.2 Choosing $s_{k}$ using a modified-Newton strategy

In this subsection, we show the results of tests similar to those in §2.6.1 except that now we compute the descent direction $s_{k}$ by a modified-Newton approach. Specifically, we compute $s_{k}$ as the unique vector satisfying $B_{k}s_{k}=-g(x_{k})$ with

[TABLE]

where $I$ is the identify matrix and $\delta_{k}$ is the smallest nonnegative real number such that $B_{k}$ is positive definite with a condition number less than or equal to $10^{8}$.

The results from solving our CUTEst test set with this choice of $B_{k}$ using the stopping condition (20) are presented in Figure 2. We can observe that, overall, Algorithm 2($s_{k},d_{k}$) outperforms Algorithm 2($s_{k}$) in this experiment. Neither algorithm consistently outforms the other in terms of their final objective values, but Algorithm 2($s_{k},d_{k}$) typically requires fewer iterations (see Figure 2(b)) and objective function evaluations (see Figure 2(c)).

3 Stochastic Optimization

Let us now consider the problem to minimize a twice continuously differentiable and bounded below (by $f_{\inf}\in\mathbb{R}$) objective function $f:\mathbb{R}^{n}\to\mathbb{R}$ defined by the expectation, in terms of the distribution of a random variable $\xi$ with domain $\Xi$, of a stochastic function $F:\mathbb{R}^{n}\times\Xi\to\mathbb{R}$, namely,

[TABLE]

In this context, we expect that, at an iterate $x_{k}$, one can only compute stochastic gradient and Hessian estimates. We do not claim that we are able to prove convergence guarantees to second-order stationarity as in the deterministic case. That said, we are able to present a two-step method with convergence guarantees to first-order stationarity whose structure motivates a dynamic method that we show can offer beneficial practical performance by exploring negative curvature.

3.1 Two-Step Method: Stochastic Gradient/Newton with “Curvature Noise”

At an iterate $x_{k}$, let $\xi_{k}$ be a random variable representing a seed for generating a vector $s_{k}\in\mathbb{R}^{n}$. For example, if $f$ is the expected function value over inputs from a dataset, then $\xi_{k}$ might represent sets of points randomly drawn from the dataset. With $\mathbb{E}_{\xi_{k}}[\cdot]$ denoting expectation taken with respect to the distribution of $\xi_{k}$ given the current iterate $x_{k}$, we require the vector $s_{k}$ to satisfy

[TABLE]

for some $\delta\in(0,1]$, $\eta\in[1,\infty)$, and $(M_{s1},M_{s2})\in(0,\infty)\times(1,\infty)$ that are all independent of $k$. For example, as in a stochastic gradient method, these conditions are satisfied if $s_{k}$ is an unbiased estimate of $\nabla f(x_{k})$ with second moment bounded as in (26c). They are also satisfied in the context of a stochastic Newton method wherein a stochastic gradient estimate is multiplied by a stochastic inverse Hessian estimate, assuming that the latter is conditionally uncorrelated with the former and has eigenvalues contained within an interval of the positive real line uniformly over all $k\in\mathbb{N}_{+}$.

Let us also define $\xi_{k}^{H}$, conditionally uncorrelated with $\xi_{k}$ given $x_{k}$, as a random variable representing a seed for generating an unbiased Hessian estimate $H_{k}$ such that $\mathbb{E}_{\xi_{k}^{H}}[H_{k}]=\nabla^{2}f(x_{k})$. We use $H_{k}$ to compute a direction $d_{k}$. For the purpose of ideally following a direction of negative curvature (for the true Hessian), we ask that $d_{k}$ satisfies a curvature condition similar to that used in the deterministic setting. Importantly, however, the second moment of $d_{k}$ must be bounded similar to that of $s_{k}$ above. Overall, with $\lambda_{k}$ being the left-most eigenvalue of $H_{k}$, we set $d_{k}\leftarrow 0$ if $\lambda_{k}\geq 0$, and otherwise require the direction $d_{k}$ to satisfy

[TABLE]

for some $\gamma\in(0,1]$ and $(M_{d1},M_{d2})\in(0,\infty)\times(1,\infty)$ that are independent of $k$. One manner in which these conditions can be satisfied is to compute $d_{k}$ as an eigenvector corresponding to the left-most eigenvalue $\lambda_{k}$, scaled such that $\|d_{k}\|_{2}^{2}$ is bounded above in proportion to the squared norm $\|s_{k}\|_{2}^{2}$ where $s_{k}$ satisfies (26).

Note that our conditions in (27) do not involve expected descent with respect to the true gradient at $x_{k}$. This can be viewed in contrast to (2), which involves (2b). The reason for this is that, in a practical setting, such a condition might not be verifiable without computing the exact gradient explicitly, which might be intractable or prohibitively expensive. Instead, without this restriction, a critical component of our algorithm is the generation of an independent random scalar $\omega_{k}$ uniformly distributed in $[-1,1]$. With this choice, one finds $\mathbb{E}_{(\xi_{k}^{H},\omega_{k})}[\omega_{k}d_{k}]=0$ such that, effectively, $\omega_{k}d_{k}$ merely adds noise to $s_{k}$ in a manner that ideally follows a negative curvature direction. This leads to Algorithm 3 below.

Algorithm 3 maintains the convergence guarantees of a standard stochastic gradient (SG) method. To show this, let us first prove the following lemma.

Lemma 2

For all $k\in\mathbb{N}_{+}$, it follows that

[TABLE]

Proof

From Lipschitz continuity of $\nabla f$, it follows that

[TABLE]

Taking expectations with respect to the distribution of the random quantities $(\xi_{k},\xi_{k}^{H},\omega_{k})$ given $x_{k}$ and using (26)–(27), it follows that

[TABLE]

Rearranging, we reach the desired conclusion. ∎

From this lemma, we obtain a critical bound similar to one that can be shown for a generic SG method; e.g., see Lemma 4.4 in BottCurtNoce16 . Hence, following analyses such as that in BottCurtNoce16 , one can show that the total expectation for the gap between $f(x_{k})$ and a lower bound for $f$ decreases. For instance, we prove the following theorem using similar proof techniques as for Theorems 4.8 and 4.9 in BottCurtNoce16 .

Theorem 3.1

Suppose that Algorithm 3 is run with $\alpha_{k}=\beta_{k}=\bar{\alpha}$ for all $k\in\mathbb{N}_{+}$ where

[TABLE]

Then, for all $K\in\mathbb{N}_{+}$, one has that

[TABLE]

On the other hand, if Algorithm 3 is run with $\{\alpha_{k}\}$ satisfying

[TABLE]

and $\{\beta_{k}\}=\{\chi\alpha_{k}\}$ for some $\chi\in(0,\infty)$, then it holds that

[TABLE]

from which it follows, with $A_{K}:=\sum_{k=1}^{K}\alpha_{k}$, that

[TABLE]

which implies that $\liminf_{k\to\infty}\mathbb{E}[\|\nabla f(x_{k})\|_{2}^{2}]=0$.

Proof

First, suppose that Algorithm 3 is run with $\alpha_{k}=\beta_{k}=\bar{\alpha}$ for all $k\in\mathbb{N}_{+}$. Then, taking total expectation in (28), it follows with (29) and since $\tfrac{1}{2}>\tfrac{1}{6}$ that

[TABLE]

Summing both sides of this inequality for $k\in\{1,\dots,K\}$ yields

[TABLE]

Rearranging and dividing further by $K$ yields the desired conclusion.

Now suppose that Algorithm 3 is run with $\{\alpha_{k}\}$ and $\{\beta_{k}\}=\{\chi\alpha_{k}\}$ such that the former satisfies (30). Since the second condition in (30) ensures that $\{\alpha_{k}\}\searrow 0$, we may assume without loss of generality that, for all $k\in\mathbb{N}_{+}$,

[TABLE]

Thus, taking total expectation in (28) leads to

[TABLE]

Summing both sides for $k\in\{1,\dots,K\}$ yields

[TABLE]

from which it follows that

[TABLE]

The second of the conditions in (30) implies that the right-hand side here converges to a finite limit when $K\to\infty$. Then, the rest of the desired conclusion follows since the first of the conditions in (30) ensures that $A_{K}\to\infty$ as $K\to\infty$. ∎

3.2 Dynamic Method

Borrowing ideas from the two-step deterministic method in §2.1, the dynamic deterministic method in §2.2, and the two-step stochastic method in §3.1, we propose the dynamic method for stochastic optimization presented as Algorithm 4 below. After computing a stochastic Hessian estimate $H_{k}\in\mathbb{R}^{n\times n}$ and an independent stochastic gradient estimate $g_{k}\in\mathbb{R}^{n}$, the algorithm employs the conjugate gradient (CG) method NoceWrig06 for solving the linear system $H_{k}s=-g_{k}$ with the starting point $s\leftarrow 0$. If $H_{k}\succ 0$, then it is well known that this method solves this system in at most $n$ iterations (in exact arithmetic). However, if $H_{k}\not\succ 0$, then the method might encounter a direction of nonpositive curvature, say $p\in\mathbb{R}^{n}$, such that $p^{T}H_{k}p\leq 0$. If this occurs, then we terminate CG immediately and set $d_{k}\leftarrow p$. This choice is made rather than spend any extra computational effort attempting to approximate an eigenvector corresponding to the left-most eigenvalue of $H_{k}$. Otherwise, if no such direction of nonpositive curvature is encountered, then the algorithm chooses $d_{k}\leftarrow 0$. In either case, the algorithm sets $s_{k}$ as the final CG iterate computed prior to termination. (The only special case is when nonpositive curvature is encountered in the first iteration; then, $s_{k}\leftarrow-g_{k}$ and $d_{k}\leftarrow 0$.)

For setting values for the dynamic parameters $\{L_{k}\}$ and $\{\sigma_{k}\}$, which in turn determine the stepsize sequences $\{\alpha_{k}\}$ and $\{\beta_{k}\}$, the algorithm employs stochastic function value estimates. In this manner, the algorithm can avoid computing the exact objective value at any point (which is often not tractable). In the statement of the algorithm, we use $f_{k}:\mathbb{R}^{n}\to\mathbb{R}$ to indicate a function that yields a stochastic function estimate during iteration $k\in\mathbb{N}_{+}$.

As previously mentioned, we do not claim convergence guarantees for this method. However, we believe that it is well motivated by our previous algorithms. One should also note that the per-iteration cost between this method and an inexact Newton-CG method is negligible since any computed $d_{k}\neq 0$ comes essentially for free from the CG routine. The only significant extra cost might come from the stochastic function estimates, though these can be made cheaper than any stochastic gradient estimate and might even be ignored completely if one is able to tune fixed values for $L$ and $\sigma$ that work well for a particular application.

3.3 Numerical Experiments

We implemented Algorithm 4 in Python 2.7.13. As a test problem, we trained a convolutional neural network to classify handwritten digits in the well-known mnist dataset; see LeCuBottBengHaff98 . Our neural network, implemented using tensorflow,111https://www.tensorflow.org/ is composed of two convolutional layers followed by a fully connected layer. In each iteration, we computed stochastic gradient and Hessian estimates using independently drawn mini-batches of size 500. By contrast, the entire training dataset involves 60,000 feature/label pairs. The testing set involves 10,000 feature/label pairs.

In each iteration, our implementation runs the CG method for at most 10 iterations. If a direction of nonpositive curvature is found within this number of iterations, then it is employed as $d_{k}$; otherwise, $d_{k}\leftarrow 0$. In preliminary training runs, we occasionally witnessed instances in which a particularly large stochastic gradient led to a large step, which in turn spoiled all previous progress made by the algorithm. Hence, to control the iterate displacements, we scaled $s_{k}$ and/or $d_{k}$ down, when necessary, to ensure that $\|s_{k}\|\leq 10$ and $\|\beta_{k}d_{k}\|\leq 0.2\|\alpha_{k}s_{k}\|$. We witnessed similar poor behavior when the Lipschitz constant estimates were initialized to be too small; hence, for our experiments, we chose $L_{1}\leftarrow 80$ and $\sigma_{1}\leftarrow 100$ as the smallest values that yielded stable algorithm behavior for both algorithms. If stochastic function evaluations suggested an objective increase, then an estimate was increased by a factor of 1.2; see Steps 8 and 15 in Algorithm 4. The implementation never decreases these estimates. In addition, while the implementation always takes the step $\alpha_{k}s_{k}$ (i.e., it does not follow the optional Step 9), it does reset $x_{k+1}\leftarrow\hat{x}_{k}$ (recall Step 16) if/when the stochastic function estimates predict an increase in $f$ due to the step $\beta_{k}d_{k}$.

For comparison purposes, we ran the algorithm twice using the same starting point and initial seeds for the random number generators: once with $\beta_{k}$ being reset to zero for all $k\in\mathbb{N}_{+}$ (so that no step along $d_{k}$ is ever taken) and once with it set as in Algorithm 4. We refer to the former algorithm as SG since it is a stochastic-gradient-type method that does not explore negative curvature. We refer to the latter as NC since it attempts to exploit negative curvature.

The training losses as a function of the iteration counter are shown in Figure 3. As can be seen in the plot, the performance of the two algorithms is initially very similar. However, after some initial iterations, following the negative curvature steps consistently offers additional progress, allowing NC to reduce the loss and increase both the training and testing accuracy more rapidly than SG. Eventually, the plots in each of the figures near each other (after around two epochs, when the algorithms were terminated). This should be expected as both algorithms eventually near stationary points. However, prior to this point, NC has successfully avoided the early stagnation experienced by SG.

4 Conclusion

We have confronted the question of whether it can be beneficial for nonconvex optimization algorithms to compute and explore directions of negative curvature. We have proposed new algorithmic frameworks based on the idea that an algorithm might alternate or choose between descent and negative curvature steps based on properties of upper-bounding models of the objective function. In the case of deterministic optimization, we have shown that our frameworks possess convergence and competitive complexity guarantees in the pursuit of first- and second-order stationary points, and have demonstrated that instances of our framework outperform descent-step-only methods in terms of finding points with lower objective values typically within fewer iterations and function evaluations. In the case of stochastic optimization, we have shown that an algorithm that employs “curvature noise” can outperform a stochastic-gradient-based approach.