Nonlinear conjugate gradient methods: worst-case convergence rates via computer-assisted analyses

Introduction

We consider the standard unconstrained convex minimization problem

where $f$ is $L$-smooth (i.e., it has an $L$-Lipschitz gradient) and $\mu$-strongly convex. We study the worst-case performances of a few famous variants of nonlinear conjugate gradient methods (NCGMs) for solving (1). More specifically, we study Polak-Ribière-Polyak (PRP) and Fletcher-Reeves (FR) schemes with exact line search. With exact line search, many other NCGMs such as the Hestenes and Stiefel method , the conjugate descent method due to Fletcher , and the Dai and Yuan method reduce to either PRP or FR. Under exact line search, PRP and FR can be presented in the following compact form:

where PRP and FR are respectively obtained by setting $\eta=1$ and $\eta=0$. NCGMs have a long history (see, e.g., the survey and monograph ), but are much less studied compared to their many first-order competitors. For instance, even though FR is generally considered the first NCGM [7, §1], we are not aware of non-asymptotic convergence results for it. On a similar note, some variants of NCGMs are known for their generally good empirical behaviors (which we illustrate in Figure 1) with little of them being backed-up by classical complexity analyses. In this work, we apply the performance estimation approach to ($\mathcal{M}$) for filling this gap by explicitly computing some worst-case convergence properties of PRP and FR with exact line search. This work focuses on exact line search, as it is arguably the most logical starting point to understand the non-asymptotic convergence behavior of NCGMs. In certain cases, the minimizer associated with exact line search has an analytical form, while in others, it can be computed efficiently [11, $\mathsection$9.7.1]. However, in many practical implementations, inexact line searches are employed that try to either approximately minimize $f(x_{k}-\gamma\,d_{k})$ or even just reduce $f$ enough along the ray $x_{k}-\gamma d_{k}$. These inexact methods can be either monotone, which ensures a decrease in $f$ but converges slowly, or nonmonotone, which may allow faster convergence but risks nonrobust tuning [8, $\mathsection$1.2]. Examples of notable monotone inexact line search schemes include backtracking , Goldstein , and Wolfe line searches , and their variants . Nonmonotone schemes include and many others; see [8, pp. 10-14] for a brief review. Despite the computational differences between the two types of line searches, both aim to emulate the exact line search method. Consequently, when using an inexact line search process, any convergence guarantees—defined in terms of iteration numbers—are likely to be worse compared to the exact line search (neglecting the cost of performing exact line search).

The contribution of this paper is twofold. First, we compute worst-case convergence bounds and counter-examples for PRP and FR. These bounds are obtained by formulating the problems of computing worst-case scenarios as nonconvex quadratically constrained quadratic optimization problems (QCQPs), and then by solving them to global optimality. Second, these computations enable us to construct mathematical proofs that establish an improved non-asymptotic convergence bound for PRP, and, to the best of our knowledge, the first non-asymptotic convergence bound for FR. Furthermore, the worst-case bounds for PRP and FR obtained numerically reveal that there are simple adversarial examples on which these methods do not perform better than gradient descent with exact line search (GDEL), leaving very little room for improvements on this class of problems. Since we demonstrate that the convergence results of NCGMs associated with exact line search are already disappointing, we conclude that inexact line searches, which approximate exact line search, are unlikely to offer improvement.

From a methodological point of view, our approach of computing worst-case scenarios and bounds through optimization is part of what is often referred to as performance estimation. While these problems are usually amenable to convex semidefinite programs , this is generally not the case for adaptive first-order methods such as PRP and FR . To study these methods, we evaluate the worst-case performances of ($\mathcal{M}$) by solving nonconvex QCQPs, extending the standard SDP-based approach from developed for non-adaptive methods. This contribution aligns with the spirit of , developed for devising optimal (but non-adaptive) first-order methods.

The paper is organized as follows. In Section 2, we establish non-asymptotic convergence rates for PRP and FR by viewing the search direction $d_{k}$ in ($\mathcal{M}$) as an approximate gradient direction. In Section 3, we compute the exact numerical values of the worst-case $\nicefrac{{f(x_{N})-f_{\star}}}{{f(x_{0})-f_{\star}}}$ and $\nicefrac{{f(x_{k+N})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ for PRP and FR by formulating the problems as nonconvex QCQPs and then solving them to certifiable global optimality using a custom spatial branch-and-bound algorithm. The solutions to these QCQPs allow us to construct low-dimensional (dimension 4) counter-examples indicating that there is essentially no room for further improvement of the rates that we provide.

All the numerical results in this paper were obtained on the MIT Supercloud Computing Cluster with Intel-Xeon-Platinum-8260 processor with 48 cores and 128 GB of RAM running Ubuntu 18.04.6 LTS with Linux 4.14.250-llgrid-10ms kernel . We used JuMP—a domain specific modeling language for mathematical optimization embedded in the open-source programming language Julia —to model the optimization problems. To solve the optimization problems, we use the following solvers: Mosek 9.3 , KNITRO 13.0.0 , and Gurobi 10.0.0, which are free for academic use. The relative feasibility tolerance and relative optimality tolerance of all the solvers are set at 1e-6. We validated the “bad” worst-case scenarios produced by our methodology using the PEPit package , which is an open-source Python library allowing to use the PEP framework.

The code used to generate and validate the results in this paper is available at:

https://github.com/Shuvomoy/NCG-PEP-code.

Related works

Conjugate gradient (CG) methods are particularly popular choices for solving systems of linear equations and quadratic minimization problems; in this context, they are known to be information-optimal in the class of first-order methods [33, Chapter 12 & Chapter 13] or [34, Chapter 5]. There are many extensions beyond quadratics, commonly referred to as nonlinear conjugate gradient methods (NCGMs). They are discussed at length in the textbooks [35, Chapter 5 & Chapter 7] and [36, Chapter 5] and in the nice survey . In particular, when exact line searches are used, many variants become equivalent and can be seen as instances of quasi-Newton methods, see [35, Chapter 7, §“Relationship with conjugate gradient methods”] or [36, Chapter 5, §5.5]. For instance, it is well known that standard variants such as Hestenes-Stiefel and Dai-Yuan are equivalent to ($\mathcal{M}$) when exact line searches are used, while being different in the presence of more popular line search procedures (such as Wolfe’s [35, Chapter 3]). Beyond quadratics, obtaining convergence guarantees is often reduced to the problem of ensuring the search direction to be a descent direction, see for instance [34, §5.5 “Extensions to non-quadratic problems”] or . Without exact line searches, even when $f$ is strongly convex, there are counter-examples showing that even popular variants may not generate descent directions . Note that NCGMs are often used together with restart strategies, which we do not consider here; see, e.g., and the references therein. Also, in [40, §5], the authors empirically demonstrate that NCGMs work very well in training deep learning problems.

In this work, we use the performance estimation framework . This methodology is essentially mature for analyzing “fixed-step” (i.e., non-adaptive) first-order methods (and for methods whose analyses are amenable to those of fixed-step methods), whose stepsizes are essentially chosen in advance. This type of methods include many common first-order methods and operator splitting schemes, including the heavy-ball method and Nesterov’s accelerated gradient . Only very few adaptive methods were studied using the PEP methodology, namely gradient descent with exact line searches , greedy first-order methods , and Polyak stepsizes . A premise to the study of NCGMs using PEPs was done in [26, §4.5.2], where an upper bound on the worst-case $\nicefrac{{(f(x_{2})-f_{\star})}}{{(f(x_{0})-f_{\star})}}$ of FR was numerically computed for two iterations and two condition number values, $q=0.1$ and $q=0.01$, where $q\triangleq\nicefrac{{\mu}}{{L}}$. This was achieved by numerically solving an SDP relaxation through a grid search on $\beta_{k}$. In comparison, we compute the worst-case $\nicefrac{{(f(x_{N})-f_{\star})}}{{(f(x_{0})-f_{\star})}}$ by solving the nonconvex PEPs associated with both FR and PRP to global optimality across a broader range of condition numbers over $q\in$ for $N=1,2,3,4$. Furthermore, for both methods, we also compute “Lyapunov”-type bounds on $\nicefrac{{(f(x_{k+N})-f_{\star})}}{{(f(x_{k})-f_{\star})}}$ that holds for any $k$ for $N=1,2,3,4$, and also establish their analytical complexity bounds offering a more comprehensive understanding of their performance. Our work is also closely related in spirit with the technique developed in for optimizing coefficients of fixed-step first-order methods using nonconvex optimization.

Preliminaries

In this section, we recall the definition and a result on smooth strongly convex functions, as well as a base result on steepest descent with an exact line search.

We simply denote $f\in\mathcal{F}_{\mu,L}$ when the dimension is either clear from the context or unspecified. We also denote by $q\triangleq\frac{\mu}{L}$ the inverse condition number. For readability, we do not explicitly treat the (trivial) case $L=\mu$.

Smooth strongly convex functions satisfy many inequalities, see e.g., [44, Theorem 2.1.5]. For the developments below, we need only one specific inequality characterizing functions in $\mathcal{F}_{\mu,L}$. The following result can be found in [10, Theorem 4] and is key in our analysis.

Another related result from [45, §2.1] that we record next involves constructing a strongly-convex smooth function from a given set of triplets.

Consider a function $f\in\mathcal{F}_{\mu,L}$ and the approximate steepest descent method:

where the search direction $d_{k}$ satisfies a relative error criterion:

Note that the relative tolerance $\epsilon$ needs to satisfy $\epsilon\in[0,1)$ for ($\mathcal{ASD}$) to converge. If $\epsilon\geqslant 1$, then $d_{k}=0$ becomes feasible and ($\mathcal{ASD}$) cannot be guaranteed to converge anymore, because in such a case we can select $d_{k}$ to be orthogonal to $\nabla f(x_{k})$ in practice [43, $\mathsection$5].

The iterates of ($\mathcal{ASD}$) satisfies the following two necessary (weaker) conditions for $x_{k+1}$ to follow ($\mathcal{ASD}$):

where the first condition follows from optimality of $\gamma_{k}$ in the line search condition in ($\mathcal{ASD}$) as follows

and the second condition comes from putting $d_{k}=\nicefrac{{(x_{k}-x_{k+1})}}{{\gamma_{k}}}$ in (4).

We will use the following result in our analysis. Note that similar results (without line searches) to that of Theorem 1.3 can be found in , which might help in future analyses of NCGMs without exact line searches.

where $q_{\epsilon}\triangleq\nicefrac{{\mu(1-\epsilon)}}{{L(1+\epsilon)}}$.

Next, we show that the relative error criterion ($\mathcal{REC}$) can be interpreted in simple geometric fashion in the context of exact line searches in ($\mathcal{ASD}$). In short, by letting $\theta_{k}$ be the angle between $\nabla f(x_{k})$ and $d_{k}$, ($\mathcal{REC}$) is equivalent to requiring $|\sin\theta_{k}|\leqslant\epsilon$. This fact is relatively simple to show and is stated without a proof in [43, §5]. Because we use this equivalence in the next section, we provide an elementary proof for completeness.

where $q_{\epsilon}\triangleq\nicefrac{{\mu(1-\epsilon)}}{{L(1+\epsilon)}}$.

Without loss of generality we can let $\theta_{k}$ to be acute, because $|\sin\theta_{k}|=|\sin(\pi-\theta_{k})|$. Now, consider the following method, where the search direction $d_{k}$ in ($\mathcal{ASD}$) is scaled by some factor $\alpha\neq 0$ with the scaled search direction denoted by $d_{k}^{\prime}=\alpha d_{k}$:

and we denote $\theta_{k}^{\prime}$ to be the angle between $\nabla f(x_{k})$ and $d_{k}^{\prime}$. We now show that ($\mathcal{ASD}$) and ($\mathcal{ASD}_{\textup{scaled}}$) are equivalent in the sense that they generate an identical sequence of iterates $x_{k},x_{k+1}$ along $|\sin\theta_{k}|=|\sin\theta_{k}^{\prime}|$. This is so because

i.e., the optimal stepsize $\gamma_{k}^{\prime}$ in ($\mathcal{ASD}_{\textup{scaled}}$) is the optimal step-size $\gamma_{k}$ in $\mathcal{ASD}$ scaled by $\nicefrac{{1}}{{\alpha}}$, leading to

Hence to establish our convergence result (6), we can work with ($\mathcal{ASD}_{\textup{scaled}}$). Next, we carefully select a nonzero $\alpha$ that ensures $\left\langle d_{k}^{\prime};\,d_{k}^{\prime}-\nabla f(x_{k})\right\rangle=0$, i.e., $d_{k}^{\prime}-\nabla f(x_{k})$ would be perpendicular to $d_{k}^{\prime}$ (see Figure 2); this yields $\alpha=\nicefrac{{\left\langle\nabla f(x_{k});\,d_{k}\right\rangle}}{{\|d_{k}\|^{2}}}$ which is nonzero because $\epsilon\in[0,1)$ implies $\left\langle\nabla f(x_{k});\,d_{k}\right\rangle\neq 0$. For this value of $\alpha$, we have $|\sin\theta_{k}|=\nicefrac{{\|d_{k}^{\prime}-\nabla f(x_{k})\|}}{{\|\nabla f(x_{k})\|}}$, which can be shown geometrically in Figure 2 in the right triangle (colored red) involving $\nabla f(x_{k})$, $d_{k}^{\prime}$, and $d_{k}^{\prime}-\nabla f(x_{k})$.

Now we are given that $|\sin\theta_{k}|\leqslant\epsilon$, hence setting $\alpha=\nicefrac{{\left\langle\nabla f(x_{k});\,d_{k}\right\rangle}}{{\|d_{k}\|^{2}}}$ ensures that the relative error criterion $\nicefrac{{\|d_{k}^{\prime}-\nabla f(x_{k})\|}}{{\|\nabla f(x_{k})\|}}\leqslant\epsilon$ is satisfied for ($\mathcal{ASD}_{\textup{scaled}}$). Finally by applying Theorem 1.3 to ($\mathcal{ASD}_{\textup{scaled}}$), we arrive at (6). ∎

Base descent properties of NCGMs

In this section, we analyze NCGMs as approximate steepest descent methods satisfying ($\mathcal{ASD}$) through a computer-assisted approach, where only the generated search directions matter, and not their magnitudes. This renders the analysis somewhat simpler, and we argue that this is a reasonable setting for improving the analysis and understanding of NCGMs.

This section builds on the intuition that when $|\sin\theta_{k}|$, where $\theta_{k}$ is the angle between the gradient and the search direction $d_{k}$ at iteration $k$, is upper bounded in an appropriate fashion, one can use Theorem 1.3 for obtaining convergence guarantees. In particular, we get nontrivial convergence guarantees as soon as $\theta_{k}$ can be bounded away from $\pm\frac{\pi}{2}$, i.e., $\sin\theta_{k}$ should be bounded away from $1$ for ensuring that $d_{k}$’s are descent directions. Of course, viewing NCGMs as approximate steepest descent methods is adversarial by nature, as it misses the point that the directions of NCGMs are meant to be better than those of vanilla gradient descent, while such analyses can only provide worse rates. Additionally, in Section (2.1), we provide additional justification behind analyzing NCGMs as approximate steepest descent methods through the lens of performance estimation problem (PEP), where we formulate the process of computing the worst-case $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ as optimization problems.

Albeit being pessimistic by construction, the analyses of this section are, to the best of our knowledge, novel for FR (for which we provide the first non-asymptotic convergence bound) and significantly better than the state-of-the-art bound for PRP. Furthermore, we show in Section 3.2 and Section 3.3 that there is actually nearly no room for improving those analyses.

Before going into the detailed approach, let us review a few properties of the iterates of ($\mathcal{M}$). Note that the iterates of ($\mathcal{M}$) satisfy the following equalities:

where the first two equalities are the same as ($\mathcal{ASD}_{\textup{relaxed}}$) following from exact line search. The last equality in (7) follows from applying the first equality to

Combining (8) with $\langle\nabla f(x_{k});d_{k}\rangle=\|\nabla f(x_{k})\|\|d_{k}\|\cos\theta_{k}$, we obtain that $\nicefrac{{\|\nabla f(x_{k})\|}}{{\|d_{k}\|}}=\cos\theta_{k}$, thereby reaching $\sin^{2}\theta_{k}=1-\nicefrac{{\|\nabla f(x_{k})\|^{2}}}{{\|d_{k}\|^{2}}}$. If we have $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}\leqslant c$ ($c\geqslant 1$ due to (7)), then $\sin^{2}\theta_{k}=1-\left(\nicefrac{{\|\nabla f(x_{k})\|^{2}}}{{\|d_{k}\|^{2}}}\right)\leqslant 1-\left(\nicefrac{{1}}{{c}}\right)$, yielding

The first two equations of (7), in conjunction with (9), satisfied by NCGMs, correspond to the same set of conditions required to apply Theorem 1.3. Thus, if we can establish an upper bound for the ratio $\nicefrac{{\|d_{k}\|}}{{||\nabla f(x_{k})||}}$ in the context of NCGMs, we can translate this into their worst-case convergence rates using Theorem 1.3.

In Section (2.1), we provide PEP-based perspective behind analyzing NCGMs as methods satisfying ($\mathcal{ASD}$). Section 2.2, first frames the problems of computing the worst-case $\nicefrac{{\|d_{k}\|}}{{\|\nabla f(x_{k})\|}}$ for PRP and FR as optimization problems for obtaining the desired bounds measuring the quality of the angle $\theta_{k}$ as PEPs. These PEPs are nonconvex but practically tractable QCQPs and can be solved numerically to certifiable global optimality using spatial branch-and-bound algorithms (detailed in Appendix D), which allows (i) to construct “bad” functions serving as counter-examples on which the worst-case $\nicefrac{{\|d_{k}\|}}{{\|\nabla f(x_{k})\|}}$ for PRP and FR is achieved, and (ii) to identify closed-form solutions to the PEPs leading to proofs that can be verified in a standard and mathematically rigorous way. The convergence rates for PRP and FR are provided and proved in Section 2.3.

A PEP perspective behind viewing NCGMs as approximate steepest descent method

In this section, we provide a PEP-based perspective behind analyzing NCGMs as approximate steepest descent methods satisfying ($\mathcal{ASD}$) through the lens of PEP. In this PEP approach, we formulate the problems of computing the worst-case ratios of $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ as the following optimization problem:

In Section (3.1), we will illustrate how we can formulate and solve (10) by casting it as a nonconvex QCQP. Note that in (10), the second constraint corresponds to third equation of (7) and the third constraint $\|d_{k}\|^{2}\leqslant c\|\nabla f(x_{k})\|^{2}$ models that if $\nabla f(x_{k})=0$ then $d_{k}=0$ for ($\mathcal{M}$). Note that $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}\geqslant 1$ because $\|\nabla f(x_{k})\|^{2}\leqslant\|d_{k}\|^{2}$, which follows from applying Cauchy–Schwarz inequality to (8).

While solving the nonconvex QCQPs equivalent to (10) for different values of $c$, $\mu$, and $L$, we found that the worst-case $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ is strictly monotonically increasing in $c$. Naturally, assigning an arbitrary value to $c$ would not reasonable to get the best bound, because the search direction generated by ($\mathcal{M}$) may not admit such a value. For example, for PRP, $c$ is always upper bounded by $1+\nicefrac{{L^{2}}}{{\mu^{2}}}$ as $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}\leqslant 1+\nicefrac{{L^{2}}}{{\mu^{2}}}$ for PRP [1, Theorem 2]. As we are interested in obtaining the tightest upper bound on $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$, the natural question is: What is the smallest admissible value of $c$, i.e., what is the least upper bound on the ratio $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}$ generated by ($\mathcal{M}$)? To that end, we numerically computed the least upper bound on $c$ by solving a problem similar to (10), except we replaced the objective $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ with $\nicefrac{{\|d_{k+1}\|^{2}}}{{\|\nabla f(x_{k+1})\|^{2}}}$ and then replaced the indices $k,k+1$ with $k-1,k$, respectively. In Section 2.2, we provide the details on formulating the problems of computing the worst-case ratios of $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}$ as nonconvex QCQPs. After we computed the least upper bound on $c$ numerically, we put them in (10). We then solved the associated nonconvex QCQP to global optimality, which numerically provided us with the tightest upper bound on worst-case $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$. Remarkably, at this stage, we found that these numerically computed worst-case $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ for ($\mathcal{M}$) exactly matched the analytical bound prescribed in Corollary 1.1. This PEP-based observation provides us a justification for analyzing NCGMs as approximate steepest descent methods and demonstrates the viability of this approach.

Computing worst-case search directions

In this section, we formulate the problems of computing the worst-case ratios of $\nicefrac{{\|d_{k}\|}}{{\|\nabla f(x_{k})\|}}$ as optimization problems. Following the classical steps introduced in , we show that it can be cast as a nonconvex QCQP.

For doing that, we assume that at iteration $k-1$ the NCGM has not reached optimality, so $\nabla f(x_{k-1})\neq 0.$ Because $\|\nabla f(x_{k-1})\|^{2}\leqslant\|d_{k-1}\|^{2}$ (follows from applying Cauchy–Schwarz inequality to (8)), without loss of generality we define the ratio $c_{k-1}\triangleq\nicefrac{{\|d_{k-1}\|^{2}}}{{\|\nabla f(x_{k-1})\|^{2}}}$ where $c_{k-1}\geqslant 1$. Then, denoting by $c_{k}$ the worst-case ratio $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}$ arising when applying ($\mathcal{M}$) to the minimization of an $L$-smooth $\mu$-strongly convex function, we will compute $c_{k}$ as a function of $L$, $\mu$, and $c_{k-1}$. In other words, we use a Lyapunov-type point of view and take the stand of somewhat forgetting about how $d_{k-1}$ was generated (except through the fact that it satisfies (7)). Then, we compute the worst possible next search direction $d_{k}$ that the algorithm could generate given that $d_{k-1}$ satisfies a certain quality. Thereby, we obtain an upper bound on the evolution of the quality of the search directions (quantified by $c_{k}$) obtained throughout the iterative procedure. Formally, we compute

For computing $c_{k}(\mu,L,c_{k-1})$, we reformulate (11) as follows. Denote $I\triangleq\{k-1,k\}$. An appropriate sampling of the variable $f$ (which is inconveniently infinite-dimensional) allows us to cast (11) as:

Using Theorem 1.1, the existence constraint can be replaced by a set of linear/quadratic inequalities (2) for all pairs of triplets in $\{(x_{i},g_{i},f_{i})\}_{i\in I}$ without changing the objective value. So, applying Theorem 1.1 to (12) followed by an homogeneity argument and a few substitutions based on (7), we arrive at:

Note that without the variable $n$ this problem is amenable to a finite-dimensional nonconvex QCQP (see Appendix B). Fortunately standard arguments (e.g., [10, Theorem 5]) allows setting $n=4$ without changing the optimal value of this problem, thereby discarding this dimension issue; we show this in Appendix B. We can then solve the QCQP equivalent to ($\mathcal{D}$) to certifiable global optimality using a custom branch-and-bound algorithm. Reformulation details are provided in Appendix B, whereas a description of the custom spatial branch-and-bound algorithm is given in Appendix D.

Finally, we recall that numerical solutions to ($\mathcal{D}$) correspond to worst-case functions that can be obtained through the reconstruction procedure from Theorem 1.2. In addition, numerical solutions can serve as inspirations for devising rigorous mathematical proofs, as presented next.

Worst-case bounds for PRP and FR

In this section, we provide explicit solutions to ($\mathcal{D}$) for PRP and FR. Those results are then used for deducing simple convergence bounds through a straightforward application of Theorem 1.3.

Solving ($\mathcal{D}$) with $\eta=1$ to global optimality allows obtaining the following worst-case bound for PRP quantifying the quality of the search direction with respect to the gradient direction.

with $q\triangleq\nicefrac{{\mu}}{{L}}$. Equivalently, $|\sin\theta_{k}|\leqslant\epsilon$ holds, where $\theta_{k}$ is the angle between $\nabla f(x_{k})$ and $d_{k}$ and $\epsilon=\nicefrac{{(1-q)}}{{(1+q)}}$.

Recall that $x_{k}=x_{k-1}-\gamma_{k-1}\,d_{k-1}$ and $d_{k}=\nabla f(x_{k})+\beta_{k-1}d_{k-1}$. The proof consists of the following weighted sum of inequalities:

optimality condition of the line search, with weight $\lambda_{1}=-\beta_{k-1}^{2}\frac{1+q}{L\gamma_{k-1}q}$:

smoothness and strong convexity of $f$ between $x_{k-1}$ and $x_{k}$, with weight $\lambda_{2}=\frac{\beta_{k-1}^{2}(1+q)^{2}}{L\gamma_{k-1}^{2}(1-q)q}$:

smoothness and strong convexity of $f$ between $x_{k}$ and $x_{k-1}$, with weight $\lambda_{3}=\lambda_{2}$:

definition of $\beta_{k-1}$ with weight $\lambda_{4}=\frac{\beta_{k-1}(1+q)}{L\gamma_{k-1}q}$:

which can be reformulated exactly as (expand both expressions and observe that all terms match)

thereby arriving at (13). Finally, using (9), we have $|\sin\theta_{k}|\leqslant\epsilon$ where $\epsilon=\nicefrac{{(1-q)}}{{(1+q)}}$. ∎

The following rate is a direct consequence of Lemma 2.1 and Theorem 1.3. Perhaps surprisingly, the following guaranteed convergence rate for PRP corresponds to that of gradient descent with an exact line search (Theorem 1.3 with $\epsilon=0$) when the condition number is squared.

The desired claim is a direct consequence of Corollary 1.1 with $\epsilon=\frac{1-q}{1+q}$. That is, the PRP scheme can be seen as a descent method with direction $d_{k}$ satisfying $\|d_{k}-\nabla f(x_{k})\|\leqslant\epsilon\|\nabla f(x_{k})\|$. ∎

As a take-away from this theorem, we obtained an improved bound on the convergence rate of PRP, but possibly not in the most satisfying way: this analysis strategy does not allow beating steepest descent. Furthermore, this bound is tight for one iteration assuming that the current search direction satisfies $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}=\nicefrac{{(1+q)^{2}}}{{4q}}$. However, it does not specify whether such an angle can be achieved on the same worst-case instances as those where Theorem 1.3 is achieved. In other words, there might be no worst-case instances where the bounds (6) and (13) are tight simultaneously, possibly leaving room for improvement in the analysis of PRP. We show in Section 3 that we could indeed slightly improve this bound by taking into account the history of the method in a more appropriate way by examining multiple iterations of ($\mathcal{M}$) rather than a single one.

The only worst-case complexity result that we are aware of in the context of PRP for smooth strongly convex problems was provided by Polyak in [1, Theorem 2]:

Figure 3 shows that the upper bound on $\nicefrac{{f(x_{k+1})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ for PRP (for different values of $q$) provided by [1, Theorem 2] is significantly worse compared to that of Theorem 2.1. From what we can tell, this is due to two main weaknesses in the proof of Polyak [1, Theorem 2]: a weaker analysis of gradient descent, and a weaker analysis of the direction (and in particular that $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}\leqslant 1+\nicefrac{{1}}{{q^{2}}}$). That is, whereas gradient descent with exact line searches is guaranteed to achieve an accuracy $f(x_{k})-f_{\star}\leqslant\varepsilon$ in $O(\nicefrac{{1}}{{q}}\log\nicefrac{{1}}{{\varepsilon}})$, our analysis provides an $O(\nicefrac{{1}}{{q^{2}}}\log\nicefrac{{1}}{{\varepsilon}})$ guarantee for PRP, where Polyak’s guarantee for PRP is $O(\nicefrac{{1}}{{q^{3}}}\log\nicefrac{{1}}{{\varepsilon}})$. As a reference, note that the lower complexity bound (achieved by a few methods, including many variations of Nesterov’s accelerated gradients) is of order $O(\sqrt{\nicefrac{{1}}{{q}}}\log\nicefrac{{1}}{{\varepsilon}})$.

A worst-case bound for Fletcher-Reeves (FR)

Similar to the obtaining of the bound for PRP, our bound for FR follows from solving ($\mathcal{D}$) (for $\eta=0$) in closed-form. We start by quantifying the quality of the search direction with respect to the steepest descent direction. Unlike PRP, where the worst-case ratio $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}$ depends only on the condition number $q$, in FR, the ratio $\nicefrac{{\|d_{k}\|^{2}}}{{\|\nabla f(x_{k})\|^{2}}}$ depends also on the previous ratio $\nicefrac{{\|d_{k-1}\|^{2}}}{{\|\nabla f(x_{k-1})\|^{2}}}$. To show this dependence, we first establish the following bound on the FR update parameter $\beta_{k-1}$ in terms of $\nicefrac{{\|d_{k-1}\|^{2}}}{{\|\nabla f(x_{k-1})\|^{2}}}$ and $q$.

where $q\triangleq\nicefrac{{\mu}}{{L}}$.

First, note that $\beta_{k-1}\geqslant 0$ by definition. The other part of the proof consists of the following weighted sum of inequalities:

relation between $\nabla f(x_{k-1})$ and $d_{k-1}$ with weight $\lambda_{1}=\gamma_{k-1}(L+\mu)-\frac{2\sqrt{\beta_{k-1}}}{\sqrt{(c_{k-1}-1)c_{k-1}}}$:

optimality condition of the line search with weight $\lambda_{2}=\frac{2}{c_{k-1}}-\gamma_{k-1}(L+\mu)$:

definition of $\beta_{k-1}$ with weight $\lambda_{3}=\frac{\sqrt{c_{k-1}-1}}{\sqrt{\beta_{k-1}c_{k-1}}}$:

initial condition on the ratio $\frac{\|d_{k-1}\|^{2}}{\|\nabla f(x_{k-1})\|^{2}}$ with weight $\lambda_{4}=-\gamma_{k-1}^{2}L\mu+\frac{\sqrt{\beta_{k-1}}}{c_{k-1}\sqrt{(c_{k-1}-1)c_{k-1}}}$ :

smoothness and strong convexity of $f$ between $x_{k-1}$ and $x_{k}$, with weight $\lambda_{5}=L-\mu$:

smoothness and strong convexity of $f$ between $x_{k}$ and $x_{k-1}$, with weight $\lambda_{6}=\lambda_{5}$:

which can be reformulated exactly as (expand the expressions and observe that all terms match):

Because, $-c_{k-1}\gamma_{k-1}^{2}L\mu+\gamma_{k-1}(L+\mu)-1$ is a concave function in $\gamma_{k-1},$ its maximum can be achieved by differentiating the term with respect to $\gamma_{k-1},$ equating it to , and then solving for $\gamma_{k-1}$. The corresponding maximum value is equal to $\nicefrac{{(L+\mu)^{2}}}{{4c_{k-1}L\mu}}-1$ and achieved at $\gamma_{k-1}=\nicefrac{{L+\mu}}{{2c_{k-1}L\mu}}$. Hence, the last inequality becomes:

Thereby, squaring both sides (which are nonnegative) of the last inequality and then through some algebra, we reach

As $\beta_{k-1}\geqslant 0$ by definition, we have thus proven the desired statement. ∎

Next, we prove a bound quantifying the quality of the search directions of FR.

Equivalently, $|\sin\theta_{k}|\leqslant\epsilon$ holds, where $\theta_{k}$ is the angle between $\nabla f(x_{k})$ and $d_{k}$ holds with $\epsilon=\sqrt{1-\nicefrac{{1}}{{c_{k}}}}$.

The proof consists of the following weighted sum of inequalities:

optimality condition of the line search with weight $\lambda_{1}=2\beta_{k-1}$:

the quality of the search direction with weight $\lambda_{2}=\beta_{k-1}^{2}$:

definition of $\beta_{k-1}$ with weight $\lambda_{3}=-c_{k-1}\beta_{k-1}$:

where in the last line we have used the upper bound on $\beta_{k-1}$ from (14). This gives us (15). Finally, using (9), we have $|\sin\theta_{k}|\leqslant\epsilon$, where $\epsilon=\sqrt{1-\nicefrac{{1}}{{c_{k}}}}$. ∎

That being said, this bound only allows obtaining unsatisfactory convergence results for FR, although not letting much room for improvements, as showed in the next sections.

with $\epsilon_{k}=\sqrt{\nicefrac{{(1-q)^{2}(k-1)^{2}}}{{4q+(1-q)^{2}(k-1)^{2}}}}$.

The desired claim is a direct consequence of Corollary 1.1 with Lemma 2.3. Indeed, it follows from

(the guarantee from Lemma 2.3 for the quality of the search direction) which we can rewrite as

with $c_{0}-1=0$, thereby arriving to $c_{k}\leqslant 1+k^{2}\nicefrac{{(1-q)^{2}}}{{4q}}$ by recursion. For applying Theorem 1.3, we compute $\epsilon_{k}=\sqrt{1-\nicefrac{{1}}{{c_{k}}}}\leqslant\sqrt{\nicefrac{{(1-q)^{2}k^{2}}}{{4q+(1-q)^{2}k^{2}}}}$ and reach the desired statement. ∎

It is clear that the statement of Theorem 2.2 is rather disappointing, as the convergence rate of the FR variation can become arbitrarily close to 1. While this guarantee clearly does not give a total and fair picture of the true behavior of FR in practice, it seems in line with the practical necessity to effectively restart the method as it runs .

The next section is devoted to studying the possibilities for obtaining tighter guarantees for PRP and FR beyond the simple single-iteration worst-case analyses of this section (which are tight for one iteration, but not beyond), showing that we cannot hope to improve the convergence rates for those methods without further assumptions on the problems at hand.

Obtaining better worst-case bounds for NCGMs

In the previous section, we established closed-form bounds on ratios between consecutive function values for NCGMs by characterizing worst-case search directions. Albeit being tight for the analysis of NCGMs for one iteration, the bounds that we obtained are disappointingly inferior to those of the vanilla gradient descent. In this section, we investigate the possibility of obtaining better worst-case guarantees for NCGMs. For doing this using our framework, one natural possibility for us is to go beyond the study of a single iteration (since our results appear to be tight for this situation). Therefore, in contrast with the previous section, we now proceed only numerically and provide worst-case bounds without closed-forms.

More precisely, we solve the corresponding PEPs in two regimes. In short, the difference between the two regimes resides in the type of bounds under consideration.

The first type of bounds can be thought to as a “Lyapunov” approach which studies $N$ iterations of ($\mathcal{M}$) starting at some iterate $(x_{k},d_{k})$ (for which we “neglect” how it was generated). In this first setup, we numerically compute worst-case bounds on $\nicefrac{{f(x_{k+N})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ for different values of $N$ (namely $N=1,2,3,4$). As for the results of Section 2, we quantify the quality of the couple $(x_{k},d_{k})$ by requiring that $\|d_{k}\|^{2}\leqslant c_{k}\|\nabla f(x_{k})\|^{2}$. When $N=1$, this setup corresponds to that of Section 2. Stemming from the fact the worst-case behaviors observed for $N=1$ might not be compatible between consecutive iterations, we expect the quality of the bounds to improve with $N$. Of course, the main weakness of this approach is the fact that we neglect how $(x_{k},d_{k})$ was generated.

As a natural complementary alternative, the second type of bounds studies $N$ iterations of ($\mathcal{M}$) initiated at $x_{0}$ (with $d_{0}=\nabla f(x_{0})$). Whereas the first type of bounds is by construction more conservative, it has the advantage of being recursive: it is valid for all $k\geqslant 0$. On the other side, the second type of bounds is only valid for the first $N$ iterations (the bound cannot be used recursively), but it cannot be improved at all. That is, we study exact worst-case ratio $\nicefrac{{f(x_{N})-f_{\star}}}{{f(x_{0})-f_{\star}}}$ for a few different values of $N$ (namely $N\in\{1,2,3,4\}$). In this setup, we obtain worst-case bounds that are only valid close to initialization. However, it has the advantage of being unimprovable, as we do not neglect how the search direction is generated.

This section is organized as follows. First, in Section 3.1 we present the performance estimation problems for ($\mathcal{M}$) specifically for computing the worst-case ratios $\nicefrac{{f(x_{k+N})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ and $\nicefrac{{f(x_{N})-f_{\star}}}{{f(x_{0})-f_{\star}}}$. Then, Section 3.2 and Section 3.3 presents our findings for respectively PRP and FR. Details on how we managed to solve the resulting nonconvex QCQPs numerically are provided in Appendix C. In Section 3.4, we discuss how to generate the counter-examples from the solutions to the nonconvex QCQPs.

Computing numerical worst-case scenarios

Similar to (11), the problem of computing the worst-case ratio $\nicefrac{{f(x_{k+N})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ is framed as the following nonconvex maximization problem (for $c\geqslant 1$ and $q\triangleq\nicefrac{{\mu}}{{L}}$):

We proceed similarly for $\nicefrac{{f(x_{N})-f_{\star}}}{{f(x_{0})-f_{\star}}}$:

Obviously, $\rho_{N}(q,c)\geqslant\rho_{N,0}(q)$ for any $c\geqslant 1$. We solve ($\mathcal{B}_{\textup{Lyapunov}}$) and ($\mathcal{B}_{\textup{exact}}$) numerically to high precision (details in Appendix C) for $N\in\{1,2,3,4\}$ and report the corresponding results in what follows. In the numerical experiments, we fix the values of $c$ using Lemma 2.1 for PRP in ($\mathcal{B}_{\textup{Lyapunov}}$), thereby computing $\rho_{N}\left(q,\nicefrac{{(1+q)^{2}}}{{4q}}\right)$ whose results are provided in Figure 4. For FR, $c$ can become arbitrarily bad and we therefore only compute $\rho_{N,0}(q)$ via ($\mathcal{B}_{\textup{exact}}$). The numerical values for $\rho_{N,0}(q)$ respectively PRP and FR are provided in Figure 5 and Figure 6. The next sections discuss and draw a few conclusions from the numerical worst-case convergence results for PRP and FR.

Improved worst-case bounds for PRP

Figure 4 reports the worst-case values of the “Lyapunov” ratio $\nicefrac{{f(x_{k+N})-f_{\star}}}{{f(x_{k})-f_{\star}}}$ as a function of the inverse condition number $q\triangleq\nicefrac{{\mu}}{{L}}$ and for $c=\nicefrac{{(1+q)^{2}}}{{4q}}$ and $N=1,2,3,4$. This worst-case ratio seems to improve as $N$ grows, but does not outperform gradient descent with exact line search (GDEL). The diminishing improvements with $N$ also suggests the worst-case performance of PRP in this regime might not outperform GDEL even for larger values of $N\geqslant 4$, albeit probably getting close to the same asymptotic worst-case convergence rate.

As a complement, Figure 5 shows how PRP’s worst-case ratio $\nicefrac{{f_{N}-f_{\star}}}{{f_{0}-f_{\star}}}$ evolves as a function of $q$ for $N=1,2,3,4$. The worst-case performance of PRP in this setup seems to be similar to that of GDEL. Further, for small $q$ (which is typically the only regime of interest for large-scale optimization), PRP’s worst-case performance seems to be slightly better than than of GDEL. On the other hand, for larger $q$, PRP performs slightly worse than GDEL.

As a conclusion, we believe there is no hope to prove uniformly better worst-case bounds for PRP than those for GDEL for base smooth strongly convex minimization. However, we might be able to prove improvements for small values of $q$ at the cost of possibly very technical proofs. As for the Lyapunov approach, the numerical results from this section could be improved by further increasing $N$, but we believe that the transient does not suggest this direction to be promising. We recall that we computed the bounds by solving an optimization problem whose feasible points correspond to worst-case examples. Therefore, the numerical results provided in this section are backed-up by numerically constructed examples on which PRP behaves “badly” (more details in Appendix C).

Improved worst-case bounds for FR

Figure 6 reports the worst-case values for the ratio $\nicefrac{{f_{N}-f_{\star}}}{{f_{0}-f_{\star}}}$ as a function of $q$, for $N\in\{1,2,3,4\}$. The convergence bounds appears to be marginally better than GDEL for some sufficiently small inverse condition numbers. Further, the range of values of $q$ for which there is an improvement appears to be decreasing with $N\geqslant 2$. Beyond this range, the worst-case values become significantly worse than that of GDEL. Though apparently not as dramatic as the worst-case bound from Theorem 2.2, the quality of the bound appears to be decreasing with $N$, which stands in line with the practical need to restart the method .

As in the previous section, we recall that those curves were obtained by numerically constructing “bad” worst-case counter-examples satisfying our assumptions. In other words, there is no hope to obtain better results without adding assumptions or changing the types of bounds under consideration.

Constructing counter-examples

Once we have solved the nonconvex QCQP associated with ($\mathcal{B}_{\textup{Lyapunov}}$) or ($\mathcal{B}_{\textup{exact}}$), we can construct the associated triplets $\{x_{i},g_{i},f_{i}\}_{i\in I_{N}^{\star}}$ and then apply Theorem 1.2 to construct the corresponding “bad” function. This “bad” function serves as a counter-example, illustrating scenarios where ($\mathcal{M}$) performs poorly. One can access the numerically constructed triplets $\{x_{i},g_{i},f_{i}\}_{i\in I_{N}^{\star}}$ associated with the counter-examples by following the instructions provided in our github repository. Next, we provide a concrete example of how to construct a “bad” function for ($\mathcal{B}_{\textup{Lyapunov}}$) from our provided code and datasets located in the folder titled ‘Code_for_NCG_PEP’ of the github repository. Constructing counter-examples for ($\mathcal{B}_{\textup{exact}}$) is analogous.

we can run ‘1.Example_Julia.ipynb’ with the input parameters and generate the function by solving the nonconvex QCQP directly and generate the triplets, or

we can directly access the triplets from the saved datasets in the folder Saved_Output_Files with instructions provided in the file ‘2.Using_the_saved_datasets_Julia.ipynb’.

For the sake of completeness, we provide the numerical values of $\{x_{i}\}_{i\in\{\star,0,1,2\}},\{g_{i}\}_{i\in\{\star,0,1,2\}},$ and $\{f_{i}\}_{i\in\{\star,0,1,2\}}$ of the function in this setup in Table 1, Table 2, and Table 3, respectively. From the numerical values of the triplets, we can construct the “bad” function using Theorem 1.2. For this constructed function, we have the performance guarantee $\nicefrac{{f(x_{k+2})-f_{\star}}}{{f(x_{k})-f_{\star}}}\geqslant 0.056104$, which closely matches the bound provided in Figure 4. Additionally, this guarantee can be verified through other (independent) open-source software ; we provide code for this independent verification in the file called ‘3.PEPIt_verification_Python.ipynb’.

Conclusion

This works studies the iteration complexity of two variants of nonlinear conjugate gradients, namely the Polak-Ribière-Polyak (PRP) and the Fletcher-Reeves (FR) methods. We provide novel complexity bounds for both those methods, and show that albeit unsatisfying, not much can a priori be gained from a worst-case perspective, as both methods appear to behave similar or worse to regular steepest descent in the worst-case. Further, those results suggest that explaining the good practical performances of NCGMs might be out of reach for traditional worst-case complexity analyses on classical classes of problems.

This work considers only somewhat “ideal” variants of nonlinear conjugate gradient methods, as we make explicit use of exact line search procedures. However, there is a priori no reason to believe that inexact line search procedures would improve the possibly bad worst-case behaviors. Further, the performance estimation methodology allows taking such inexact line search procedures into account, so the same methodology could be applied for tackling those questions. We leave such investigations for future work.

Acknowledgments

S. Das Gupta and R. M. Freund acknowledge support by AFOSR Grant No. FA9550-22-1-0356. A. Taylor acknowledges support from the European Research Council (grant SEQUOIA 724063). This work was partly funded by the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

The authors thank Nizar Bousselmi and Ian Ruffolo for careful reading of the manuscript and constructive feedback.

References

Organization of the appendix

In what follows, we report detailed numerical results and computations that are not presented in the core of the paper. Table 4 details the organization of this additional material.

Appendix A Tightness of the worst-case search directions

Figure 7 and Figure 8 illustrate the tightness of the bounds (13) and (15) for PRP and FR respectively. That is, we compare the numerical bounds (discrete points) with closed-forms (continuous lines) for a few different values of $q$ and $c_{k-1}$. Numerical bounds are obtained by solving ($\mathcal{D}$) with $\eta=1$ for PRP and $\eta=0$ for FR. These numerical examples strongly suggest that our bounds cannot be improved in general. Absolute relative differences between closed-form expressions and numerical ratios is less than $1\textrm{e}-6$ in all cases.

Appendix B Nonconvex QCQP reformulation of (𝒟𝒟\mathcal{D})

To reformulate ($\mathcal{D}$) as a nonconvex QCQP, we introduce the following Grammian matrices that is a common step in performance estimation literature :

This ensures that $x_{i}=H\mathbf{x}_{i}$, $g_{i}=H\mathbf{g}_{i}$, $d_{i}=H\mathbf{d}_{i}$, $f_{i}=F\mathbf{f}_{i},$ for all $i,j\in I$. Next, for appropriate choices of matrices $A_{i,j}$, $B_{i,j}$, $C_{i,j}$, $\widetilde{C}_{i,j}$, $D_{i,j}$, $\widetilde{D}_{i,j}$, $E_{i,j}$, and vector $a_{i,j}$, we can ensure that the following reformulations hold for all $i,j\in I$:

where $G$, $F$, $H$, $\Theta$, $\{\Theta_{i,j}\}_{i,j\in I}$, $\gamma_{k-1}$, $\beta_{k-1}$ are the decision variables. This nonconvex QCQP can be solved to certifiable global optimality using a custom spatial branch-and-bound algorithm described in Appendix D.

Similar to the reformulations from ($\mathcal{D}$), ($\mathcal{B}_{\textup{Lyapunov}}$) and ($\mathcal{B}_{\textup{exact}}$) can be cast as nonconvex QCQPs, where the number of nonconvex constraints grows quadratically with $N$. Thereby, solving them to global optimality in reasonable time for $N=3,4$ is already challenging.

Therefore, rather than solving the nonconvex QCQP reformulations of ($\mathcal{B}_{\textup{Lyapunov}}$) and ($\mathcal{B}_{\textup{exact}}$) directly, we compute upper bounds and lower bounds by solving more tractable nonconvex QCQP formulations. We then show that the relative gap between the upper and lower bounds is less than $10\%$ which thereby indicates that there is essentially no room for further improvement.

This section presents our upper bound $\overline{\rho}_{N}(q,c)$ and lower bound $\underline{\rho}_{N}(q,c)$ on $\rho_{N}(q,c)$.

Using (7), we have the following relaxation of ($\mathcal{B}_{\textup{Lyapunov}}$), which provides upper bounds on ${\rho}_{N}(q,c)$:

Using the notation $g_{i}\triangleq\nabla f(x_{i})$ and $f_{i}\triangleq f(x_{i})$ again, and then applying an homogeneity argument, we write (21) as:

where $f,n,\{x_{k+i}\}_{i\in[0:N]},\{d_{k+i}\}_{i\in[0:N]}$ are the decision variables. Define $I_{N}^{\star}=\{\star,k,k+1,\ldots,k+N\}$. Next, note that the equation $d_{k+i+1}=g_{k+i+1}+\beta_{k+i}d_{k+i}$ for $i\in[0:N-2],$ can be written equivalently as the following set of equations:

where we have introduced the intermediate variables $\chi_{j,i}$, which will aid us in formulating (22) as a nonconvex QCQP down the line. In absence of these intermediate variables in (23), the resultant constraints in the final optimization problem will involve polynomials of degree three or more in the decision variables, and such optimization problems present a significantly greater challenge in solving to global optimality compared to a QCQP. Next, using (23) and Theorem 1.1, we can equivalently write (22) as:

where $\{x_{k+i},g_{k+i},f_{k+i}\}_{i},n,\{d_{k+i}\}_{i},\{\beta_{k+i}\}_{i},\{\chi_{j,i}\}_{j,i}$ are the decision variables. Note that we have set $g_{\star}=0$, $x_{\star}=0$, and $f_{\star}=0$ without loss of generality, because both the objective and the function class are closed and invariant under shifting variables and function values. We introduce Grammian matrices again:

where $F,G,H,\{\chi_{j,i}\}_{j,i},\{\beta_{k+i}\}_{i}$ are the decision variables.

We now discuss how we can calculate $\underline{\rho}_{N}(q,c)$ and construct the corresponding “bad” function. This function serves as a counter-example, illustrating scenarios where ($\mathcal{M}$) performs poorly. Once we have solved (27), it provides us with the corresponding CG update parameters, which we denote by $\overline{\beta}_{i}$. If we can solve ($\mathcal{B}_{\textup{Lyapunov}}$) with the CG update parameters fixed to the $\overline{\beta}_{i}$ found from (27), then it will provide us with the lower bound $\underline{\rho}_{N}(\mu,L,c)$. This process also yields a “bad” function that acts as a counter-example, which we explain next. Using the notation $g_{i}\triangleq\nabla f(x_{i})$ and $f_{i}\triangleq f(x_{i})$, then applying the homogeneity argument, we can compute $\underline{\rho}_{N}(q,c)$ by finding a feasible solution to the following optimization problem:

where $f,n,\{x_{k+i}\},\{d_{k+i}\}_{i},\{\gamma_{k+i}\}_{i}$ are the decision variables. Next, note that the NCGM iteration scheme in (28) can be equivalently written as:

where we have introduced intermediate variables $\chi_{j,i}$ and $\alpha_{i,j}$ which will aid us in formulating (28) as a nonconvex QCQP. Define $I_{N}^{\star}=\{\star,k,k+1,\ldots,k+N\}$. Now using (29), Theorem 1.1, and (7), we can equivalently write (22) as:

where $\{x_{k+i},g_{k+i},f_{k+i}\}_{i},n,\{\gamma_{k+i}\}_{i},\{\chi_{j,i}\}_{j,i},\{\alpha_{i,j}\}_{i,j}$ are the decision variables. We introduce the Grammian transformation:

where $G,F,\{\Theta_{i,j}\}_{i,j\in I_{N}^{\star}},H,\gamma,\alpha,\chi$ are the decision variables. Note that $\{\Theta_{i,j}\}_{i,j\in I_{N}^{\star}}$ is introduced as a separate decision variable to formulate the cubic constraints arising from $B_{i,j}$ as quadratic constraints. Also, to compute $\underline{\rho}_{N}(q,c)$, it suffices to find just a feasible solution to (33), in Appendix D we will discuss how to do so using our custom spatial branch-and-bound algorithm.

Now we discuss how we compute the upper bound $\overline{\rho}_{N,0}(q)$ and lower bound $\underline{\rho}_{N,0}(q)$ to $\rho_{N,0}(q)$ defined in ($\mathcal{B}_{\textup{exact}}$). The bound computation process is very similar to that of ($\mathcal{B}_{\textup{Lyapunov}}$). Observe that, in ($\mathcal{B}_{\textup{Lyapunov}}$), if we remove the constraint $\|d_{k}\|^{2}\leqslant c\|\nabla f(x_{k})\|^{2},$ set $k\triangleq 0$ , and then add the constraint $d_{0}=\nabla f(x_{0})$, then it is identical to ($\mathcal{B}_{\textup{exact}}$) (the constraint $\left\langle\nabla f(x_{0});\,d_{0}\right\rangle=\|\nabla f(x_{0})\|^{2}$ in ($\mathcal{B}_{\textup{Lyapunov}}$) is a valid but redundant constraint for ($\mathcal{B}_{\textup{exact}}$)).

To compute the lower bound $\underline{\rho}_{N,0}(q)$, we follow the same set of changes described in the last paragraph but to (28) in Section C.1.2.

The relative gap between the lower bounds and upper bounds

Tables 5, 6, 7 record the relative gap between lower bounds and upper bounds for a few representative values of $q$ obtained by solving the aforementioned nonconvex QCQPs associated with ($\mathcal{B}_{\textup{Lyapunov}}$) and ($\mathcal{B}_{\textup{exact}}$) using our custom spatial branch-and-bound algorithm described in Appendix D. Note that the tables contain a few negative entries close to zero which are due to the absolute gap being of the same order as the accuracy of the solver ($1\textrm{e}-6$). For the full list for all values, we refer to our open-source code in Section 4, which also allows for computing these bounds for a user-specified value of $q$ as well. In all cases, the relative gap is less than $10\%$. In most cases, it is significantly better.

Appendix D Custom spatial branch-and-bound algorithm

This section discusses implementation details for solving the nonconvex QCQPs of this paper (namely (20), (27), or (33)) using a custom spatial branch-and-bound method. This strategy proceeds in three stages, as follows.

Stage 1: Compute a feasible solution. First, we construct a feasible solution to the nonconvex QCQP. We do that by generating a random $\mu$-strongly convex and $L$-smooth quadratic function, and by applying the corresponding nonlinear conjugate gradient method on it. The corresponding iterates, gradient and function values correspond to a feasible point for the nonconvex QCQPs under consideration.

Stage 2: Compute a locally optimal solution by warm-starting at Stage 1 solution. Stage 2 computes a locally optimal solution to the nonconvex QCQPs using an interior-point algorithm, warm-starting at the feasible solution produced by Stage 1. When a good warm-starting point is provided, interior-point algorithms can quickly converge to a locally optimal solution under suitable regularity conditions , [50, §3.3]. In the situation where the interior-point algorithm fails to converge, we go back to the feasible solution from Stage 1. We have empirically observed that Stage 2 consistently provides a locally optimal solution.

Stage 3: Compute a globally optimal solution by warm-starting at Stage 2 solution. Stage 3 computes a globally optimal solution to the nonconvex QCQP using a spatial branch-and-bound algorithm , warm-starting at the locally-optimal solution produced by Stage 2. For details about how spatial branch-and-bound algorithm works, we refer the reader to [27, §4.1].

In stage 3, the most numerically challenging nonconvex quadratic constraint in (20), (27) or (33) is $G=PP^{\top}$. To solve those problems in reasonable times, we use the lazy constraints approach, [27, §4.2.5].

In short, we replace the constraint $G=PP^{\top}$ by the infinite set of linear constraints $\mathop{\bf tr}\left(Gyy^{\top}\right)\geqslant 0$ for all $y$, which we then sample to obtain a finite set of linear constraints (we recursively add additional linear constraints afterwards if need be). More precisely, we use

where the initial $Y$ is generated randomly as a set of unit vectors following the methodology described in [53, $\mathsection$5.1]. By replacing $G=PP^{\top}$ by (34) we obtain a simpler (but relaxed) QCQP. Then, we update the solution $G$ lazily by repeating the following steps until $G\succcurlyeq 0$ is satisfied subject to a termination criterion. Practically speaking, our termination criterion is that the minimal eigenvalue of $G$ is larger than $\epsilon\approx-1\textrm{e}-6$; until then, we repeat the following procedure:

Solve the relaxation of the nonconvex QCQPs, where (34) is used instead of $G=PP^{\top}$, which provides us an upper bound on the original nonconvex QCQP.

Compute the minimal eigenvalue $\textup{eig}_{\textup{min}}(G)$ and the corresponding eigenvector $u$ of $G$. If $\textup{eig}_{\textup{min}}(G)\geq 0$, we reached an optimal solution to the nonconvex QCQP and we terminate.

If $\textup{eig}_{\textup{min}}(G)<0$, we add a constraint $\mathop{\bf tr}(Guu^{\top})\geqslant 0$ lazily, which makes the current $G$ infeasible for the new relaxation. We use the lazy constraint callback interface of JuMP to add constraints lazily, which means that after adding one additional linear constraint, updating the solution in step 1 is efficient since Gurobi and all modern solvers based on the simplex algorithm can quickly update a solution when only one linear constraint is added [54, pp. 205-207].