Complexity analysis of second-order line-search algorithms for smooth nonconvex optimization

Introduction

We consider the unconstrained optimization problem

where $\epsilon_{g},\epsilon_{H}\in(0,1)$ are (typically small) prescribed tolerances. Numerous algorithms have been proposed in recent years for finding points that satisfy (2), each with a complexity guarantee, which is an upper bound on an index $k$ that satisfies (2), in terms of $\epsilon_{g}$, $\epsilon_{H}$, and other quantities. We summarize below the main results.

Classical second-order convergent trust-region schemes can be shown to satisfy (2) after at most $\mathcal{O}\left(\max\left\{\epsilon_{g}^{-2}\,\epsilon_{H}^{-1},\epsilon_{H}^{-3}\right\}\right)$ iterations . Cubic regularization methods in their basic form have better complexity bounds than trust-region schemes, requiring at most $\mathcal{O}\left(\max\left\{\epsilon_{g}^{-2},\epsilon_{H}^{-3}\right\}\right)$ iterations. The difference can be explained by the restriction enforced by the trust-region constraint on the norm of the steps. Recent work has shown that it is possible to improve the bound for trust-region algorithms using specific definitions of the trust-region radius . The best known iteration bound for a second-order algorithm (that is, an algorithm relying on the use of second-order derivatives and Newton-type steps) is $\mathcal{O}\left(\max\left\{\epsilon_{g}^{-3/2},\epsilon_{H}^{-3}\right\}\right)$. This bound was established originally (under the form of a global convergence rate) in , by considering cubic regularization of Newton’s method. The same result is achieved by the adaptive cubic regularization framework under suitable assumptions on the computed step . Recent proposals have shown that the same bound can be attained by algorithms other than cubic regularization. A modified trust-region method , a variable-norm trust-region scheme , and a quadratic regularization algorithm with cubic descent condition all achieve the same bound.

When $\epsilon_{g}=\epsilon_{H}=\epsilon$ for some $\epsilon\in(0,1)$, all the bounds mentioned above reduce to $\mathcal{O}(\epsilon^{-3})$. It has been established that this order is sharp for the class of second-order methods , and it can be proved for a wide range of algorithms that make use of second-order derivative information; see . Setting $\epsilon_{H}=\epsilon^{1/2}$ and $\epsilon_{g}=\epsilon$ for some $\epsilon>0$ yields bounds varying between $\mathcal{O}(\epsilon^{-3})$ and $\mathcal{O}(\epsilon^{-3/2})$, the latter being again optimal within the class of second-order algorithms .

A new trend in complexity analyses has emerged recently, that focuses on measuring not just the number of iterations to achieve (2) but also the computational cost of the iterations. Two independent proposals, respectively based on adapting accelerated gradient to the nonconvex setting and approximately solving the cubic subproblem , require $\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\epsilon^{-7/4}\right)$ operations (with high probability, showing only dependency on $\epsilon$) to find a point $x_{k}$ that satisfies

with $L_{H}$ being a Lipschitz constant of the Hessian. The difference factor of $\epsilon^{-1/4}$ by comparison with the complexities of the previous paragraph is due to the cost of computing a negative eigenvalue of $\nabla^{2}f(x_{k})$ and/or the cost of solving the linear system. A later proposal focuses on solving cubic subproblems via gradient descent, together with an inexact eigenvalue computation: It satisfies (3) in at most $\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\epsilon^{-2}\right)$ with high probability. Another technique requires only gradient computations, with noise being added to some iterates. It reaches with high probability a point satisfying (3) in at most $\mathcal{O}\left(\log^{4}\left(\frac{1}{\epsilon}\right)\epsilon^{-2}\right)$ iterations. Up to the logarithmic factor, this bound is characteristic of gradient-type methods, but classical work establishes only first-order guarantees . Although this setting is not explicitly addressed in the cited papers, it appears that to reach an iterate satisfying (2) with $\epsilon_{g}=\epsilon_{H}=\epsilon$, the methods studied in would require $\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\epsilon^{-7/2}\right)$ iterations, while the methods described in and could require $\mathcal{O}\left(\log\left(\frac{1}{\epsilon}\right)\epsilon^{-3}\right)$ and $\mathcal{O}\left(\log^{4}\left(\frac{1}{\epsilon}\right)\epsilon^{-3}\right)$ iterations, respectively. Although these bounds look worse than those of classical nonlinear optimization schemes, they are more informative, in that they not only account for the number of outer iterations of the algorithm, but also for the cost of performing each outer iteration (often measured in terms of the number of inner iterations, each of which has similar cost). We note, however, that unlike the classical complexity results, the newer procedures make use of randomization, so the bounds typically hold only with high probability.

Our goal in this paper is to describe an algorithm that achieves optimal complexity, whether measured by the number of iterations required to satisfy the condition (2) or by an estimate of the number of fundamental operations required (gradient evaluations or Hessian-vector multiplications). Each iteration of our algorithm takes the form of a step calculation followed by a backtracking line search. (To our knowledge, ours is the first line-search algorithm that is endowed with a second-order complexity analysis.) The “reference” version of our algorithm is presented in Section 2, along with its complexity analysis. In this version, we assume that two key operations — solution of the linear equations to obtain Newton-like steps and calculation of the most negative eigenvalue of a Hessian — are performed exactly. In Section 3, we refine our study by introducing inexactness into these operations, and adjusting the complexity bounds appropriately. Finally, we discuss the established results and their practical connections in Section 4.

A Line-Search Algorithm Based on Exact Step Computations

We now describe an algorithm based on exact computation of search directions, in particular, the Newton-like search directions and the eigenvector that corresponds to the most negative eigenvalue of the Hessian.

We use a standard line-search framework [18, Chapter 3]. Starting from an initial iterate $x_{0}$, we apply an iterative scheme of the form $x_{k+1}=x_{k}+\alpha_{k}d_{k}$, where $d_{k}$ is a chosen search direction and $\alpha_{k}$ is a step length computed by a backtracking line-search procedure.

Algorithm 1 defines our method. Each iteration begins by evaluating the gradient, together with the curvature of the function along the gradient direction. This information determines whether the negative gradient direction is a suitable choice for search direction $d_{k}$, and if so, what scaling should be applied to it. If not, we compute the minimum eigenvalue of the Hessian. The corresponding eigenvector is used as the search direction whenever the eigenvalue is sufficiently negative. Otherwise, we compute a Newton-like search direction, adding a regularization term if needed to ensure sufficient positive definiteness of the coefficient matrix. There are a total of five possible choices for the search direction $d_{k}$ (including two different scalings of the negative gradient). Table 1 summarizes the various steps that can be performed and the conditions under which those steps are chosen.

Once a search direction has been selected, a backtracking line search is applied with an initial choice of $1$. A sufficient condition related to the cube of the step norm must be satisfied; see (7). Such a condition has been instrumental in the complexity analysis of recently proposed Newton-type methods achieving the best known iteration complexity rates .

At most one eigenvector computation and one linear system solve are needed per iteration of Algorithm 1, along with a gradient evaluation and the Hessian-vector multiplication required to calculate $R_{k}$.

The algorithm contains two tests for termination, with the option of switching to a “Local Phase” instead of terminating at a point that satisfies approximate second-order conditions. The Local Phase aims for rapid local convergence to a point satisfying second-order necessary conditions for a local solution; it is detailed in Algorithm 2. Termination (or switch to the Local Phase) occurs at an iteration $k$ at which an $(\epsilon_{g},\epsilon_{H})$-approximate second-order critical point is reached, according to the following definition:

where $g_{k}=\nabla f(x_{k})$, etc. As we see below, the quantity $\min\left\{\|g_{k}\|,\|g_{k+1}\|\right\}$ arises naturally in the decrease formula we establish for the steps computed by Algorithm 1. In fact, for the methods we reviewed in introduction, one observes that the decrease formulas obtained for their steps either involve only $\|g_{k}\|$ , only $\|g_{k+1}\|$ , or the minimum of the two quantities . The later case appears due to the presence of both gradient-type (see Lemma 2.3) and Newton-type steps (see Lemmas 2.5 and 2.7).

The main convergence results of this section are complexity results on the number of iterations or function evaluations required to satisfy condition (8) for the first time. (Algorithm 2 makes provision for re-entering the main algorithm, if the approximate second-order conditions are violated at any point. This re-entry feature is not covered by our complexity analysis.)

2 Iteration Complexity

We now establish a complexity bound for Algorithm 1, in the form of the maximum number of iterations that may occur before the Termination conditions are satisfied for the first time. To this end, we provide guarantees on the decrease that can be obtained for each of the possible choices of search direction.

In the rest of this paper, we make the following assumptions.

The level set $\mathcal{L}_{f}(x_{0})=\{x|f(x)\leq f(x_{0})\}$ is a compact set.

The function $f$ is twice Lipschitz continuously differentiable on an open neighborhood of $\mathcal{L}_{f}(x_{0})$, and we denote by $L_{g}$ and $L_{H}$ the respective Lipschitz constants for $\nabla f$ and $\nabla^{2}f$ on this set.

We point out that the choice $U_{H}=L_{g}$ is a valid one for theoretical purposes. However, $U_{H}$ will serve as an explicit parameter of our inexact method in Section 3, so we use separate notation, to allow $U_{H}$ to be an overestimate of $L_{g}$.

An immediate consequence of these assumptions is that for any $x$ and $d$ such that Assumption 2 is satisfied at $x$ and $x+d$, we have

The following four technical lemmas derive bounds on the decrease obtained from each type of step. The proofs are rather similar to each other, and follow the usual template for backtracking line-search methods.

We begin with negative curvature directions, showing that our choices for initial scaling yield a decrease proportional to the cube of the (negative) curvature in that direction.

Under Assumption 2, suppose that the search direction for the $k$-th iteration of Algorithm 1 is chosen either as $d_{k}=\frac{R_{k}}{\|g_{k}\|}g_{k}$ with $R_{k}<-\epsilon_{H}$ in Step 1 or $d_{k}=v_{k}$ in Step 2. Then the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$ with $j_{k}\leq j_{e}+1$, where

and the decrease in the function value resulting from the chosen step length satisfies

For the direction $d_{k}=R_{k}g_{k}/\|g_{k}\|$, we have

For the other choice $d_{k}=v_{k}$, we have $d_{k}^{T}\nabla^{2}f(x_{k})d_{k}=\lambda_{k}^{3}=-\|d_{k}\|^{3}$, so that in both cases we have

Thus, if the unit value $\alpha_{k}=1$ is accepted by (7), the result (12) holds trivially.

Suppose now that the unit step length is not accepted. Then the choice $\alpha=\theta^{j}$ does not satisfy the decrease condition (7) for some $j\geq 0$. Using (10) and the definition of $d_{k}$, we obtain

where the last line follows from (13). Therefore, we have

which holds only if $j\leq j_{e}$ by definition of $j_{e}$. Thus, the line search must terminate with (7) being satisfied for some value $j_{k}\leq j_{e}+1$. Because the line search did not stop with step length $\theta^{j_{k}-1}$, we must have

As a result, the decrease satisfied by the step $\alpha_{k}d_{k}=\theta^{j_{k}}d_{k}$ is such that

This inequality, together with the analysis for the case of $\alpha_{k}=1$, establishes the desired result.

The second result concerns use of the step $d_{k}=-g_{k}/\|g_{k}\|^{1/2}$ in the case in which the curvature of the function along the gradient direction is small.

Let Assumptions 1 and 2 hold. Then, if at the $k$-th iteration of Algorithm 1, the search direction is $d_{k}=-g_{k}/\|g_{k}\|^{1/2}$, the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$, with $j_{k}\leq j_{g}+1$, where

and the resulting step length $\alpha_{k}$ is such that

Recall that the choice $d_{k}=-g_{k}/\|g_{k}\|^{1/2}$ is adopted only when $\|g_{k}\|>\epsilon_{g}$ and $|R_{k}|\leq\epsilon_{H}$. If the unit step length $\alpha_{k}=1$ is accepted, we have

satisfying (16). Otherwise, it means that there exists $j\geq 0$ for which the decrease condition (7) is not satisfied using the step size $\theta^{j}$. For such $j$, we have from (10) that

Therefore, at least one of the two terms between brackets must be nonnegative. If

we have $\theta^{j}\geq\frac{5}{3}\|g_{k}\|^{1/2}\epsilon_{H}^{-1}$. On the other hand, if

then $\theta^{j}\geq\sqrt{\frac{1}{L_{H}+\eta}}$. Putting the two bounds together, we have that

Since $j>j_{g}$ contradicts (18c), the line search terminates with (7) being satisfied for some value $j_{k}\leq j_{g}+1$. Since (7) did not hold for $\alpha=\theta^{j_{k}-1}$, we have from (18b) that

The decrease obtained by the step length $\alpha_{k}=\theta^{j_{k}}$ thus satisfies

Thus (16) is also satisfied in the case of $\alpha_{k}<1$, completing the proof.

Lemma 2.3 describes the reduction that can be achieved along the negative gradient direction when the curvature of the function in this direction is modest. When this curvature is significantly positive (or when this curvature is slightly positive but the gradient is small), we compute the minimum Hessian eigenvalue (Step 2) and consider other options for the search direction.

Our next result concerns the decrease that can be guaranteed by the Newton step, when it is computed.

Let Assumptions 1 and 2 hold. Suppose that the Newton direction $d_{k}=d^{n}_{k}$ is used at the $k$-th iteration of Algorithm 1. Then the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$, with $j_{k}\leq j_{n}+1$, where

Note first that the Newton direction $d_{k}=d^{n}_{k}$ is computed only when $\nabla^{2}f(x_{k}))\succ\epsilon_{H}I$, so we have

Suppose first that the step length $\alpha_{k}=1$ satisfies the decrease condition (7). Then from (5) and (10), we have

We thus have the following bound on the decrease obtained with the unitary Newton step:

Suppose now that the unit step length does not allow for a sufficient decrease as measured by (7). Then this condition must fail for $\alpha_{k}=\theta^{j}$ for some $j\geq 0$. For this value, we have from (10) that

where we used $\nabla^{2}f(x_{k})\succeq\epsilon_{H}I$ for the final inequality. This relation holds in particular for $j=0$, in which case it gives

leading to the following lower bound on the norm of the Newton step:

More generally, for any integer $j$ such that the decrease condition is not satisfied, we have from (24) that

For any $j>j_{n}$, the last inequality is violated since

where we used (22) for the final inequality. This proves that the condition (7) will be satisfied by some $j_{k}\leq j_{n}+1$. Since $\alpha=\theta^{j_{k}-1}$ does not fulfill the decrease requirement, it follows from (26) that

By substituting this lower bound into the sufficient decrease condition, and then using (25), we obtain

where the final inequality is from (25). We obtain the required result by combining this inequality with the bound (23) for the case of $\alpha_{k}=1$.

Our last intermediate result addresses the case of a regularized Newton step.

Let Assumptions 1 and 2 hold. Suppose that $d_{k}=d^{r}_{k}$ at the $k$-th iteration of Algorithm 1. Then the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$, with $j_{k}\leq j_{r}+1$, where

Note first that the regularized Newton step is taken only when $\nabla^{2}f(x_{k})\succeq-\epsilon_{H}I$. Thus the minimum eigenvalue of the coefficient matrix in (6) is $\lambda_{k}+2\epsilon_{H}\geq\epsilon_{H}$, and we have

Suppose first that the unit step is accepted. Then the gradient norm at the new point satisfies

By treating the left-hand side as a quadratic in $\|d_{k}\|$, and applying Lemma A.1 with $a=2$, $b=2L_{H}$ and $t=\|\nabla f(x_{k}+d_{k})\|/\epsilon_{H}^{2}$, we obtain from this bound that

Therefore, if the unit step is accepted, we have

If the unit step does not yield a sufficient decrease, there must be a value $j\geq 0$ such that (7) is not satisfied for $\alpha=\theta^{j}$. For such $j$, and using again (10), we have

Thus, for any $j\geq 0$ for which sufficient decrease is not obtained, one has

Meanwhile, we have from the definition of $j_{r}$ that

using the upper bound (29). By comparing this bound with (32), we deduce that the backtracking line-search procedure terminates with $j_{k}\leq j_{r}+1$, where $j_{k}\geq 1$ by our earlier assumption. Thus, since (32) is satisfied for $j=j_{k}-1$, we have

By combining this bound with (31), obtained for the unit-step case, we obtain the result.

By combining the estimates of function decrease proved in the lemmas above, we bound the number of iterations needed by Algorithm 1 to satisfy the approximate second-order optimality conditions (8).

Let Assumptions 1 and 2 hold. Then Algorithm 1 reaches an iterate that satisfies (8) in at most

Suppose $l$ is an iteration at which the conditions for termination are not satisfied. We consider in turn the various types of steps that could have been taken at iteration $l$, and obtain a lower bound on the amount of decrease obtained from each. Table 1 is helpful in working through the various cases. We consider two main cases, and several subcases.

From Table 1, we see that in this case, the search direction is either a scaling of $-g_{k}$, or the most-negative-curvature direction $v_{k}$. When $R_{l}<-\epsilon_{H}$, we have $d_{l}=\frac{R_{l}}{\|g_{l}\|}g_{l}$, and Lemma 2.1 indicates the following bound on function decrease:

When $R_{l}\in[-\epsilon_{H},\epsilon_{H}]$ and $\|g_{l}\|>\epsilon_{g}$, we have $d_{l}=-g_{l}/\|g_{l}\|^{1/2}$. Thus, using Lemma 2.3, we have

For the remaining cases of “$\|g_{l}\|\leq\epsilon_{g}$ and $R_{l}\in[-\epsilon_{H},\epsilon_{H}]$” and “$\|g_{l}\|>\epsilon_{g}$ and $R_{l}>\epsilon_{H}$”, the search direction is necessarily $v_{l}$. We have from Lemma 2.1 that

Case 2: $\lambda_{l}\geq-\epsilon_{H}$, $\|g_{l}\|>\epsilon_{g}$, and $\|g_{l+1}\|>\epsilon_{g}$.

In this case, we have three possible choices for the search direction. The first one is $d_{l}=-g_{l}/\|g_{l}\|^{1/2}$, in which case we have from Lemma 2.3 that

The second possible choice is the Newton direction $d_{l}=d^{n}_{l}$. Using Lemma 2.5, we obtain

The third choice is the regularized Newton direction $d_{l}=d^{r}_{l}$, for which Lemma 2.7 yields

By putting all these bounds together, we obtain the following lower bound on the decrease in $f$ on iteration $l$:

where $c$ is defined in (34). Consequently, summing across all iterations up to $k$ yields

which implies that $k$ is bounded above by (33). Therefore, there must exist a finite index $k_{\epsilon}$ such that (8) is satisfied. For this index, the bound (33) applies, hence the result.

We now look further into the various components of the bound established in Theorem 2.9.

The result (33) makes explicit the variation of the bound with respect to the two tolerances. As this result differs from those in the literature, we follow two usual approaches to ease the comparison with other methods. Letting $\epsilon_{g}=\epsilon$ and $\epsilon_{H}=\sqrt{\epsilon}$ for some $\epsilon\in(0,1)$ allows to equate all components of the maximum term in (33); indeed,

and therefore our bound is $\mathcal{O}(\epsilon^{-3/2})$. On the other hand, the choice $\epsilon_{g}=\epsilon_{H}=\epsilon$, that puts first- and second-order requirement on an equal footing, leads to a bound in $\mathcal{O}(\epsilon^{-3})$. Both match the optimal bounds known for second-order globally convergent methods in terms of iteration count.

Dependencies on problem-algorithmic constants

Although our main goal is to analyze dependencies with respect to the tolerances, our bounds can also reflect dependencies on problem-dependent quantities, namely, the initial function value discrepancy $f(x_{0})-f_{\mbox{\rm\scriptsize low}}$ and the Lipschitz constants $L_{g}$ and $L_{H}$. It can be seen from the lemmas of this subsection that

As a result, the iteration complexity of our method is in

3 Evaluation/Inner Iteration Complexity

We now discuss the function evaluation complexity of Algorithm 1, which counts the number of function calls required by the algorithm before its termination conditions are satisfied. We need to refine the iteration complexity analysis of Section 2.2 to take into account the function evaluations associated with the backtracking line-search process.

Suppose that Assumptions 1 and 2 hold. The number of function evaluations required by Algorithm 1 prior to reaching a point that satisfies (8) is at most

and $\mathcal{C}$ is defined as in Theorem 2.9.

Theorem 2.9 gives a bound on the number of iterations. By Lemmas 2.1–2.7, a bound on the corresponding number of function evaluations is

Using the definitions of $j_{e}$, $j_{g}$, $j_{n}$, and $j_{r}$ from Lemmas 2.1, 2.3, 2.5, and 2.7, respectively, and the fact that $\epsilon_{g},\epsilon_{H}\in(0,1)$ yields the result.

With our specific choices of $\epsilon_{g}$ and $\epsilon_{H}$ mentioned in the previous section, the evaluation complexity bounds are $\mathcal{O}\left(\log(\frac{1}{\epsilon})\epsilon^{-3/2}\right)$ and $\mathcal{O}\left(\log(\frac{1}{\epsilon})\epsilon^{-3}\right)$, respectively. We can also derive a bound that includes dependencies on problem constants; for instance, the bound corresponding to $\epsilon_{g}=\epsilon_{H}=\epsilon$ is

4 Local Convergence

In the previous sections, we have derived global complexity guarantees for Algorithm 1. We now aim to show rapid local convergence for the variant of the algorithm that invokes the Local Phase, Algorithm 2, rather than terminating as soon as the conditions (8) are satisfied. We note that local convergence results like the one we prove here have in the past gone hand-in-hand with global convergence results in smooth nonconvex optimization (see for example ). More recently, several works in the optimization literature have established rapid local convergence alongside global complexity guarantees .

For this section, we will make the following additional assumption.

The sequence of iterates generated by Algorithm 1 in conjunction with Algorithm 2 converges to a local minimizer, that is, a point $x^{*}$ at which $\nabla f(x^{*})=0$ and $\nabla^{2}f(x^{*})\succ 0$.

Under this assumption, the following result is immediate.

Note that the conditions on $k_{0}$ in Lemma 2.13 are such that the combined strategy of Algorithm 1-Algorithm 2 will have entered the Local Phase (Algorithm 2) before iteration $k_{0}$, and will stay in this phase at all subsequent iterations.

We now establish a local quadratic convergence result.

Suppose that Assumptions 1, 2, and 3 are satisfied, and let $\mu$ and $k_{0}$ be as defined in Lemma 2.13. Then for every $k\geq k_{0}$, the method always takes the Newton direction with a unit step length, and we have

Let $k\geq k_{0}$, so that we are in the Local Phase (Algorithm 2) at iteration $k$. By Lemma 2.13, the Hessian at $\nabla^{2}f(x_{k})$ is positive definite, with smallest eigenvalue bounded below by $\mu>0$. Thus Algorithm 2 computes the Newton direction $d_{k}=d^{n}_{k}$, and we have

Thus if the sufficient decrease condition $f(x_{k}+d_{k})-f(x_{k})\leq-\frac{\eta}{6}\|d_{k}\|^{3}$ is not satisfied for the unit step, we must have

which by the bound $\|d_{k}\|\leq\|g_{k}\|/\mu$ can be true only if

which contradicts (38). Thus the unit Newton step is taken, and we have

A Variant with Inexact Directions

In Section 2, we have assumed that certain linear-algebra operations in Algorithm 1 — the linear system solves of (5) and (6) and the eigenvalue / eigenvector computation of (4) — are performed exactly. In a large-scale setting, the cost of these operations can be prohibitive, so iterative techniques that perform these operations inexactly are of interest. In this section, we describe inexact methods for these key operations, and examine their consequences for the complexity analysis.

The problem of finding the minimum eigenvalue of the matrix in (4) and its associated eigenvector can be reformulated as one of finding the maximum eigenvalue and eigenvector of a positive semidefinite matrix. The Lanczos algorithm with a random starting vector is an appealing option for the latter problem, yielding an $\epsilon$-approximate eigenvector in $\mathcal{O}\left(\log(n/\delta)\epsilon^{-1/2}\right)$ iterations, with probability at least $1-\delta$ . This fact has been used in several methods that achieve fast convergence rates . In order to apply this method to a matrix that is not positive definite, one must make use of a bound on the Hessian norm. For sake of completeness, we spell out the procedure in the following lemma.

Let $H$ be a symmetric matrix satisfying $\|H\|\leq M$ for some $M>0$. Suppose that the Lanczos procedure is applied to find the largest eigenvalue of $MI-H$ starting at a random vector uniformly distributed over the unit sphere. Then, for any $\varepsilon>0$ and $\delta\in(0,1)$, there is a probability at least $1-\delta$ that the procedure outputs a unit vector $v$ such that

After at most $n$ iterations, the procedure obtains a unit vector $v$ such that $v^{\top}Hv=\lambda_{\mbox{\rm\scriptsize{min}}}(H)$, with probability $1$.

By definition, the matrix $H^{\prime}=MI-H$ it is a symmetric positive semidefinite matrix with its spectrum lying in $[0,2M]$. Applying the Lanczos procedure to this matrix from a starting point drawn randomly from the unit sphere yields a unit vector $v$ such that

in no more than $\min\left\{n,\frac{\ln(n/\delta^{2})}{4\sqrt{\varepsilon/(2M)}}\right\}$ iterations with probability at least $1-\delta$. (This result is from [15, Theorem 4.2] extended by a continuity argument from the positive definite case to the positive semidefinite case; see [15, Remark 7.5].) Moreover, using (42), we have

Lemma 3.1 admits the following variant, for the case in which we fix the number of Lanczos iterations.

Let $H$ be a symmetric matrix with $\|H\|\leq M$. Suppose that $q$ iterations of the Lanczos procedure are applied to find the largest eigenvalue of $MI-H$ starting at a random vector uniformly distributed over the unit sphere. Then for any $\varepsilon>0$, the procedure outputs a unit vector $v$ such that $v^{\top}Hv\leq\lambda_{\min}(H)+\varepsilon$ with probability at least

We point out that the choice $\delta=0$ (or, equivalently, $q=n$) is possible, that is, after $n$ iterations, the Lanczos procedure started with a random vector uniformly generated over the unit sphere returns an approximate eigenvector with probability one [15, Theorem 4.2 (a)].

2 Inexact Newton and Regularized Newton Directions: Conjugate Gradient Method

Here we describe the use of the conjugate gradient (CG) algorithm to solve the symmetric positive definite linear systems (5) or (6) — the Newton and regularized Newton equations, respectively. The conjugate gradient method is the most popular iterative method for positive definite linear systems, due to its rich convergence theory and strong practical performance. It has also been popular in the context of nonconvex smooth minimization; see . It requires only matrix-vector operations involving the coefficient matrix (often these can be found or approximated without explicit knowledge of the matrix) together with some vector operations. It does not require knowledge or estimation of the extreme eigenvalues of the matrix.

We apply CG to a system $Hd=-g$ where there are positive quantities $m$ and $M$ such that $mI\preceq H\preceq MI$, so that the condition number $\kappa$ of $H$ is bounded above by $M/m$. Standard convergence theory indicates that CG outputs a vector $d$ such that $\|Hy+g\|\leq\zeta\|g\|$ (for $\zeta\in(0,1)$) in

with $\kappa$ being the condition number of $H$ (we obtain the result as a Corollary of Lemma 3.4 below). We use a different stopping criterion, namely

for some $\zeta\in(0,1)$. This criterion is stronger than the one typically used in truncated Newton-Krylov methods, in that we require the residual norm to be bounded by a multiple of the norm of the approximate direction, as well as being bounded by a specified fraction of the initial residual norm. The extra criterion resembles the so-called s-condition arising in cubic regularization techniques, where the approximate minimizer $s_{k}$ of the cubic model $m_{k}$ is required to satisfy

This property provides a lower bound on $\|s_{k}\|$, that is instrumental in obtaining the optimal complexity order of $\mathcal{O}(\epsilon_{g}^{-3/2})$ for first-order convergence . Our condition replaces $\|s_{k}\|^{2}$ by $m\|d_{k}\|$, but serves a similar purpose.

The next lemma establishes a bound on the number of CG iterations needed to reach the desired accuracy.

Let $Hd=-g$ be a linear system with $H$ symmetric and $mI\preceq H\preceq MI$, where $m\in(0,1)$, $M>0$, and $\|g\|>0$. Then the conjugate gradient algorithm computes a vector $d$ such that (44) holds for some $\zeta\in(0,1)$ in at most

Let $d^{(q)}$ be the iterate obtained at the $q$-th iteration of the conjugate gradient method applied to $Hd=-g$, with $d^{(0)}=0$. The classical bound on the behavior of the conjugate gradient residual [18, Section 5.1] yields

where $\|x\|_{H}=\sqrt{x^{\top}Hx}$. From this definition and the bounds on the spectrum of $H$, we have

By substituting these bounds into (47), we obtain the following relation:

Thus, as long as our stopping criterion is not satisfied, we have

Furthermore, defining $r^{(q)}=Hd^{(q)}+g$, we have

for all $q\geq 1$, where we used the fact that using the facts that $r^{(0)}=g$ and that in CG, the residuals are orthogonal: $(r^{(i)})^{T}r^{(j)}=0$ for $i\neq j$. Using this bound within (49), we obtain

By taking logarithms on both sides, we arrive at

where the bound $\ln(1+\frac{1}{t})\geq\frac{1}{t+1/2}$ was used to obtain the last inequality.

3 Complexity Analysis Based on Inexact Computations

We present a variant of our main algorithm, specified as Algorithm 3, in which computation of approximate eigenvectors and linear system solves are performed inexactly by the means described above. Algorithm 3 requires two parameters not used in Algorithm 1: the upper bound $U_{H}$ on the Hessian norms, defined in (9), and a probability threshold $\delta$. As we expect only to recover inexact global complexity guarantees, the method does not exploit a local phase.

When Algorithm 3 terminates, condition (8) must hold. At termination, we have $\min(\|g_{k}\|,\|g_{k+1}\|)\leq\epsilon_{g}$ and $\lambda^{i}_{k}\geq-\tfrac{1}{2}\epsilon_{H}$. With high probability, $\lambda^{i}_{k}$ is within $\tfrac{1}{2}\epsilon_{H}$ of $\lambda_{\mbox{\rm\scriptsize{min}}}(\nabla^{2}f(x_{k}))$, so we must have $\lambda_{\mbox{\rm\scriptsize{min}}}(\nabla^{2}f(x_{k}))\geq-\epsilon_{H}$, thus satisfying (8).

Table 2 shows a summary of the possible choices for the search direction. It shows the same number of cases as Table 1, with the context now determined by the eigenvalue estimate $\lambda_{k}^{i}$, with one exception. There is an extra row for the case $\|g_{k}\|\leq\epsilon_{g}$, $R_{k}\in[-\epsilon_{H},\epsilon_{H}]$, $\lambda^{i}_{k}>\tfrac{3}{2}\epsilon_{H}$, because of possible (but low-probability) failure of the randomized Lanczos process to detect the smallest eigenvalue of $\nabla^{2}f(x_{k})$ to the required accuracy. Table 2 mentions two additional lemmas, that respectively replace Lemmas 2.5 and 2.7 in order to take inexactness into account. We state and prove these results next.

Let Assumptions 1 and 2 hold. Suppose that an inexact Newton direction $d_{k}=d^{in}_{k}$ is computed at the $k$-th iteration of Algorithm 3. Then with probability at least $1-\delta$, the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$, with $j_{k}\leq j_{inr}+1$, where

We observe first that when the Newton step is computed in Algorithm 3, we have from (40) that

with probability $1-\delta$. Suppose first that the step length $\alpha_{k}=1$ satisfies the decrease condition (53). Then, defining

and using the inexactness criterion for the inexact Newton step $d_{k}$, we find that the gradient at the next point $x_{k}+d_{k}$ satisfies

We obtain a lower bound on $\|d_{k}\|$ by taking the root of the above quadratic and applying Lemma A.1 with $a=\zeta\epsilon_{H}/2$, $b=2L_{H}\epsilon_{H}^{2}$, and $t=\|\nabla f(x_{k}+d_{k})\|/\epsilon_{H}^{2}$ to obtain

Therefore, taking the inexact Newton step with a unit step length guarantees

so the inequality (55) is satisfied in the case of a unit step $\alpha_{k}=1$.

To complete the proof, consider the case in which the unit step length does not lead to sufficient decrease. In that case, for any value $j\geq 0$ such that (53) is not satisfied, we have

Thus, for any $j\geq 0$ for which sufficient decrease is not obtained, we have

In particular, since (58) holds for $j=0$, we have

By the definitions of $d_{k}$ and of $r_{k}$ in (56), we also have the following upper bound on its norm:

using again the fact that $g_{k}$ and $r_{k}$ are orthogonal (by the properties of the CG algorithm), as well as the criterion (51) and the bound (9).

As a result, (58) is violated for $j>j_{inr}$, which means that the line search must terminate with a step length $\alpha_{k}=\theta^{j_{k}}$ satisfying (53), with $1\leq j_{k}\leq j_{inr}+1$. Since the index $j=j_{k}-1\geq 0$ satisfies (58), we have

and from the sufficient decrease condition, we have

where the second inequality follows from (60) and the third inequality follows from (59) (using the fact that $\theta\in(0,1)$). Hence, the claim (55) is satisfied in the case of non-unit step length $\alpha_{k}$ too, and the proof is complete.

Let Assumptions 1 and 2 hold. Suppose that an inexact regularized Newton direction $d_{k}=d^{ir}_{k}$ is computed at the $k$-th iteration of Algorithm 3. Then with probability at least $1-\delta$, the backtracking line search terminates with step length $\alpha_{k}=\theta^{j_{k}}$, with $j_{k}\leq j_{inr}+1$, where $j_{inr}$ is defined as in (54), and we have

The inexact regularized Newton step is computed only when $-\tfrac{1}{2}\epsilon_{H}\leq\lambda^{i}_{k}\leq\tfrac{3}{2}\epsilon_{H}$, so from (40) with $\varepsilon=\frac{1}{2}\epsilon_{H}$, we have

with probability at least $1-\delta$. Suppose first that the step length $\alpha_{k}=1$ satisfies the decrease condition (53). Then, defining $r_{k}=(\nabla^{2}f(x_{k})+2\epsilon_{H}I)d_{k}+g_{k}$, we have that

Reasoning as in (57), with $\tfrac{4+\zeta}{2}$ replacing $\tfrac{\zeta}{2}$, we obtain the following lower bound on $\|d_{k}\|$:

Therefore, taking the unit regularized Newton step guarantees

so the result of the theorem holds in the case in which the unit step satisfies the sufficient decrease condition.

To complete the proof, we consider the case in which $\alpha_{k}<1$. In that case, for any value $j\geq 0$ such that (53) is not satisfied, we have from the definition of $r_{k}$, the bound on $\|r_{k}\|$ in the definition of $d^{ir}_{k}$, and (62) that

Thus, for any $j\geq 0$ for which sufficient decrease is not obtained, one has

In particular, setting $j=0$ in this expression, we obtain

The right-hand side of (64) is bounded below, since

where we used again the orthogonality of $g_{k}$ and $r_{k}$ (from the properties of conjugate gradient) as well as the condition (52). For any $j>j_{inr}$, we have

where the last inequality follows from (66). Therefore, (64) is violated for $j>j_{inr}$, which means that the line search must terminate with a step length $\alpha_{k}=\theta^{j_{k}}$, with $1\leq j_{k}\leq j_{inr}+1$. The previous index $j=j_{k}-1$ satisfies (64), so we have

where the final inequality follows from (65), using the fact that $\theta\in(0,1)$. Thus, condition(61) also holds in the case of $\alpha_{k}<1$, and the proof is complete.

Let Assumptions 1 and 2 hold. Then, Algorithm 3 returns a point $x_{k}$ satisfying (2) in at most

with probability at least $1-\hat{K}\delta$. The constants $c_{e}$, $c_{g}$, $c_{in}$, and $c_{ir}$ are defined in Lemmas 2.1, 2.3, 3.6, and 3.8, respectively.

For any iteration $l$ such that $x_{l}$ does not satisfy (8), we must have that either $\min(\|g_{l}\|,\|g_{l+1}\|)>\epsilon_{g}$ or $\lambda_{\mbox{\rm\scriptsize{min}}}(\nabla^{2}f(x_{l}))<-\epsilon_{H}$, where the latter implies that $\lambda^{i}_{l}<-\tfrac{1}{2}\epsilon_{H}$. Thus, similarly to the proof of Theorem 2.9, we can consider the following two cases.

Case 1: $\lambda^{i}_{l}<-\tfrac{1}{2}\epsilon_{H}$.

From Table 2, we see that the same three choices for $d_{l}$ as in the exact version are possible. if $d_{l}=\frac{R_{l}}{\|g_{l}\|}g_{l}$, we have exactly as in Lemma 2.1 that

When $d_{l}=-g_{l}/\|g_{l}\|^{1/2}$, we have from Lemma 2.3 that

The remaining case corresponds to the choice $d_{l}=v^{i}_{l}$. Since

with probability at least $1-\delta$ in that case, one can use the result of Lemma 2.1 to deduce that

again with probability at least $1-\delta$.

Case 2: $\lambda^{i}_{l}\geq-\tfrac{1}{2}\epsilon_{H}$, $\|g_{l}\|>\epsilon_{g}$, and $\|g_{l+1}\|>\epsilon_{g}$.

In this situation, we have three possible choices of search direction $d_{l}$. If $d_{l}=-g_{l}/\|g_{l}\|^{1/2}$, we have again from Lemma 2.3 that (68) holds. If the inexact Newton direction is taken, we obtain by Lemma 3.6 that

Finally, if the search direction is the inexact regularized Newton direction, that is, $d_{l}=d^{ir}_{l}$, we have from Lemma 3.8 that

By putting all these bounds together, as in the proof of Theorem 2.9, we obtain that the number of iterations before reaching a point satisfying (8) is bounded above by $\hat{K}$ defined in the statement of the theorem.

Recalling that for each of these iterations, there is a probability $\delta$ that the randomized Lanczos iteration in (50) will fail, we bound the probability of failure during the course of the algorithm by $\hat{K}\delta$.

Note that if $\delta$ is chosen large enough such that $1-\hat{K}\delta<0$, Theorem 3.10 is not informative. The same remark holds for the corollary below, that makes use of the results from Sections 3.1 and 3.2 to obtain a bound on the total number of Hessian-vector multiplications and gradient evaluations needed by the procedure (assuming that these operations cost roughly the same).

Suppose the assumptions of Theorem 3.10 hold, and let $\delta\in(0,1)$ be given. Then the total number of gradient evaluations and Hessian-vector multiplications required by Algorithm 3 to reach an iterate satisfying (8) is satisfied is bounded by

The proof follows directly from Lemmas 3.1 and 3.4, setting $M=U_{H}+2$ and $\varepsilon=\epsilon_{H}/2$, noting that for both Newton and regularized Newton steps, the condition number of the respective coefficient matrices can be bounded by $(U_{H}+2)/\epsilon_{H}$.

As in Section 2.2, we can particularize this result to a specific choice of tolerances.

Suppose that the assumptions of Theorem 3.10 hold, and let $\delta\in(0,1)$ be given. Define $\epsilon_{g}=\epsilon$ and $\epsilon_{H}=\sqrt{\epsilon}$, for some $\epsilon\in(0,1)$. Then the number of gradient evaluations and Hessian-vector products needed to Algorithm 3 to satisfy (8) is bounded by

where $\hat{C}$ is defined as in Theorem 3.10, with probability at least $1-\hat{C}\epsilon^{-3/2}\delta$.

This result is meaningful when $\delta\ll\epsilon^{3/2}$. In terms of the complexity bound, such a choice is not prohibitively small, because $\delta$ enters into the bound (70) only inside a $\log$ term.

We can obtain a bound for the case of $\delta=0$ (that is, almost certainty), at the cost of taking $n$ Lanczos iterations whenever the smallest eigenvalue is needed (see Lemma 3.1). In this case, the bound (70) either becomes $\mathcal{O}\left(\left(n+\ln\left(\epsilon^{-1}\right)\right)\epsilon^{-7/4}\right)$ or $\mathcal{O}(n\epsilon^{-3/2})$, depending on which term dominates in the quantity corresponding to conjugate gradient iterations.

For very large $n$ and $\delta>0$, we can consider that the term involving $\epsilon$ is smaller than $n$ in both minimum expressions in Corollary 3.14. In this case, the bound is

This complexity matches other recent findings .

In terms of dependencies with respect to problem constants, we can reproduce the analysis from Section 2.2, replacing $c_{n}$ and $c_{r}$ by $c_{in}=\mathcal{O}(L_{H}^{-3})$ and $c_{ir}=\mathcal{O}(L_{H}^{-3})$, respectively. For instance, the bound from the previous paragraph is in

We point out that the dependency on $L_{H}$ of our bound is worse than those of , due to the lack of explicit use of this constant within our algorithm. Still, we believe our dependency to match that of other Newton-type methods (although those are not enlightened in the related literature), and we consider such schemes as being more amenable to highly nonlinear settings where estimating such a constant would likely be impractical.

As a final note, we observe that one could also include the number of line-search iterations into our complexity bound. However, this cost is essentially logarithmic in $1/{\epsilon}$, therefore it is dominated by the cost of the linear algebra techniques.

Discussion

Among the many algorithmic frameworks that have been proposed for smooth nonconvex optimization with second-order complexity guarantees, it can be difficult to determine the algorithmic features that affect the complexity analysis, or to understand how the guarantees provided by different algorithms relate to one another. We have presented a second-order complexity analysis of a framework that is based exclusively on line searches along certain directions. It does not require solution of cubic-regularized or trust-region subproblems, or minimization of convexified functions — operations that are needed by other approaches. Our search directions are of several types — gradient, negative-curvature, Newton, and regularized Newton — and we presented a variant of our method that allows inexact direction computation using iterative methods. We believe that ours is the first approach of line-search type to achieve known optimal complexity, among methods that identify points that satisfy approximate second-order necessary conditions.

In addition to the results of this paper, we observe that it is possible to modify our algorithms to attain points that satisfy termination conditions of the form (2) (rather than (8)) by continuing to iterate in the situation in which $\|g_{k+1}\|\leq\epsilon_{g}$ but $\lambda_{min}(\nabla^{2}f(x_{k+1}))<-\epsilon_{H}$. Step $k+1$ then yields a decrease that is a multiple of $\epsilon_{H}^{3}$ (per Lemma 2.1), so the overall complexity estimates are preserved, even if step $k$ in this situation fails to produce a significant decrease in $f$.

In designing the framework of Algorithms 1 and 3, we have made some choices to give preference to one direction choice over another, and we have also incorporated several types of steps. Given the recent literature in this area, our proposed scheme is actually one particular instance of a broader class of methods with similar complexity guarantees but possibly diverse practical performance. An implementation of our approach would raise several delicate issues, for example, issues associated with failure of the randomized Lanczos procedure for obtaining an estimate of the smallest eigenvalue. An incorrect estimate here could lead to the conjugate gradient method subsequently being applied to an indefinite matrix; a robust implementation would need to detect and recover from such an occurrence. Additionally, the choice of suitable values for the bound on the Hessian norm is likely to be of critical importance. Addressing these concerns in the aim of developing a practical algorithm with good complexity guarantees is the subject of ongoing research.

Acknowledgments

We are grateful to the anonymous referees and associate editor of the original version of the paper, whose construtive comments led to numerous improvements.

Appendix A Technical Result

We prove a technical result that is used in several proofs, including that of Lemma 2.7.

For positive scalars $a$ and $b$, and $t\geq 0$, we have

so the result holds in this case. For $t\in(0,1)$, we need to show that

This claim follows from the following chain of equivalences:

References