Nonlinear conjugate gradient for smooth convex functions^††thanks: Supported in part by a grant from the Natural Sciences and Engineering Research Council (NSERC) of Canada.

The problem under consideration is minimizing an unconstrained $L$-smooth convex function $f:\mathbb{R}^{n}\rightarrow\mathbb{R}$. Recall that a convex function is $L$-smooth if it is differentiable and for all ${\bm{x}},\bm{y}\in\mathbb{R}^{n}$,

Note that convexity implies that the quantity on the left-hand side of this inequality is nonnegative. We will also sometimes consider $\ell$-strongly convex functions, which are defined to be those functions satisfying

for all ${\bm{x}},\bm{y}\in\mathbb{R}^{n}$ for some modulus $\ell\geq 0$. (This definition presumes differentiability. The class of strongly convex functions includes nondifferentiable functions as well, but those functions are not smooth and are therefore not relevant to this paper.) Note that $\ell=0$ means simply that $f$ is differentiable and convex.

Nonlinear conjugate gradient (NCG) [21] is perhaps the most widely used first-order method for unconstrained optimization. It is an extension of linear conjugate gradient (LCG) introduced by Hestenes and Stiefel [13] to minimize convex quadratic functions. The smoothness modulus $L$ for a quadratic function $f({\bm{x}})={\bm{x}}^{T}A{\bm{x}}/2-\bm{b}^{T}{\bm{x}}$, where $A\in\mathbb{R}^{n\times n}$ is positive semidefinite, is the maximum eigenvalue of $A$. LCG is optimal for quadratic functions in a strong sense that $f({\bm{x}}_{k})$, where $k$ is the $k$th iterate, is minimized over all choices of ${\bm{x}}_{k}$ lying in the $k$th Krylov space leading to the estimate of $k=O(\ln(1/\epsilon)\sqrt{L/\ell})$ iterations to reduce the residual of $f$ by a factor $\epsilon$.

Stated precisely, Daniel [9, Theorem 1.2.2] proved the following bound for the $k$th iteration of LCG:

Here and in the remainder of the paper, $f^{*}$ stands for $\min_{{\bm{x}}\in\mathbb{R}^{n}}f({\bm{x}})$ and is assumed to be attained. Therefore, if we define the residual reduction

then a rearrangement and application of simple bounds in (3) shows that $\epsilon_{1}$ reduction is achieved after $k=\ln(4/\epsilon_{1})\cdot\sqrt{L/\ell}/2$ iterations.

However, if the class of problems is enlarged to include all smooth, strongly convex functions, then NCG fails to satisfy this complexity bound and could perform even worse than steepest descent as shown by Nemirovsky and Yudin [18]. Nonetheless, in practice, nonlinear conjugate gradient is difficult to beat even by algorithms with better theoretical guarantees, as observed in our experiments and by others, e.g., Carmon et al. [4].

Several first-order methods have been proposed for smooth strongly convex functions that achieve the running time $O(\ln(1/\epsilon)\sqrt{L/\ell})$ including Nemirovsky and Yudin’s [18] method, which requires solution of a two-dimensional subproblem on each iteration, Nesterov’s accelerated gradient (AG) [19], and Geometric Descent (GD) [2]. In the case of GD, the bound is on iterations rather than gradient evaluations since each iteration requires additional gradient evaluations for the line search. It is known that $\Omega(\ln(1/\epsilon)\sqrt{L/\ell})$ gradient evaluations are required to reduce the objective function residual by a factor of $\epsilon$; see, e.g., the book of Nesterov [20] for a discussion of the lower bound result. Since the algorithms mentioned in this paragraph match the lower bound, these algorithms are optimal up to a constant factor. In the case of smooth convex (not strongly convex) functions, AG requires $O(\epsilon^{-1/2})$ iterations.

The precise bound for accelerated gradient for strongly convex functions is stated below in (14), and this bound can be rearranged and estimated as follows. Define the residual reduction

This reduction can be attained after no more than $k=\ln(1/\epsilon_{2})\cdot\sqrt{L/\ell}$ iterations.

A technical issue in comparing the bounds is the different definitions of $\epsilon_{1}$ and $\epsilon_{2}$. A standard result of extremal eigenvalues states that

Therefore, it follows that $\epsilon_{2}\leq 2\epsilon_{1}\leq(L/\ell)\epsilon_{2}.$ Thus, $\epsilon_{1}$ tightly controls $\epsilon_{2}$, whereas $\epsilon_{2}$ only loosely controls $\epsilon_{1}$, and therefore, the LCG bound is slightly better than the corresponding AG bound in this technical sense.

In practice, conjugate gradient for quadratic functions is typically superior to AG applied to quadratic functions. The reason appears to be not so much the technical reason mentioned in the previous paragraph but the fact that LCG is optimal in a stronger sense than AG: the iterate ${\bm{x}}_{k}$ is the true minimizer of $f({\bm{x}})$ over the $k$th Krylov space, which is defined to be $\{\bm{b},A\bm{b},\ldots,A^{k-1}\bm{b}\}$ assuming ${\bm{x}}_{0}=\bm{0}$. For quadratic problems in which the eigenvalues lie in a small number of clusters, one expects convergence much faster than $O(\ln(1/\epsilon)\sqrt{L/\ell})$. But even in the case of spread-out eigenvalues, conjugate gradient often bests accelerated gradient. Table 1 shows the results for these algorithms for three $1000\times 1000$ matrices:

two of which have clustered eigenvalues, and the third has spread-out eigenvalues.

Note that the three conjugate gradient methods require roughly the same number of iterations, all much less than AG. The function evaluations are: 1 per LCG iteration (i.e., a single matrix-vector multiplication), 2 per C+AG iteration, approximately 2 per NCG iteration, and 1 per AG iteration. The requirement for more than one function evaluation per NCG iteration stems from the need to select a step-size. We return to this point below. Table 1 shows results for the adaptive versions of C+AG and AG for the case that $L$ is not known a priori, discussed in Section 6. Thus, further function-gradient evaluations are needed to estimate $L$, and the number of function-gradient evaluations is increased accordingly. The reader will note from Table 1 that C+AG required only $2$ iterations for $A_{2}$ but $27$ function evaluations due to the overhead of initially estimating $L$, and a similar pattern is noted for $A_{3}$. The adaptive estimation doubles the function-gradient evaluations for AG from 1 to 2 per iteration.

In the previous paragraph and for the remainder of this paper, we speak of a “function-gradient evaluation”, which means, for a given ${\bm{x}}$, the simultaneous evaluation of $f({\bm{x}})$ and $\nabla f({\bm{x}})$. We regard function-gradient evaluation as the atom of work. This metric may be inaccurate in the setting that $f({\bm{x}})$ can be evaluated much more efficiently than $\nabla f({\bm{x}})$ for an algorithm that separately uses $f$ and $\nabla f$. However, in most applications evaluating $f,\nabla f$ together is not much more expensive than evaluating either separately. Theoretically, gradient evaluation is never more than a constant factor larger than function evaluation due to the existence of reverse-mode automatic differentiation; see, e.g., [21].

It is also known that NCG is fast for functions that are nearly quadratic. In particular, a classic result of Cohen [5] shows that NCG exhibits $n$-step quadratic convergence near the minimizer provided that the objective function is nearly quadratic and strongly convex in this neighborhood.

Thus, there is seemingly a gap in the menu of available algorithms because NCG, though optimal for quadratic functions and fast for nearly quadratic functions, is not able to achieve reasonable complexity bounds for the larger class of smooth, strongly convex functions. On the other hand, the methods with the good complexity bounds for this class are all suboptimal for quadratic functions.

We propose C+AG to close this gap. The method is based on the following ideas:

1. 1.

   A line-search for NCG described in Section 2 is used that is exact for quadratic functions but could be poor for general functions. Failure of the line-search is caught by the progress measure described in the next item.

2. 2.

   A progress measure described in Section 3 for nonlinear conjugate gradient based on Nesterov’s estimate sequence is checked. The progress measure is inexpensive to evaluate on every iteration.

3. 3.

   Restarting in the traditional sense of NCG (i.e., an occasional steepest descent step) is used, as described in Section 4

4. 4.

   On iterations when insufficient progress is made, C+AG switches to AG to guarantee reduction in the progress measure as explained in Section 5.

5. 5.

   C+AG switches from AG back to NCG after a test measuring nearness to a quadratic is satisfied. When CG resumes, the progress measure may be modified to a more elaborate procedure discussed in Appendix A, although this elaborate progress measure is disabled by default.

When applied to a quadratic function, the progress measure is always satisfied (this is proved in Section 3), and hence all steps will be NCG (and therefore LCG) steps. The fact that the C+AG satisfies the $O(\ln(1/\epsilon)\sqrt{L/\ell})$ running-time bound of AG for general convex functions follows because of the progress measure. The conjectured final theoretical property of the method, namely, $n$-step quadratic convergence, is described in the Appendix B. Computational testing in Section 8 indicates that the running time of the method measured in function-gradient evaluations is roughly equal to whichever of AG or NCG is better for the problem at hand, and in some cases it outperforms both of them.

Both AG and C+AG use prior knowledge of $L,\ell$. However, $L$ can be estimated by both methods, and $\ell$ may be taken to be 0 if no further information is available. Thus, both algorithms are useable, albeit slower, in the absence of prior knowledge of these parameters. We return to this topic in Section 6.

We conclude this introductory section with a few remarks about our previous unpublished manuscript [16]. That work explored connections between AG, GD, and LCG, and proposed a hybrid of NCG and GD. However, that work did not address the number of function-gradient evaluations because the the evaluation count for the line search was not analyzed. Closer to what is proposed here is a hybrid algorithm from the first author’s PhD thesis [15]. Compared to the PhD thesis, the proposed algorithm herein has several improvements including the line search and the method for switching from AG to CG as well as additional analysis.

The full C+AG algorithm is presented in Section 7. In this section, we present the NCG portion of C+AG, which is identical to most other NCG routines and is as follows

For $\beta_{k+1}$ we use the Hager-Zhang formula [12], which is as follows:

For $\alpha_{k}$, we use the following formulas.

The rationale for this choice of $\alpha_{k}$ is as follows. In the case that $f({\bm{x}})={\bm{x}}^{T}A{\bm{x}}/2-\bm{b}^{T}{\bm{x}}$ (quadratic), $\nabla f({\bm{x}}_{k})=A{\bm{x}}_{k}-\bm{b}$ and $\nabla f(\tilde{\bm{x}})=A\tilde{\bm{x}}-\bm{b}=A{\bm{x}}_{k}-\bm{b}+A\bm{p}_{k}/L$, and therefore $\bm{s}=A\bm{p}_{k}$ and $\bm{p}_{k}^{T}\bm{s}=\bm{p}_{k}^{T}A\bm{p}_{k}$. In this case, $\alpha_{k}$ is the exact minimizer of the univariate quadratic function $\alpha\mapsto f({\bm{x}}_{k}+\alpha\bm{p}_{k})$.

In the case of a nonquadratic function, this choice of $\alpha_{k}$ could be inaccurate, which could lead to a poor choice for ${\bm{x}}_{k+1}$. The progress measure in the next section will catch this failure and select a new search direction $\bm{p}_{k}$ in this case.

As for the choice of $\beta_{k}$, we mention that Dai and Yuan [7] carried out an extensive convergence study of different formulas for $\beta_{k}$ (not including the Hager-Zhang formula, but covering many others), and in the end, proposed a certain hybrid formula for $\beta_{k}$. Note that Dai and Kou [8] propose another family with the interesting property [8, Theorem 4.2] that convergence is guaranteed for smooth, strongly convex functions, although no rate is derived.

Our experiments (not reported here) suggest that the choice of $\beta_{k}$ is not critical to the performance of C+AG. The Dai-Yuan paper points out that certain choices of $\beta_{k}$ will lead to nonconvergence if the line-search fails to satisfy the strong Wolfe conditions. However, our line-search does not satisfy even the weak Wolfe conditions. In other words, the criteria used in the previous literature to distinguish among the choices of formulas for $\beta_{k}$ may not be applicable to C+AG.

The progress measure is taken from Nesterov’s [20] analysis of the AG method. The basis of the progress measure is a sequence of strongly convex quadratic functions $\phi_{0}({\bm{x}}),\phi_{1}({\bm{x}}),\ldots$ called an “estimate sequence” in [20].

The sequence is initialized with

where $\bm{v}_{0}:={\bm{x}}_{0}$. Then $\phi_{k}$ for $k\geq 1$ is defined in terms of its predecessor via:

Here, $\theta_{k}\in(0,1)$ and $\bar{\bm{x}}_{k}\in\mathbb{R}^{n}$ are detailed below. Note that the Hessian of $\phi_{0}$ and of the square-bracketed quadratic function in (8) are both multiples of $I$, and therefore inductively the Hessian of $\phi_{k}$ is a multiple of $I$ for all $k$. In other words, for each $k$ there exists $\gamma_{k},\bm{v}_{k},\phi^{*}_{k}$ such that $\phi_{k}({\bm{x}})=\phi_{k}^{*}+(\gamma_{k}/2)\|{\bm{x}}-\bm{v}_{k}\|^{2}$, where $\gamma_{k}>0$. The formulas for $\gamma_{k+1}$, $\bm{v}_{k+1}$, $\phi_{k+1}^{*}$ are straightforward to obtain from (8) given $\gamma_{k}$, $\bm{v}_{k}$, $\phi_{k}^{*}$ as well as $\theta_{k}$ and $\bar{\bm{x}}_{k}$. These formulas appear in [20, p. 73] and are repeated here for the sake of completeness:

We call the sequence $\bar{\bm{x}}_{0},\bar{\bm{x}}_{1},\ldots$ the gradient sequence because the gradient is evaluated at these points.

To complete the formulation of $\phi_{k}(\cdot)$, we now specify $\theta_{k}$ and $\bar{\bm{x}}_{k}$. Scalar $\theta_{k}$ is selected as the positive root of the quadratic equation

The fact that $\theta_{k}\in(0,1)$ is easily seen by observing the sign-change in this quadratic over $[0,1]$. A consequence of (12) and (9) is the identity

Finally, C+AG has three cases for selecting gradient sequence $\bar{\bm{x}}_{k}$. The first case is $\bar{\bm{x}}_{k}:={\bm{x}}_{k}$. The second and third cases are presented later.

The principal theorem assuring progress is Theorem 2.2.2 of Nesterov [20], which we restate here for the sake of completeness. The validity of this theorem does not depend on how the gradient sequence $\bar{\bm{x}}_{k}$ is selected in the definition of $\phi_{k}(\cdot)$.

The sufficient-progress test used for CG is the condition of the theorem, namely, $f({\bm{x}}_{k})\leq\phi_{k}^{*}$. That is, in the case of CG, the witness sequence is also the main sequence computed by CG. As long as this test holds, the algorithm continues to take CG steps, and the right-hand side of (14) assures us that the optimal complexity bound is achieved (up to constants) in the case of both smooth convex functions and smooth, strongly convex functions.

The main result of this section is the following theorem, which states that if the objective function is quadratic, then the progress measure is always satisfied by conjugate gradient.

Line (18) follows because, according to the recursive formula (10), $\bm{v}_{k}$ lies in the affine space ${\bm{x}}_{0}+\mathop{\rm span}\{\nabla f({\bm{x}}_{0}),\cdots,\nabla f({\bm{x}}_{k-1})\}$. Similarly, ${\bm{x}}_{k}$ lies in this space. Therefore, $\bm{v}_{k}-{\bm{x}}_{k}$ lies the $k$th Krylov subspace $\mathop{\rm span}\{\nabla f({\bm{x}}_{0}),\cdots,\nabla f({\bm{x}}_{k-1})\}$. It is well known (see, e.g., [10]) that $\nabla f({\bm{x}}_{k})$ is orthogonal to every vector in this space.

Line (19) is a standard inequality for $L$-smooth convex functions; see, e.g., p. 57 of [20]. Finally (20) follows because ${\bm{x}}_{k}-\nabla f({\bm{x}}_{k})/L$ lies in the $k+1$-dimensional affine space ${\bm{x}}_{0}+\mathop{\rm span}\{\nabla f({\bm{x}}_{0}),\ldots,\nabla f({\bm{x}}_{k})\}$ while ${\bm{x}}_{k+1}$ minimizes $f$ over the same affine space, another well known property of CG. $\square$

We note that this proof yields another proof of Daniel’s result that conjugate gradient converges at a linear rate with a factor $1-\mbox{const}\sqrt{\ell/L}$ per iteration. Specifically, the two bounds are

We see that the reduction factor is better in Daniel’s bound because of the exponent $2k$ rather than $k$, which means approximately reduction of $1-2\sqrt{\ell/L}$ per iteration (assuming $\ell\ll L$) instead of $1-\sqrt{\ell/L}$. Furthermore, the multiplicative constant in Daniel’s bound is better as noted earlier.

A limitation of this proof is that it requires $\bm{v}_{0}={\bm{x}}_{0}$. In the case that the CG iterations follow a sequence of AG iterations, this equality will not hold. We return to this point in Appendix A.

The classical method for coping with failure to make progress in nonlinear CG is ‘restarting’, which means taking a step using the steepest descent direction instead of the conjugate gradient direction. Our proposed C+AG method, when it detects insufficient progress on the $k$th iteration (i.e., the inequality $f({\bm{x}}_{k})\leq\phi_{k}^{*}$ fails to hold), attempts two fall-back solutions. The first fall-back is restarting CG in the traditional manner (steepest descent). It should be noted that this restart does not play a role in the theoretical analysis of the algorithm, but we found that in practice it improves performance (versus immediately proceeding to the second fall-back of the AG method).

The issue of when to restart has received significant attention in the previous literature. One guideline is to restart every $n+1$ iterations. In more modern codes, restarting is done every $kn+1$ iterations, where $k\geq 1$ is a constant, e.g., Hager-Zhang’s CG-Descent restarts every $6n+1$ iterations. We also have implemented this restart in our code, although it is usually not activated because some other test fails first.

Three other conditions for restarting from the previous literature, collectively termed the “Beale-Powell” conditions by Dai and Yuan [6], are

where $c_{1}\in(0,\infty)$, $c_{2}\in[0,1]$, and $c_{3}\in(1,\infty)$. It should be noted that in exact arithmetic on a quadratic objective function, none of these tests will ever be activated. It should also be noted that Powell [24] as well as Dai and Yuan consider other restart directions (which we have not implemented) besides the steepest descent direction.

We have found that including these tests (BP1)–(BP3) caused the code to take more iterations in the majority of tests (although the difference in iteration counts was never large), so we have turned them off by default. In other words, the default values are $c_{1}=\infty$, $c_{2}=0$ (i.e., we restart if $\bm{p}_{k}$ is not a descent direction), and $c_{3}=\infty$. Presumably this is because there are cases when progress is possible according to our progress measure even though a BP condition fails.

One additional restart condition in our code is that if the denominator of (6) is nonpositive, then we restart because this means that calculation of $\alpha_{k}$ has failed.

A recent paper by Buhler et al. [3] considers nonlinear conjugate gradient for problems in machine learning and in particular the role of (BP) restarts. The motivation for their work is somewhat different from ours: they consider general objective functions, not necessarily smooth convex functions, and therefore accelerated gradient is not available for their class of problems. Their first main concern, costly line search, is not an issue for C+AG since its line search requires only one additional gradient per iteration. Their second main concern, nonconjugacy in search directions, is also less of a concern for C+AG since conjugacy is not an ingredient of the C+AG progress measure. They conclude that (BP1)–(BP3) (and (BP1) in particular) are important for a successful nonlinear CG implementation. The fact that their conclusions differ from ours arises from the different aims of the papers.

If the inequality $f({\bm{x}}_{k})\leq\phi_{k}^{*}$ fails to hold also for the steepest descent direction described in the previous section, then the second fall-back is to switch to a sequence of AG steps.

AG steps continue until the following “almost-quadratic” test is satisfied. Let

The termination test for AG is

The rationale for this formula is as follows. One checks with a few lines of algebra that if $f$ were quadratic, then the identity $f({\bm{x}}_{k+1})=f(\bar{\bm{x}}_{k})+q_{k}$ holds, assuming that ${\bm{x}}_{k+1}$ is defined as $\bar{\bm{x}}_{k}-\nabla f(\bar{\bm{x}}_{k})/L$ as in AG (see l. 46 of Algorithm C+AG, Part II below ). Thus, (24) tests whether $f$ in the neighborhood of $\bar{\bm{x}}_{k}$ is behaving approximately as a quadratic function would behave, which indicates that conjugate gradient might succeed.

Note that (24) involves an extra function-gradient evaluation since an ordinary AG iteration evaluates $f,\nabla f$ at $\bar{\bm{x}}_{k}$ but not at ${\bm{x}}_{k+1}$. Therefore, (24) is checked only every 8 iterations to reduce the overhead of the test.

Overall, the C+AG algorithm requires two function-gradient evaluations for a CG step, one to obtain the gradient at ${\bm{x}}_{k}$ and the second to evaluate $\nabla f(\tilde{\bm{x}})$ in (5). Similarly, the first fall-back requires two function-gradient evaluations. Finally, the accelerated gradient step requires another gradient evaluation. Therefore, the maximum number of function-gradient evaluations on a single iteration before progress is made is five. Therefore, Theorem 14 applies to our algorithm except for a constant factor of 5. In practice, on almost every step, either the CG step succeeds or AG is used, meaning the actual number of function-gradient evaluations per iteration lies between 1 and 2.

The C+AG method as described thus far requires knowledge of $L$ and $\ell$, the moduli of smoothness and strong convexity respectively. In case these parameters are not known, the method estimates $L$ and simply assumes $\ell=0$. For many problems (including our examples) $\ell\ll L$, meaning that convergence is governed by the second parenthesized of factor of (14). Thus, taking $\ell=0$ does not appear to significantly impede convergence.

As for estimating $L$, we use the following procedure outlined in §10.4.2 of Beck [1]. Starting from an initial guess of $L=1$, we first decrease $L$ by multiplying by $1/\sqrt{2}$ or increase by multiplying by $\sqrt{2}$ until the condition

is satisfied. This means that C+AG and AG require $O(|\log(L)|)$ iterations for the initial estimate. On subsequent iterations, we re-estimate $L$ on the first iteration of consecutive sequences of CG iterations and on every iteration of AG. On iterations after the first, $L$ cannot decrease; it will either stay the same or increase. Beck proves that using this procedure does not worsen the iteration-complexity of AG when compared to the case that $L$ is known a priori. However, this procedure increases the number of function-gradient evaluations per AG step from one to approximately two. In the case of the conjugate gradient steps, the impact is negligible since only the first in a sequence requires the extra function-gradient evaluation.

It is not necessary to re-estimate $L$ on every CG iteration because the CG step uses $L$ only as the finite-difference step size in its procedure to determine $\alpha_{k}$ as in (4)–(6). This dependence is mild in the sense that, assuming $f$ is quadratic, (4)–(6) will compute the same (correct) value of $\alpha_{k}$ regardless of the step size. All computational tests in Tables 1 and 2 use the adaptive version of C+AG described in this section.

There are procedures in the literature for estimating $\ell$ also, for example, the procedure by Carmon et al. [4]. However, Carmon et al.’s procedure is costly as it involves several trial steps and a guess-and-update procedure for $\ell$. Even in the pure quadratic case, estimating $\ell$ apparently requires multiple steps of the Lanczos method [10] and would be costly. We have not implemented it since, as mentioned above, taking $\ell=0$ appears to be satisfactory.

In this section we present pseudocode for the C+AG algorithm. The pseudocode includes the procedure for estimating $L$. It does not, however, include the elaborate procedure for finding witnesses after CG is restarted as described in Appendix A. We did not include this in the pseudocode because it is disabled by default as discussed in the appendix. See the concluding remarks of the paper for a downloadable implementation in Matlab, which does include that procedure as an option. The input arguments are $f$ and $\nabla f$, the objective function and gradient, ${\bm{x}}_{0}$, the initial guess, gtol, the termination criterion, and $L$, $\ell$, the smoothness modulus and the modulus of strong convexity respectively. Note that $\ell=0$ is a valid argument, and $L=\mathrm{NaN}$ (not a number) is a signal to the algorithm to estimate $L$.