A single potential governing convergence of conjugate gradient, accelerated gradient and geometric descent
Sahar Karimi Stephen Vavasis
math.OC

[
    [
    December 30, 2023
=====================


Nesterov's accelerated gradient (AG)
method for minimizing a smooth strongly convex function f is
known to reduce f(_k)-f(^*) by a factor 
of ∈(0,1) after k=O(√(L/ℓ)log(1/)) iterations, where
ℓ,L are the two parameters of smooth strong convexity. Furthermore,
it is known that this is the best possible complexity in the function-gradient oracle
model of computation.  Modulo a line search, the geometric descent (GD)
method of Bubeck, Lee and Singh has the same bound for this class of functions.
The method of linear conjugate gradients (CG) 
also satisfies the same
complexity bound in the special case of strongly convex quadratic functions,
but in this special case it can be faster than the AG and GD methods.

Despite similarities in the algorithms and their
asymptotic convergence rates, the conventional analysis of the
running time of CG is mostly disjoint
from that of AG and GD.  The analyses of 
the AG and GD methods are also rather distinct.

Our main result is analyses of the three methods that share several
common threads: all three analyses show
a relationship to a certain “idealized algorithm”, all three
establish the convergence rate
through the use of the Bubeck-Lee-Singh geometric lemma, and all three
have the same potential
that is computable at run-time and exhibits decrease
by a factor of 1-√(ℓ/L) or better
per iteration.

One application of these analyses is that they open the possibility of hybrid
or intermediate algorithms.  One such algorithm is proposed herein
and is shown to perform well in  computational tests.








§ FIRST-ORDER METHODS FOR STRONGLY CONVEX FUNCTIONS

Three methods for minimizing smooth, strongly
convex functions are considered in this work, conjugate gradient,
accelerated gradient, and geometric descent. 
CG is the oldest and perhaps best known of the methods.  It
was introduced by Hestenes and Stiefel <cit.>
for 
minimizing strongly convex quadratic functions
of the form f()=^TA/2-^̱T, where A is
a symmetric positive definite matrix.

There is a significant body of work on
gradient methods for more general
smooth, strongly convex functions.  We say that a 
differentiable convex function f:^n→ is
smooth, strongly convex 
<cit.>
if there exist two scalars L≥ℓ>0 such that
for all ,∈^n,

    ℓ‖-‖^2/2≤ f()-f()-∇ f()^T(-)≤ L‖-‖^2/2.

This is equivalent to assuming convexity and lower and upper
Lipschitz constants on the gradient:

    ℓ‖-‖≤‖∇ f()-∇ f()‖≤
    L‖-‖.


Nemirovsky and Yudin <cit.>
proposed
a method for minimizing smooth strongly convex
functions requiring
k=O(√(L/ℓ)log(1/)) iterations to produce an iterate 
_k such that f(_k)-f(^*)≤(f(_0)-f(^*)), where
^* is the optimizer (necessarily unique under the assumptions made).  
A drawback of their method is that it requires a two-dimensional
optimization on each iteration that can
be cumbersome to implement (to the best of our knowledge,
the algorithm was not ever widely adopted).
Nesterov <cit.> proposed
another method, nowadays
known as the “accelerated gradient” (AG) method, 
which achieves the same optimal complexity that requires a single
function and gradient evaluation on each iteration.

In the special case of strongly convex quadratic functions, 
the parameters ℓ and L appearing in eq:strconvdef
correspond to λ_min(A) and λ_max(A), the extremal
eigenvalues of A.  
The conjugate
gradient method has already been known to satisfy the asymptotic
iteration bound k=O(√(L/ℓ)log(1/)) since the work
of Daniel (1967) described below.

Although the two methods satisfy the same asymptotic bound,
the analyses of the two methods are completely
different.  In the case of AG, there are
two analyses by Nesterov <cit.>.
In our own previous work <cit.>, we provided a third
analysis based on another potential.

In the case of linear
conjugate gradient, we are aware of no direct analysis of the algorithm
prior to our own previous work <cit.>.
By “direct,” we mean an analysis of f(_k)-f(^*) using the
recurrence inherent in CG.  Instead, the standard
analysis introduced by Daniel, whose theorem is
stated precisely below, proves that another iterative method, for
example Chebyshev iteration <cit.> or the heavy-ball iteration
<cit.> achieves reduction of 
(1-O(√(ℓ/L))) per iteration.  Then one appeals
to the optimality of the CG iterate in the Krylov space generated
by all of these methods to claim that the CG iterate must be at
least as good as the others.

Recently, Bubeck, Lee and Singh <cit.> proposed the
geometric descent (GD) algorithm for
analysis of a variant of accelerated
gradient <cit.>, called “geometric descent” (GD).
As presented by the authors, the algorithm requires an exact
line-search on each iteration, although it is possible that similar theoretical
guarantees could be established for an approximate
line search.  Under the assumption
that the line-search requires a constant number of function and 
gradient evaluations, then GD also requires
k=O(√(L/ℓ)log(1/)) iterations.

We propose analyses of these three algorithms that share several 
common features.  
First, all three algorithms 
can be analyzed using the geometric lemma
of Bubeck, Lee and Singh, which is presented in Section <ref>.
Second, all three are related to an “idealized” unimplementable
algorithm which is described and analyzed in Section <ref>.  
Finally, the convergence
behavior for all three of them is governed by a potential σ̃_k,
which has the following three properties:


  * There exists an auxiliary sequence of vectors _0,_1,… such that

    σ̃_k^2 ≥‖_k-^*‖^2 +2(f(_k)-f(^*))/ℓ,

for k=0,1,2…,

  * σ̃_k+1^2 ≤(1-√(ℓ/L))σ̃_k^2, and

  * σ̃_k is computable on each iteration (assuming prior knowledge of
  ℓ,L) in O(n) operations.  Our definition of “computable” is
  explained in more detail in Section <ref>.

These results are established for the GD algorithm in Section <ref>,
for the CG algorithm in
Section <ref>,
and for the 
AG algorithm in Section <ref>.

The relationship between IA and the three algorithms is explained in detail
as follows.
Section <ref> shows that CG exactly implements IA
for quadratic objective functions even though IA is in general unimplementable.
On the other hand, GD (analyzed in Section <ref>) and AG
(analyzed in Section <ref>)
both simulate IA in the sense
that they produce optimal iterates given  partial information
about the objective in the current iterate.

Because the three algorithms each compute a scalar σ̃_k satisfying
the above properties, it becomes straightforward to create hybrids.  In other
words, the above analysis treats all three algorithms as essentially 1-step
processes as opposed to long inductive chains.  In Section <ref> we
propose a hybrid CG algorithm that performs well in computational tests, which
are described in Section <ref>.
The reason from making a hybrid CG algorithm is that the performance
of linear conjugate gradient on specific instances can be much better than the
worst-case bound given by Daniel's theorem; the performance on specific
instances is highly governed by the eigenvalues of A.  Therefore, using
a conjugate-gradient-like algorithm for a nonlinear problem may also perform
better than the (1-√(ℓ/L))^k worst-case convergence bound.  This is
also the motivation for traditional nonlinear conjugate gradient, as we
discuss below.

We conclude this introductory section with a few remarks about our
previous related manuscript <cit.>.  In that work, we established
that a potential defined by

    Ψ_k=‖_k-^*‖^2 +2(f(_k)-f(^*))/ℓ

decreases by a factor (1-√(ℓ/L)) per iteration for both CG
and AG.  This Ψ_k is not computable since ^* is not known
a priori, and therefore our previous result does not have any
immediate algorithmic application.  Our notion of “computability”
is defined in more detail in Section <ref>.
The potential σ̃_k developed
herein is computable on every step and therefore may be used to guide
a hybrid algorithm.  In addition, the current work also applies to
the GD method, which was not addressed in our previous manuscript.



§ NOTATION

Define B(,r)={∈^n:‖-‖≤ r}, i.e., the
closed ball centered at ∈^n of radius r.

An affine set
is a set of the form ℳ={+̧:∈𝒲} where ∈̧^n is fixed
and 𝒲⊂^n is a linear subspace.
We write this as ℳ=+̧𝒲, a special case of a Minkowski sum.
Another notation for an affine set is
{_1,…,_k}, which is defined as
{α_1_1+⋯+α_k_k:α_1+⋯+α_k=1}.
If _1,…,_k span 𝒲, then
it is clear that 
+̧𝒲={,̧+̧_1,…,+̧_k}.

Suppose 𝒰 is an affine subset of ^n.
The set T𝒰={-:,∈𝒰} 
is called the tangent space
of 𝒰 and is a linear subspace.  If 𝒰 is presented as
𝒰=+̧𝒲, where 𝒲 is a linear subspace, then
it follows that T𝒰=𝒲.




§ PRELIMINARY LEMMAS



We start with a special case of
a lemma from Drusvyatskiy et al. <cit.>, which
is an extension of work by Bubeck et al. <cit.>:


Suppose ,∈^n.  Let
δ,ρ,σ be three nonnegative scalars
such that δ≤‖-‖.
Suppose λ∈[0,1] and 

    =(1-λ)+λ.

Then

    B(,ρ)∩ B(,σ)⊂ B(,ξ),

where

    ξ = √((1-λ)ρ^2+λσ^2-λ(1-λ)δ^2),

The
argument of the square-root in eq:xidef is guaranteed to be nonnegative
whenever B(,ρ)∩ B(,σ)∅,
or equivalently, whenever ρ+σ≥‖-‖.




We prove the second claim first.
The quantity appearing in the square root of eq:xidef is
nonnegative as the following inequalities show:

    (1-λ)ρ^2 + λσ^2 - (1-λ)λδ^2
       =(1-λ)λ(ρ^2 + σ^2 - δ^2) +
    (1-λ)^2ρ^2+λ^2σ^2 
       ≥ (1-λ)λ(ρ^2 + σ^2 - δ^2) +
    2(1-λ)λρσ
       = (1-λ)λ((ρ + σ)^2 - δ^2) 
       ≥ 0,

where the last line uses the assumptions ρ+σ≥‖-‖≥δ.

Now for the first part of the lemma,
the proof that B(,ρ)∩ B(,σ)⊂ B(,ξ)
for ξ given by eq:xidef follows from more general
analysis in Drusvyatskiy et al. <cit.>.  
Assume that ∈ B(,ρ)∩ B(,σ) so

    (-)^T(-) -ρ^2    ≤ 0, 
    
    (-)^T(-) -σ^2    ≤ 0.

For λ∈[0,1], add (1-λ) times eq:px1 to
λ times eq:py1 and rearrange to obtain a new inequality satisfied
by :

    (-)^T(-)+(1-λ)λ‖-‖^2 -(1-λ)ρ^2-
    λσ^2≤ 0,

i.e.

    ‖-‖≤((1-λ)ρ^2+λσ^2-(1-λ)λ‖-‖^2)^1/2,

where  is defined by eq:lemmaz.  By substituting the definition
δ≤‖-‖, we observe that ∈ B(,ξ), where ξ is defined
by eq:xidef.


This leads to the following, which is
a more precise statement of the geometric lemma from
Bubeck et al. <cit.>:


Let ,,ρ,σ,δ be as in the preceding lemma.
Under the assumption 
ρ+σ≥δ and the
additional assumption
δ≥√(|ρ^2-σ^2|), 
eq:xidef is minimized 
over possible choices of λ∈[0,1] by:

    λ^* = δ^2+ρ^2-σ^2/2δ^2,

yielding

    ^*=(1-λ^*)+λ^*,

in which case the minimum value of eq:xidef is,

    ξ^*=1/2√(2ρ^2+2σ^2-δ^2-
    (ρ^2-σ^2)^2/δ^2).





First, note that
|ρ^2-σ^2|/δ^2≤ 1 by the assumption made, thus
ensuring that λ^*∈[0,1].  
Therefore, it follows from the preceding lemma that
the quantity appearing in the square root of eq:xidef is nonnegative.
The previous lemma establishes that for any ∈ B(,ρ)∩ B(,σ),
and for an arbitrary λ∈[0,1],

    ‖-‖^2 ≤ (1-λ)ρ^2 + λσ^2 - λ(1-λ)δ^2.

We observe that the right-hand side is a convex quadratic in λ and
hence is minimized when the derivative with respect to λ is zero, and one
checks that this value is precisely eq:lambdastar.  Substituting 
λ=λ^* into eq:xidef yields eq:xistardef.





§ IDEALIZED ALGORITHM



We consider the following idealized algorithm for minimizing f(), where
f:^n→ is smooth, strongly convex.  As in the introduction,
let ℓ,L denote the two parameters of strong convexity.


    _0:=
       ℳ_1 := _0+{∇ f(_0)}
       k:=1,2,…
       _k:={f():∈ℳ_k}
       _k := {‖-^*‖:∈ℳ_k}
       ℳ_k+1:=_k+{_k-_k,∇ f(_k)}



This algorithm is called “idealized” because it is not implementable 
in the general case;
it requires prior knowledge of ^* in 
eq:ia2.yupd.  Nonetheless, we will argue
that CG, accelerated gradient, and geometric gradient are related
to the idealized algorithm in different ways. 

Notice that ℳ_k is an affine set that is two-dimensional
on most iterations.
Alternate notation for this set, also used herein, is
ℳ_k={_k-1,_k-1,_k-1-∇ f(_k-1)}.
Note also that by eq:ia2.xupd and eq:ia2.yupd,
_k,_k∈ℳ_k, and by eq:ia2.mupd,
_k,_k∈ℳ_k+1, and therefore ℳ_k,
ℳ_k+1 have a common 1-dimensional affine subspace.


We start with the main theorem about
IA.  For iteration k, define
a potential Ψ_k as follows:

    Ψ_k= ‖_k-^*‖^2+2(f(_k)-f(^*))/ℓ.



For Algorithm IA, for each k=1,2,…, 

    Ψ_k+1≤(1-√(ℓ/L))Ψ_k.





The proof follows closely from the analysis in <cit.>.  Define

    _k   =_k-∇ f(_k)/L, 
    _k   =_k-∇ f(_k)/ℓ.

The point _k
satisfies f(_k)≤ f(_k)-‖∇ f(_k)‖^2/(2L).
Observe that
_k∈ℳ_k+1,
so _k is a candidate
for the optimizer in eq:ia2.xupd on iteration k+1, and hence

    f(_k+1)≤ f(_k)-‖∇ f(_k)‖^2/(2L),

which is equivalent to

    2(f(_k+1)-f(^*))/ℓ≤2(f(_k)-f(^*))/ℓ-
    ‖∇ f(_k)‖^2/Lℓ.


Next, observe that a rearrangement of the definition of strong convexity yields:

    -2∇ f(_k)^T(_k-^*)/ℓ
    +‖_k-^*‖^2
    ≤-2(f(_k)-f(^*))/ℓ.

We use this result in the following:

    ‖_k-^*‖^2     = ‖_k-_k+_k-^*‖^2 
       = ‖_k-_k‖^2 + 2(_k-_k)^T(_k-^*) + ‖_k-^*‖^2 
       =‖∇ f(_k)‖^2/ℓ^2 -
    2∇ f(_k)^T(_k-^*)/ℓ +‖_k-^*‖^2 
       ≤‖∇ f(_k)‖^2/ℓ^2 
    -2(f(_k)-f(^*))/ℓ
       ≡ρ_k^2,

where we introduced ρ_k for the square root of the quantity in eq:rdef0.
Thus, ^*∈ B(_k,ρ_k).

Next, define

    σ_k=‖_k-^*‖,

so that ^*∈ B(_k,σ_k).

By the minimality
property of _k, we know that ∇ f(_k) is orthogonal to 
Tℳ_k,
which contains _k-_k, i.e.,

    ∇ f(_k)^T(_k-_k)=0.

Thus,

    ‖_k-_k‖   = ‖(_k-_k)+(_k-_k)‖
       = ‖(_k-_k)+∇ f(_k)/ℓ‖
       = √(‖_k-_k‖^2+‖∇ f(_k)/ℓ‖^2)
       ≥‖∇ f(_k)‖/ℓ,

so define

    δ_k =  ‖∇ f(_k)‖/ℓ,

to conclude that ‖_k-_k‖≥δ_k.  
We have defined δ_k,ρ_k,σ_k as in
Lemma <ref>.
We need to confirm the inequality
ρ_k+σ_k≥δ_k:

    ρ_k+σ_k   ≥‖_k-^*‖ + ‖_k-^*‖
       ≥‖_k-_k‖
       ≥δ_k.

The other inequality is derived as follows. First,
ρ_k≤δ_k since δ_k^2 is the first term in eq:rdef0.
Also, σ_k≤δ_k since

    σ_k^2    = ‖_k-^*‖^2 
       ≤‖_k-^*‖^2       
       ≤2(f(_k)-f(^*))/ℓ      
       ≤δ_k^2       .

Thus, δ_k≥max(ρ_k,σ_k) so δ_k^2≥ |ρ_k^2-σ_k^2|.

Therefore, we can conclude from  Lemma <ref> that
there exists a _k^*∈{_k,_k} (and hence in 
ℳ_k+1)
such that 

    ‖_k^*-^*‖≤ξ_k^*,

where ξ_k^* is defined by eq:xistardef for ρ_k,σ_k,δ_k 
given by eq:rdef,
eq:sdef and eq:ddef respectively.  After some simplification
and cancellation of eq:xistardef, one arrives at:

    (ξ_k^*)^2=‖_k-^*‖^2 - (f(_k)-f(^*)+‖_k-^*‖^2·ℓ/2/‖∇ f(_k)‖)^2.

Since _k+1 is the optimizer of eq:ia2.yupd, 
_k+1 is
at least as close to ^* as _k^*, and hence, 

    ‖_k+1-^*‖^2≤‖_k-^*‖^2 - (f(_k)-f(^*)+‖_k-^*‖^2·ℓ/2/‖∇ f(_k)‖)^2.

Adding this inequality to eq:fdesc yields:

    Ψ_k+1   =‖_k+1-^*‖^2 + 2(f(_k+1)-f(^*))/ℓ
       ≤‖_k-^*‖^2 - (f(_k)-f(^*)+‖_k-^*‖^2·ℓ/2/‖∇ f(_k)‖)^2 + 2(f(_k)-f(^*))/ℓ-
    ‖∇ f(_k)‖^2/Lℓ
       ≤‖_k-^*‖^2 -2[f(_k)-f(^*)+‖_k-^*‖^2·ℓ/2]/√(Lℓ)
    + 2(f(_k)-f(^*))/ℓ
       =
    [‖_k-^*‖^2+2(f(_k)-f(^*))/ℓ]
    ·(1-√(ℓ/L)) 
       =
    Ψ_k
    ·(1-√(ℓ/L)).

The third line was obtained by applying the inequality
 a^2+b^2≥ 2ab to the second and fourth terms of the second line.




§ ANALYSIS OF THE GEOMETRIC DESCENT ALGORITHM




In this section we present the
geometric descent (GD) algorithm due to
<cit.> and
an analysis of it. 
Our analysis varies slightly from the proof
due to <cit.>; in their proof, the
potential involves the term 2(f(_k)-f(^*))/ℓ
rather than 2(f(_k)-f(^*))/ℓ. 
The reason for
the change is to unify the analysis with the other algorithms
considered in order for the NCG construction in Section <ref>
to be applicable.


    _0:=
       _0:=_0 
       k:=1,2,…
       _k-1 : =_k-1-∇ f(_k-1)/L 
       _k-1 : =_k-1-∇ f(_k-1)/ℓ
       
       _k:=(1-λ_k)_k-1 + λ_k_k-1
       _k:={f():∈{_k-1,_k}}


Note: The operation in eq:gd.xupd is a line search and
requires an inner iteration to find the optimal  in the specified line.

This algorithm (which is one of several variants
of GD presented by the authors) is derived in <cit.>.
It can be regarded as
extracting the essential properties
of _k and _k used in the proof of Theorem <ref> to
obtain an implementable algorithm.
Indeed, the authors use that proof that we presented
in Section <ref> to analyze GD rather 
IA.  A more precise statement of the relationship between
GD and IA is provided in the next section.

Intuitively, the proof in the previous
section shows that _k need not be the minimizer
in eq:ia2.xupd; it suffices for _k to satisfy the two
properties eq:fdesc0 and eq:gradorth.  The Geometric
Descent algorithm satisfies these two properties with a “dogleg” step
in eq:gd.xupd
that combines a gradient step with a step toward _k.
Property eq:fdesc0 is satisfied because f(_k)≤ f(_k-1),
and property eq:gradorth is satisfied because of the minimality of
_k with respect to {_k-1,_k}.

The proof also shows that it suffices to take a _k+1 that
satisfies the inequality for _k^* in 
eq:zbound rather than solving eq:ia2.yupd.

We now turn to the computation of λ_k and the
associated issues with the radii ρ_k,σ_k.
Recall that ρ_k^ IA from eq:rdef0
and σ_k^ IA from eq:sdef both involve
^* and hence are unimplementable.  The difficulty with
σ_k is straightforward to resolve: define σ_k
to be an upper on ‖_k-^*‖ rather than its exact value,
and ensure inductively that σ_k+1 is an upper
bound on ‖_k+1-^*‖.

The difficulty with eq:rdef0 is resolved using 
offsets denoted by γ_k, 
a clever  device from <cit.>.  As in eq:rdef0 and eq:rdef,

    ‖_k-^*‖^2 
       ≤‖∇ f(_k)‖^2/ℓ^2 
    -2(f(_k)-f(^*))/ℓ
       ≡ρ_k^2,
       ≡ρ̃_k^2 - γ_k

where

    ρ̃_k    = ‖∇ f(_k)‖/ℓ, 
    γ_k    = 2(f(_k)-f(^*))/ℓ.

Let  σ_0,σ_1,…, be a sequence of positive scalars 
such that σ_k≥‖_k-^*‖ for all k=0,1,…,
and suppose that

    σ̃_k = (σ_k^2 + γ_k)^1/2.

Thus, we have the relationships:

    σ̃_k^2    = σ_k^2 + γ_k, 
    ρ̃_k^2    = ρ_k^2 + γ_k.

Note that ρ̃_k is easily computable on the kth iteration, while
σ̃_k can be updated recursively.  The rationale of these definitions is
as follows. From eq:lambdastar, one sees that
if σ^2 and ρ^2 are both incremented by the same constant
additive term γ_k, 
then λ^* is unaffected.  Also, it follows from this observation
and from eq:xistardef that if σ^2 and ρ^2 are both incremented
by γ_k, then (ξ^*)^2 is also incremented by γ_k.  Thus,
the GD algorithm works throughout with radii whose squares are incremented by γ_k.
This increment γ_k changes from one iteration to the next and hence must
be adjusted at the start of each iteration (see eq:sigmadef below).

In more detail, the sequence of computations is as follows.
We initialize the algorithm by:

    σ_0:=√(2)‖∇ f(_0)‖/ℓ.

This initialization is carried out prior to the main loop of GD.
The rationale for this formula is provided in the proof of
eq:sigma-ub below.

Assuming inductively that σ̃_k-1 is already known, 
compute as follows:

    ρ̃_k-1 := ‖∇ f(_k-1)‖/ℓ
       
       δ_k-1 := ‖_k-1-_k-1‖, 
       λ_k := 
    δ_k-1^2+ρ̃_k-1^2-σ̃_k-1^2/2δ_k-1^2, 
          
       ξ̃^*_k :=
    1/2√(2ρ̃_k-1^2+2σ̃_k-1^2-δ_k-1^2-
    (ρ̃_k-1^2-σ̃_k-1^2)^2/δ_k-1^2),
          
       σ̃_k := √((ξ̃^*_k)^2-γ_k-1+γ_k)
       
       λ_k := 0 
       σ̃_k := √(ρ̃_k-1^2 - γ_k-1+γ_k).

Although computation of γ_k alone requires prior knowledge of ^*, 
the difference γ_k-γ_k-1 appearing
in eq:sigmadef and eq:sigmadef2 does not, as is
evident from eq:gammakdef.  This is the motivation for
using σ̃_k and ρ̃_k in the computation instead
of σ_k,ρ_k.
In the theorems below
it is confirmed that the square
roots in eq:sigmadef and eq:sigmadef2 take nonnegative
arguments.

The convergence of the GD algorithm is proved via
two theorems, which are both variants of
theorems due to <cit.>.

Before stating and proving the two theorems, we establish two
inequalities.  
As noted earlier, eq:gradorth holds for GD, and
hence so does eq:ykoobx.
Then it follows from eq:ykoobx
combined with eq:gd.rpdef, eq:gd.deltadef that

    ρ̃_k-1= ‖∇ f(_k-1)‖/ℓ≤δ_k-1.


Next, regarding γ_k-γ_k-1 appearing in
eq:sigmadef and eq:sigmadef2, observe

    γ_k-γ_k-1    = 2(f(_k)-f(_k-1))/ℓ
       ≤2(f(_k-1)-f(_k-1))/ℓ
       ≤ -‖∇ f(_k-1)‖^2/Lℓ

where the second line follows because f(_k)≤ f(_k-1) by
eq:gd.xupd while the third follows by eq:fdesc0.


For all k=1,…,

    σ̃_k^2 ≥‖_k-^*‖^2 + 2(f(_k)-f(^*))/ℓ.




Note that the statement of the theorem may be equivalently written,

    σ̃_k^2 ≥‖_k-^*‖^2 + γ_k.


The proof is by induction.  
The base case is that σ̃_0^2≥‖_0-^*‖^2+γ_0.
Both terms of the right-hand side may be bounded by noting that strong convexity applied
to the two points _0,^* and rearranged may be written:

    -2/ℓ∇ f(_0)^T(^*-_0)≥‖_0-^*‖^2 +γ_0.

Again by strong convexity, ‖_0-^*‖≤‖∇ f(_0)‖/ℓ, so
we can apply this inequality and the Cauchy-Schwarz inequality on the left-hand
side to obtain

    2‖∇ f(_0)‖^2/ℓ^2≥‖_0-^*‖^2 +γ_0.

Thus, the bound
σ̃_0^2≥‖_0-^*‖^2+γ_0 is assured
by eq:sigma0.

For the induction case, assume k≥ 1 and
the induction hypothesis

    σ̃_k-1^2 ≥‖_k-1-^*‖^2 + γ_k-1.

There are two possibilities depending on the “if”-statement 
eq:gd.ifstmt.  First, suppose the condition
of eq:gd.ifstmt holds, which may
be rewritten as
σ̃_k-1^2-ρ̃_k-1^2 ≤ρ̃_k-1^2.
By eq:gd.deltarhobd, this implies  δ_k-1^2≥σ̃_k-1^2-ρ̃_k-1^2,
and we already know from eq:gd.deltarhobd that
δ_k-1^2≥ρ̃_k-1^2≥ρ̃_k-1^2-σ̃_k-1^2.  
The conclusion from all these inequalities is

    δ_k-1^2≥ρ̃_k-1^2≥ |ρ̃_k-1^2-σ̃_k-1^2|.

Referring now to Lemma <ref>, make the following identifications:

    = _k-1,
       = _k-1, 
    ρ^2     = ρ̃_k-1^2 - γ_k-1, 
    σ^2    = σ̃_k-1^2 - γ_k-1, 
    δ   = δ_k-1,

in order to apply the lemma.  The condition ρ+σ≥δ
follows immediately since we already have assumed by induction
that ^*∈ B(,σ) and it follows from eq:gd.rdef
that ^*∈ B(,ρ), thus implying that
B(,ρ)∩ B(,σ)∅.  The condition
δ^2≥ρ^2-σ^2 follows because δ_k-1≥ρ̃_k-1
as in eq:gd.chain1.  The condition 
δ^2≥σ^2-ρ^2 follows because
δ_k^2≥σ̃_k-1^2-ρ̃_k-1^2 (as established in
eq:gd.chain1).
Therefore, by the lemma,
if we define λ^* 

    λ^*=δ^2+ρ^2-σ^2/2δ^2
    =δ_k-1^2 +ρ̃_k-1^2 
    -σ̃_k-1^2/2δ_k-1^2,

i.e., the formula for λ_k in eq:gd.lambda1
(notice that the two terms γ_k-1 cancel),
and we define _k as in eq:gd.yupd, 
then

    ^*∈ B(_k,ξ)

where

    ξ^2    = 1/4(2ρ^2 + 2σ^2 - δ^2 -1/δ^2
    (ρ^2-σ^2)^2) 
        = 
    1/4(
    2ρ̃_k-1^2 
    +2σ̃_k-1^2 - 4γ_k-1 - δ_k-1^2 
    -1/δ_k-1^2(
    ρ̃_k-1^2 - σ̃_k-1^2)^2) 
       =
    (ξ̃_k^*)^2 - γ_k-1      
        =
    σ̃_k^2 - γ_k,       .

This establishes
the theorem in the first case.

If the condition in eq:gd.ifstmt fails, then
λ_k=0 as in eq:gd.lambda2, implying from
eq:gd.yupd that _k=_k-1.
Then ‖_k-^*‖^2≤ρ̃_k-1^2-γ_k-1
by 
eq:gd.rdef2, implying by eq:sigmadef2
that ‖_k-^*‖^2 ≤σ̃_k^2-γ_k,
thus establishing the theorem in the second case.



For each k=1,2,…,

    σ̃_k^2 ≤(1-√(ℓ/L))σ̃_k-1^2.




We again take two cases depending on whether
the condition in eq:gd.ifstmt holds.
If it holds, then 
eq:gd.chain1 in the preceding proof
holds.  
Observe that the function x↦ x+C/x for C>0 is unimodal on 
(0,∞)  with a minimizer at √(C), which means that if
all the other parameters are fixed, the maximizing choice for δ_k-1^2
in eq:gd.xik is δ_*^2=|ρ̃_k-1^2-σ̃_k-1^2|.
Therefore, by unimodality combined with the ordering
δ_k-1≥ρ̃_k-1≥δ_* 
(which is eq:gd.chain1),
the right-hand side of eq:gd.xik can only
increase if we replace δ_k-1 by ρ̃_k-1, thus obtaining,

    (ξ̃^*_k)^2
       =
    1/4(2ρ̃_k-1^2+2σ̃_k-1^2-δ_k-1^2-
    (ρ̃_k-1^2-σ̃_k-1^2)^2/δ_k-1^2) 
       ≤1/4(2ρ̃_k-1^2+2σ̃_k-1^2-ρ̃_k-1^2-
    (ρ̃_k-1^2-σ̃_k-1^2)^2/ρ̃_k-1^2) 
       = σ̃_k-1^2 -σ̃_k-1^4/4ρ̃_k-1^2.

Therefore,

    σ̃_k^2     = (ξ̃_k^*)^2 -γ_k-1+γ_k 
       ≤σ̃_k-1^2 -σ̃_k-1^4/4ρ̃_k-1^2
    -γ_k-1+γ_k 
        =
    σ̃_k-1^2 -σ̃_k-1^4/4‖∇ f(_k-1)‖^2/ℓ^2
    -γ_k-1+γ_k 
       ≤σ̃_k-1^2 -σ̃_k-1^4/4‖∇ f(_k-1)‖^2/ℓ^2
    -‖∇ f(_k-1)‖^2/Lℓ      
       ≤σ̃_k-1^2 -2·σ̃_k-1^2/2‖∇ f(_k-1)‖/ℓ·‖∇ f(_k-1)‖/√(Lℓ)      
        =
    σ̃_k-1^2 (1-√(ℓ/L)).


In the other case, σ̃_k-1^2/2≥ρ̃_k-1^2
so we obtain

    σ̃_k^2     = ρ̃_k-1^2 -γ_k-1+γ_k
          
       = ‖∇ f(_k-1)‖^2/ℓ^2  -γ_k-1+γ_k
       ≤‖∇ f(_k-1)‖^2/ℓ^2-
    ‖∇ f(_k-1)‖^2/Lℓ      
        = 
    ‖∇ f(_k-1)‖^2/ℓ^2(1-ℓ/L) 
        = 
    ρ̃_k-1^2(1-ℓ/L) 
       ≤σ̃_k-1^2(1-ℓ/L)/2

It is a simple matter to confirm that (1-)/2≤ (1-√())
for any ∈[0,1], thus establishing the theorem in this case.




§ RELATIONSHIP BETWEEN GD AND IA



We already observed that GD and IA both work in the same
affine subspace ℳ_k on each iteration.  In this section,
we develop further insight into their connection.
On each step IA solves two optimization
problems exactly, eq:ia2.xupd and eq:ia2.yupd,
to obtain _k and _k.  We argue
that GD computes optimal
solutions to these two problems not for the actual objective function f
but for some other objective function that agrees
with its partial information 
about f.
This is stated precisely in the following two theorems, the
first about _k and the second about _k.

Regarding the computation of _k, let us recall that
GD takes a steepest descent step from _k-1 in
line eq:gd.steep followed by a line search in
eq:gd.xupd.  The simplest line-search method for
minimizing a convex function
is bisection based on the sign of the directional derivative
of the objective along the line.  Evaluation of signs of derivatives
does not give any information other than the location of the minimizer
since univariate convex functions are unimodal.  (In other words,
bisection to find the minimizer of two convex functions f_1,f_2
with a common minimizer will evaluate the same sequence
of points for f_1 as for  f_2.)    Therefore, the information
about f used in GD to find _k
is the value of ∇ f(_k-1), the line
{_k-1,_k-1}, and the location of the minimizer on that
line.  The following theorem says that GD chooses the optimal
_k∈ℳ_k
given this partial information about the objective function
f.  (In contrast, IA chooses the optimal _k for the
true objective f.)


  Given scalar L>0, a point _k-1∈^n, n≥ 2,
  a nonzero vector ∈^n, define
  _k-1=_k-1-/L.  Also, assume we are
  given
  a line Λ⊂^n containing
  _k-1 but not _k-1, and point _k∈Λ.
  Let ℳ_k={Λ,_k-1}.
  Assume that
  
    ^T(_k-_k-1)<0.

  (Assumption eq:angassum will be explained later.)
  Define
  ℱ to be:
  
    ℱ={f:^n→:    , 
       {f():∈Λ}=_k, 
       , 
       ∇ f(_k-1)=}.

  Define q()=sup{f()-f(_k-1):f∈ℱ}.
  Then 
  
    q(_k)=min{q():∈ℳ_k}.


The class ℱ is meant to capture the set of all functions whose
partial information known to the GD algorithm agrees with the partial information
of the actual objective function, that is, ∇ f(_k-1)= and
_k is the minimizer of f on the line Λ.  The conclusion
of the theorem is that for the worst f in this class (the sup appearing
in the definition of q), _k chosen by GD is optimal over ℳ_k.
The hypotheses on the given data correspond to the induction hypotheses of GD
except for eq:angassum, which we discuss later.



  Since the algorithm is invariant under translation and rotation
  of space, without loss of generality we can 
  transform coordinates to identify ℳ_k with
  the (x_1,x_2)-plane. 
  Assume that the transformation places
  the point _k-1 at the origin (0,0),
  the line Λ on the x_1-axis, the point
  _k at (s,0).
  Denote the point _k-1 with (a,b) in this plane.
  Because of the identification of (0,0) with _k-1,
  we know that for f∈ℱ, =∇ f(_k-1)=∇ f(a,b)=(aL,bL) 
  in this rotated and translated coordinate system.  
  Since the length  is unchanged by translation and rotation,
  then (a,b) must satisfy the restriction that ((aL)^2+(bL)^2)^1/2 is the
  original length ‖‖.  Observe that
  f(0,0)-f(a,b)≤ -‖‖^2/(2L)=-(L/2)(a^2+b^2)
  by eq:fdesc0.
  For the remainder of the proof, let k_0=-(L/2)(a^2+b^2).
  Thus, f(s,0)-f(a,b)≤ k_0
  since (s,0) is the minimizer of f∈ℱ over Λ
  (hence achieves a value
  less than f(0,0)).  This inequality holds for all f∈ℱ,
  thus showing that q(_k)≤ k_0.

  In these transformed coordinates, assumption eq:angassum
  is rewritten as the inequality as<0.
  For the remainder of this proof, assume a>0 hence s<0; the
  other case is obtained by reflection of the x_1-coordinate.

  Consider the function,
  
    f̂_0(x,y) = {[ Lx^2/2 + Ly^2/2,            x≥ 0,;          Ly^2/2,            x≤ 0, ].

  The level curves of this function are semicircles in the right half-plane
  and parallel rays in the left half-plane.  It satisfies all the conditions
  for membership in ℱ except that (s,0) is not the unique minimizer
  over Λ;
  all points of the ray {(x,0):x≤ 0} are minimizers. Therefore,
  we perturb this function slightly.


  Fix a δ>0 small.
  Define scalars m=-δ (s+a)/(L-2δ) and p=δ m - δ s.
  
  
  
  
  
  
  
  Note that
  
  m,p
  
  tend to 0 as δ→ 0, so assume that δ is
  sufficiently small that 2δ<L, |m|<|s|, |m|<|a|.
  Consider the following function f̂_δ:
  
    f̂_δ(x,y) = {[ q_1(x-m)^2/2 + p(x-m)+Ly^2/2,                         x≥ m,;  q_2(x-m)^2/2+p(x-m) +Ly^2/2,                         x≤ m, ].

  where q_1=L-δ, q_2=δ.
  
  It is straightforward to check that f̂_δ∈ℱ.
  In particular, ∂f̂_δ(s,0)/∂ x=δ(s-m)+p=0 since
  p=δ(m-s).
  The choice of m ensures that ∂f̂_δ(a,b)/∂ x=aL.
  
  
  Observe that f̂_δ(s,0)-f̂_δ(a,b)→ k_0 as δ→ 0,
  and (s,0) is the minimizer of f_δ.  Therefore,  for all 
  ∈ℳ_k, f_δ()-f_δ(_k-1) is bounded
  below by k_0 plus a residual that tends to 0 as δ→ 0.
  which means
  that for all ∈ℳ_k, q()≥ k_0 (since f_δ is a candidate
  for the supremum in the definition of q(·)).
  Since we already established that q(_k)≤ k_0, this proves the theorem.



Let us know examine assumption eq:angassum, which
in transformed coordinates is written as<0.  If
as>0, then the minimizer over ℳ_k
cannot lie on the x-axis, i.e.,
ℱ is empty as the following argument shows.
(The case of as=0 needs separate
treatment, which we omit).  Consider an arbitrary f∈ℱ.
Select t to solve the following equation:

    ((s+at,bt) - (a,b))^T(a,b)=0;

one easily determines that

    t=a^2+b^2-as/a^2+b^2.

It is impossible that a^2+b^2-as<0, or equivalently, that t<0 because then

    f(s,0)   ≥ f(a,b)+∇ f(a,b)^T((s,0)-(a,b))        
        = f(a,b)+(aL,bL)^T((s,0)-(a,b)) 
       = f(a,b)+L(as-a^2-b^2) 
        > f(a,b),

contradicting the minimality of (s,0).

But there is also a contradiction when t≥ 0.
Consider the following chain of inequalities:

    f(s,0)+Lt^2‖(a,b)‖^2/2   ≥ f(s+at,bt)       
             
       ≥ f(a,b)+∇ f(a,b)^T((s+at,bt)-(a,b))
             
       =f(a,b)       
       ≥ f(0,0)+L/2‖(a,b)‖^2       
       ≥ f(s,0)+L/2‖(a,b)‖^2       .

This chain of inequalities starts and ends at the same quantity except for the
presence of t^2; thus, the inequalities can hold only if t≥ 1.  But
if as>0 and a^2+b^2-as≥ 0,
then it follows from eq:tformula that 0≤ t <1.  So we conclude
that as>0 contradicts the optimality of (s,0).
More comments on this matter appear at the end of this section.



The next theorem covers _k and requires a different construction.
Its format and interpretation are analogous to the previous
theorem.

  Suppose the following data is given: two points _k-1,_k-1 in
  ^n, n≥ 3, 
  a nonzero vector ∈^n, four positive scalars ρ,σ,δ,ℓ
  that satisfy the following conditions.  For
  the remainder of the theorem, let _k-1=_k-1-/ℓ.
  The conditions are:
  
  
  * δ=‖_k-1-_k-1‖,
  
  * ρ<‖‖/ℓ
  
  * ρ+σ≥δ,
  
  * δ≥√(|ρ^2-σ^2|), and
  
  * 
    (_k-1-_k-1)^T > ρ^2ℓ/λ^*,

    where λ^* is given by eq:lambdastar.
    (Assumption eq:yxbbcond will be explained  later.)
  
  Let ℱ be the following set of functions:
  
    ℱ={f:^n→:
           , 
       ∇ f(_k-1)=, 
       ‖‖^2/ℓ^2-2(f(_k-1)-min{f})/ℓ=ρ^2,
        
       {f}∈ B(_k-1,σ)}.

  Here, min{f} is shorthand for min{f():∈^n} and
  similarly for {f}.
  Let ℳ_k denote {_k-1,_k-1,_k-1}.
  Define
  
    q()=sup{‖-{f}‖: f∈ℱ}.

  Let _k be the point computed by the GD algorithm for this data
  using eq:gd.yupd.
  Then
  
    _k={q():∈ℳ_k}.




  By translating and rigidly rotating space, we can identify ℳ_k with
  the (x_1,x_2) plane.  
  After this coordinate transformation,
  we may assume _k-1 coincides with the origin (0,0).
  For the remainder of the discussion, remaining coordinates x_3,x_4,…,x_n  are 
  not written and are assumed to be 0's.
  Let us write _k-1=(a,b).
  Since _k-1 is at the origin, this implies =(aℓ,bℓ).
  Choose the transformation to make
  _k-1 lie on the positive x_1-axis at position (δ,0).
  (Recall that δ=‖_k-1-_k-1‖.)
  
  Let λ^* be as in eq:lambdastar. As in eq:gd.yupd,
  let  _k=(1-λ^*)_k-1+λ^*_k-1, which in
  this coordinate system is 
  _k=(δλ^*, 0).
  
  Let f be an arbitrary member of ℱ. The hypotheses
  on ℱ imply that 
  {f}∈ B(_k-1,ρ)∩ B(_k-1,σ), and therefore,
  by Lemmas <ref> 
  and <ref>, ‖{f}-_k‖≤ξ^*,
  where ξ^* is given by eq:xistardef.  This shows that
  q(_k)≤ξ^*.

  The remainder of the proof shows that
  q()>ξ^* for ∈ℳ_k-{_k}, which will
  establish eq:yksolve.   Let
  κ stand for either +1 or -1.  Embed
  ℳ_k in one higher dimension and write 3-tuples of
  coordinates (so that _k-1=(a,b,0), 
  =(aℓ,bℓ,0), _k-1=(δ,0,0) and
  so forth).  Define
  
    ^*=([ δλ^*;    0; κξ^* ]).

  Define 
  
    f()=(ℓ/2)‖-^*‖^2 +ℓ|(^*)^T(-^*)|.

  First, observe the obvious consequences of this
  formula that ^*={f} and that f(^*)=0.
  We claim that f∈ℱ.  The fact that f is 
  ℓ-strongly convex follows from the presence of the
  first term.  The second term is convex but not strongly convex.

  In order to establish the remaining conditions for
  membership in ℱ, we
  first determine which branch of the absolute value
  holds when evaluating f(_k-1); in particular,
  we establish the inequality that (^*)^T(_k-1-^*)> 0,
  i.e., that (^*)^T_k-1>‖^*‖^2.  The left-hand side
  evaluates to aδλ^*, while the right-hand side
  evaluates to (δλ^*)^2+(ξ^*)^2
  which simplifies to ρ^2 according to
  eq:lambdastar and eq:xistardef.  
  Thus, we must establish aδλ^*>ρ^2; this follows
  from eq:yxbbcond which states that
  δ aℓ>ρ^2ℓ/λ^* in the transformed coordinates.

  This inequality implies that in the neighborhood of
  _k-1, the absolute value sign appearing in the definition
  of f may be dropped.
  We now establish
  the remaining conditions for membership of f in ℱ.
  We have ∇ f(_k-1)=ℓ(_k-1-^*)+ℓ^*=ℓ_k-1=,
  so the second condition is established. 
  For the third condition, recalling that (^*)^T_k-1=aδλ^*
  while ‖^*‖^2=ρ^2, we compute,
  
    f(_k-1)    =(ℓ/2)((a-δλ^*)^2+b^2+(ξ^*)^2) + ℓ(aδλ^*-ρ^2) 
       =(ℓ/2)(a^2 + b^2 -ρ^2),

  where, to obtain the second line, we combined like terms and used
  the already established identity  (δλ^*)^2+(ξ^*)^2=ρ^2.
  Therefore,
  
    ‖‖^2/ℓ^2 - 2/ℓ(f(_k-1)-f(^*))
           = (a^2+b^2) - (a^2+b^2-ρ^2) 
       =ρ^2,

  thus establishing eq:rhocond2.  For the last condition,
  
    ‖_k-1-^*‖^2    = (δ-δλ^*)^2+(ξ^*)^2 
       = (σ^*)^2.


  Thus, membership of f∈ℱ is established.  Note that
  the two choices for κ, namely, ± 1 
  lead to two different optimizers ^* in eq:xoptdef1
  whose distance apart is
  2ξ^*.  Therefore, the midpoint of the two optimizers is the
  only point in ^n whose distance from both of them is
  bounded above by ξ^*.  This midpoint is exactly _k.  This
  proves that q(_k)≥ξ^* (we already showed that q(_k)≤ξ^*)),
  and also that for any ∈ℳ_k-{_k}, q()>ξ^*.
  This establishes eq:yksolve.
  

We now discuss the conditions imposed on the given data.
The condition δ=‖_k-1-_k-1‖ is
simply part of the construction, and the conditions
ρ<‖‖/ℓ and
ρ+σ≥δ ensure that ℱ is
nonempty.  The condition
δ≥√(|ρ^2-σ^2|) is the main case of
the two cases arising in the main loop of GD.
In particular, the condition eq:gd.ifstmt is used to establish
eq:gd.chain1.  We omit
the treatment of the minor case.

Finally, we discuss eq:yxbbcond.  If a strengthened version of this
condition fails to hold,  then
ℱ is empty.
In particular, strong convexity implies that the optimizer lies in a ball
of radius ‖‖/ℓ about _k-1, i.e.,

    ‖_k-1-^*‖^2 ≤‖‖^2/ℓ^2.

In the
transformed coordinates, this is written

    (a-δλ^*)^2+b^2+(ξ^*)^2≤ a^2+b^2,

which simplifies to

    2aδλ^*≥ρ^2,

while eq:yxbbcond is written as aδλ^*>ρ^2 in the
transformed coordinates, so the
two bounds differ by a factor of 2.  In the case that
ρ^2∈ [aδλ^*,2aδλ^*], the two possibilities are (1)
ℱ=∅, or (2) there is an f∈ℱ, but
the construction of f used in the proof
does not work.  We do not know
which possibility is correct.

We now summarize the results in this section with some observations.
We have shown that the computations of _k,_k in the GD algorithm
solve minimization problems akin to those in IA except that
the minimizer pertains to the worst case f that agrees with
the partial information that GD has about f rather than the true f.

In both cases, the theorems had some apparently extraneous assumptions,
namely eq:angassum in the first theorem, and eq:yxbbcond in the second.
However, in both these cases, the extraneous assumptions indicate
that GD may not be using the
information about f entirely.  In the case of the first theorem,
if ^T(_k-_k-1)>0
then the minimizer of f over ℳ_k cannot be on the
line searched by GD, yet GD does not use this information.
In the case of the second theorem, GD does not
use the fact that ^*∈  B(_k-1,‖∇ f(_k-1)‖/ℓ).

The main point of this section is to clarify the relationship between IA
and GD, but the preceding paragraph reveals a second point.  By constructing
these example functions to show that the GD algorithm is the best possible
in some cases but not others
given the limited information that it uses,
we also show that it may be possible to improve
on GD by making better use of the information
(such as the third ball mentioned in the
previous paragraph) in certain cases.
We do not pursue this idea here, but
see, e.g., <cit.> for results in this direction.



§ COMPUTABILITY OF THE POTENTIAL



A point to make about the GD algorithm and its
analysis is that the potential σ̃_k^2 is
computable on every iteration of the algorithm, i.e.,
it does not require prior knowledge of ^*.  (It does,
however, require prior knowledge of ℓ.)
This fact
can also be deduced from the original
presentation of <cit.>, although the computability
is not further used therein.

In this section we explain in more detail what we
mean by “computable”. In particular, a potential
σ̃_k that is an upper bound to eq:psidef
is computable if it has the
following properties.



  * It must be possible to compute the potential without
  prior knowledge of ^* or f(^*).
  This is the reason that we regard the potential
  Ψ_k itself defined in eq:psidef as
  noncomputable.

  * It must have an a posterior dependence on the actual convergence
  of the algorithm, in the sense that if
  _k,_k are very close to ^*, then it
  should be the case that σ̃_k is
  very close to zero.  In particular, this
  rules out using a completely a priori potential like
  
    σ̃_k := C‖∇ f(_0)‖^2(1-√(ℓ/L))^k

  for a fixed constant C.  Although this potential indeed
  is an upper bound on eq:psidef and is computable, it has no relationship
  to the current iterate and therefore has no algorithmic use.

  * The potential should decrease by the factor (1-√(ℓ/L))
  per iteration (or perhaps (1-·√(ℓ/L)), since a constant
  factor improvement may be possible).  Thus, although
  C‖∇ f(_k)‖^2 for a correctly chosen C is
  an upper bound on eq:psidef, it does not satisfy
  our requirement of steady decrease and in fact can be
  oscillatory.

  The significance of guaranteed decrease in the potential is twofold.
  First, the guaranteed decrease is useful for theoretical analysis
  to establish a linear convergence rate.  Indeed, it is used
  herein for this
  purpose to establish previously known convergence rates for AG
  and CG in a new manner.

  Second, steady decrease in the potential can be used in an
  algorithm to ensure that progress is being made.  To give
  one example not pursued herein, consider the problem of detecting stagnation
  due to imprecise arithmetic in linear conjugate gradient.
  Although
  ‖∇ f(_k)‖^2 is commonly used as a termination
  criterion for linear conjugate gradient, it is not suitable
  for use as a stagnation test because it can be highly oscillatory
  and therefore cannot be used to check whether a single
  iteration was successful.
  A steadily decreasing potential, however, could be used
  in a CG stagnation test.  (We have preliminary results
  on this matter that will be the subject of future work.)

  Our present algorithmic use of the potential as a measure of steady
  decrease is reported in Sections <ref>–<ref>.  
  We present an experiment to illustrate why
  ‖∇ f(_k)‖^2 is not a suitable substitute for
  σ̃_k in Section <ref>.


The potential for GD given by
eq:sigmadef and eq:sigmadef2.  As mentioned
earlier, it is a slight variant of the potential defined
by the authors of GD; theirs also has these properties.
It is somewhat surprising that the same potential also
applies to conjugate gradient and accelerated
gradient, as developed in the next few sections.
We do not know of any other computable potential with these
properties.




§ ANALYSIS OF LINEAR CONJUGATE GRADIENT



The linear conjugate gradient (CG) algorithm for 
minimizing f()=^TA/2 -^̱T, 
where
A is a symmetric positive definite matrix, 
is due to Hestenes and Stiefel
<cit.>
and is as follows.

    _0:=
       _̊0:=-̱A_0 
       _1:=_̊0 
        k:=1,2,…, 
       α_k:= _̊k-1^T_̊k-1/_k^TA_k
       _k := _k-1+α_k_k
       _̊k := _̊k-1-α_k A_k
       β_k+1:= _̊k^T_̊k/_̊k-1^T_̊k-1
       _k+1 :=  β_k+1_k+_̊k



We now show that linear CG exactly implements
Algorithm IA, and therefore also satisfies the bound of Theorem <ref>.
This is perhaps surprising because CG does not have prior
information about ^*.  The following key results about CG are from the
original paper:

(<cit.>)
Let 𝒱_k=_0+{_̊0,…,_̊k-1} in CG.
Then

(a) An equivalent
formula is 𝒱_k=_0+{_1,…,_k},

(b) _k is the minimizer of f() over 𝒱_k,

(c) _̊k=-∇ f(_k), and

(d) _k+τ_k_k is the minimizer of ‖-^*‖ over
𝒱_k, where

    τ_k = 2(f(_k)-f(^*))/‖_̊k-1‖^2.




Part (b) appears as Theorem 4:3 of <cit.>,
while parts (a) and  (c) are not stated explicitly.  All of (a)–(c)
are covered
by most textbook treatments of CG.  Part (d) appears as Theorem 6:5 and
is less well known.

The following theorem establishes the claim that Algorithm CG implements
IA.  The principal result is part (b).  The remaining parts
are necessary to support the induction proof.


Suppose IA and CG are applied to the same quadratic function with
the same _0.  Let the sequences of iterates be denoted
_k^ IA and _k^ CG respectively.
Then for each k=1,2,…,

(a) ℳ_k⊂𝒱_k,

(b) _k^ CG=_k^ IA, and

(c) _k=_k^ CG+τ_k_k.




For the k=1 case, observe that _1=-∇ f(_0) so
ℳ_1=𝒱_1. Since _1^ CG minimizes
f() over 𝒱_1 while _1^ IA minimizes
f() over ℳ_1, we conclude _1^IA=_1^ CG.
For (c), observe that _1 minimizes ‖-^*‖
over ℳ_1 by eq:ia2.yupd, while _1^ CG+τ_1_1
minimizes the same function over the same affine space, so (c)
is established.


Now assuming (a)–(c) hold for some k≥ 1, we establish them for k+1.
We will write _k for both _k^ CG and  _k^ IA
since these
are equal by induction.  For (a), we start with
ℳ_k+1=_k+{_k-_k,∇ f(_k)}.
We already know from (c) that _k-_k=τ_k_k∈ T𝒱_k
⊂ T𝒱_k+1.  Also,
∇ f(_k)=-_̊k=β_k+1_k-_k+1 (by eq:pupd), hence
∇ f(_k)∈ T𝒱_k+1
Thus, Tℳ_k+1⊂ T𝒱_k+1,
so showing ℳ_k+1⊂𝒱_k+1 is reduced to finding
a single common point, and we may take _k
to be this point.

For (b), _k+1^ CG minimizes f() over 𝒱_k+1
by Theorem <ref>(a).
We also know that _k+1^ CG∈ℳ_k+1
because 

    _k+1^ CG   =_k+α_k+1_k+1      
       =_k+α_k+1(_̊k+β_k+1_k)        
       =_k+α_k+1(-∇ f(_k) + β_k+1_k) 
          
       =_k+α_k+1(-∇ f(_k) + β_k+1(_k-_k)/τ_k)       
       ∈_k+{∇ f(_k), _k-_k}.

Since ℳ_k+1⊂𝒱_k+1 according to part (a), the
optimality of _k+1^ CG with respect to 𝒱_k+1 implies
that it is also optimal for f() with respect to ℳ_k+1,
hence _k+1^ IA=_k+1^ CG.  Thus, write _k+1 for 
both vectors for the remainder of the argument.

A similar argument establishes (c).  First, we observe
that _k+1+τ_k+1_k+1 lies in ℳ_k+1 
because _k+1+τ_k+1_k+1=_k+(α_k+1+τ_k+1)_k+1,
and the we can proceed as in the last paragraph except with 
(α_k+1+τ_k+1) playing the role of α_k+1.  Next,
_k+1+τ_k+1_k+1 minimizes ‖-^*‖ over 𝒱_k+1
by Theorem <ref>(d).  Thus, _k+1+τ_k+1_k+1 must
be the same as _k+1.  This concludes the induction.



The surprising aspect of this analysis is that CG exactly
identifies ℳ_k that appears in Algorithm IA despite
not ever computing _k^ IA.  The reason is that the line
{_k,_k^ IA} agrees with the
line  _k+{_k}, which is computed by CG.



§ A COMPUTABLE POTENTIAL FOR LINEAR CONJUGATE GRADIENT



The analysis in the preceding section shows that CG implements
the idealized algorithm.  However, the decrease in the potential
cannot be measured during the algorithm because eq:tau_k requires
prior knowledge of the optimizer.  In this section, we observe
that the GD 
potential can also be used for CG, yielding
a computable potential.
An application of
this potential will be presented in Section <ref>.   

We define an auxiliary sequence of vectors _k
using the formulas in GD.  This sequence is not the
true minimizer that occurs in IA
eq:ia2.yupd and in
Theorem <ref>, part (d). But nonetheless,
‖_k-^*‖ shrinks sufficiently fast
to establish the necessary decrease in the potential.

In particular, we exactly mimic the equations that
define the quantities
σ̃_k, ρ̃_k, λ_k, δ_k, _k
in GD, and modify only _k so that it is computed using
the CG algorithm instead of the GD algorithm.
The two same two theorems that held for GD also hold for CG:


For each k=0,1,2,…,

    ‖^*- _k‖^2 + 2(f(_k)-f(^*))/ℓ≤σ̃_k^2.





For each k=1,2,…,

    σ̃_k^2≤(1-√(ℓ/L))σ̃_k-1^2.



The following observations about _k
computed in CG show that the same proofs of the
previous theorems work for CG.

First, _k∈𝒱_k, whereas _̊k is orthogonal
to T𝒱_k (a well known property of CG), and thus 
∇ f(_k)^T(_k-_k)=0.
This is eq:gradorth, which was a necessary ingredient
in the proof of GD.

Second, eq:gammadiff still holds because 
the CG step from _k-1 to _k improves
f at least as much as the step from  _k-1
to _k-1, since _k lies in the Krylov space 
𝒱_k
where _k is optimal for f over this space.
 


§ ACCELERATED GRADIENT



The Accelerated Gradient (AG) algorithm of Nesterov is an even looser approximation
to IA than GD in the sense that there is no optimization subproblem per iteration;
instead, all step lengths are fixed.   
For this section, let us define

    κ = L/ℓ,

because this ratio, sometimes called the
condition number of f, is used often throughout the algorithm and
analysis.

The algorithm is as follows.  


    _0:=
       _0:=_0 
        k:=1,2,…, 
       _k := _k-1-∇ f(_k-1)/L 
       _k := _k + θ(_k-_k-1)
       θ=
           √(κ)-1/√(κ)+1.


Note that some versions of AG in the literature
vary the choice of θ
(e.g., see <cit.>).

For the purpose of analysis, let us define the following auxiliary sequences
of vectors and scalars:

    _k    = _k-∇ f(_k)/ℓ,        (k=0,1,…)
    
    _0    =_0, 
    _k    = _k+τ(_k-_k-1),
              (k=1,2, …) 
    σ̃_0    = √(2)‖∇ f(_0)‖ /ℓ, 
    σ̃_k+1   = 
    [(1-κ^-1/2)σ̃_k^2+ 2(f(_k+1)-f(_k))/ℓ
         +‖∇ f(_k)‖^2/Lℓ-(κ^1/2-κ^-1/2)‖_k-_k‖^2]^1/2       (k=0,1,2,…)

where

    τ= √(κ)-1.


We prove two main results about these
scalars.  The first shows that
σ̃_k is decreasing at the appropriate
rate, while the second shows that it is
an upper bound on the distance to the optimizer.


For each k=0,1,2,…,

    σ̃_k+1^2 ≤ (1-κ^-1/2)σ̃_k^2.




By squaring both sides of eq:ag.tsigma, it is apparent that
eq:agthm1 reduces to showing:

    2(f(_k+1)-f(_k))/ℓ +
    ‖∇ f(_k)‖^2/Lℓ-(κ^1/2-κ^-1/2)‖_k-_k‖^2≤ 0.

Clearly it suffices to show

    2(f(_k+1)-f(_k))/ℓ +
    ‖∇ f(_k)‖^2/Lℓ≤ 0.

This follows immediately from eq:wupd, which implies
that f(_k+1)≤ f(_k)-‖∇ f(_k)‖^2/(2L).




For each k=0,1,2…,

    ‖_k-^*‖^2 +2(f(_k)-f(^*))/ℓ≤σ̃_k^2.





The proof of eq:ag.tsigmabd1 is by induction on k.  
We start by deriving some preliminary relationships.
It is clear from eq:ag.xupd and eq:ag.ydef that _k,_k,_k
are collinear and the tangent to their common line is _k-_k-1, hence
we easily obtain from these equations:

    _k   =_k+τ/θ(_k-_k)       
       =_k+(√(κ)+1)(_k-_k)       
       =(√(κ)+1)_k - √(κ)_k.



For the k=0 case,
eq:ag.tsigmabd1 follows from the initialization in
eq:ag.tsigma0 and strong convexity.

We now assume the result eq:ag.tsigmabd1 holds for
k and establish it for k=1.
The proof relies on Lemma <ref>, so first we must
argue that _k+1 lies on the
line segment between _k and _k.  This is the content
of the following derivation:

    _k+1   = _k+1+(√(κ)-1)(_k+1-_k)       
       =√(κ)_k+1-(√(κ)-1)_k 
       =√(κ)_k -(√(κ)-1)_k- √(κ)∇ f(_k)/L 
          
       =√(κ)_k -(√(κ)-1)_k- κ^-1/2∇ f(_k)/ℓ      
       =κ^-1/2_k+(1-κ^-1/2)((√(κ)+1)_k
    -√(κ)_k) 
         
    - κ^-1/2∇ f(_k)/ℓ
       =κ^-1/2_k+(1-κ^-1/2)_k- κ^-1/2∇ f(_k)/ℓ      
        = κ^-1/2_k + (1-κ^-1/2)_k.



We take
,  appearing in  Lemma <ref> to be
_k, _k respectively.
Next we need to define δ,ρ,σ to be used
in Lemma.
In the case of ρ and 
we copy the definitions used in the analysis of IA:

    ‖_k-^*‖^2    ≤‖∇ f(_k)‖^2/ℓ^2 - 2(f(_k)-f(^*))/ℓ      
       ≡ρ^2.


For σ, we use the induction hypothesis:

    ‖_k-^*‖^2   ≤σ̃_k^2 - 2 (f(_k)-f(^*))/ℓ
       ≡σ^2.



In the case of δ, we have:

    ‖_k-_k‖^2    = 
    ‖_k-∇ f(_k)/ℓ - 
    (√(κ)+1)_k+√(κ)_k)‖^2       
        = ‖√(κ)(_k-_k)-∇ f(_k)/ℓ‖^2 
       =κ‖_k-_k‖^2 + 2√(κ)(_k-_k)^T∇ f(_k)/ℓ
         
    + ‖∇ f(_k)‖^2/ℓ^2
       ≥κ‖_k-_k‖^2 + 2√(κ)(f(_k)-f(_k))/ℓ
         
    +√(κ)‖_k-_k‖^2
    + ‖∇ f(_k)‖^2/ℓ^2      
       =(κ + √(κ))
    ‖_k-_k‖^2 + 2√(κ)(f(_k)-f(_k))/ℓ
         
    + ‖∇ f(_k)‖^2/ℓ^2
       ≡δ^2.


From eq:ag.lambdadef, λ = 1-κ^-1/2
(and hence λ(1-λ)=κ^-1/2-κ^-1).
Finally, by Lemma <ref>,

    ‖_k+1-^*‖^2    ≤κ^-1/2ρ^2+(1-κ^-1/2)σ^2
    -(κ^-1/2-κ^-1)δ^2 
       =
    ‖∇ f(_k)‖^2/Lℓ-2(f(_k)-f(^*))/ℓ
    +(1-κ^-1/2)σ̃_k^2-(κ^1/2-κ^-1/2)‖_k-_k‖^2 
       =
    σ̃_k+1^2-2(f(_k+1)-f(^*))/ℓ,

thus completing the induction.  The second line was obtained 
by substituting eq:ag.rhodef, eq:ag.sigmadef, eq:ag.deltadef
in the first line
followed by cancellation of like terms.  The third line was obtained
from eq:ag.tsigma.




§ RELATIONSHIP BETWEEN IA AND AG



The relationship between the idealized algorithm and AG is weaker than
that between IA and either GD or CG because AG does not solve any
optimization subproblems and instead takes fixed stepsizes.
Thus, at best it is an approximation to IA.  Furthermore,
the computations of _k and _k are more closely tied together,
making it unclear whether any kind of induction hypothesis can be
applied to either in isolation. For these reasons, we propose the following
theorem characterizing the AG–IA relationship.


  Suppose one is given points _k-1,_k-1∈^n,
  a nonzero vector ∈^n, n≥ 3, and
  scalars ℓ,L such that L>ℓ>0.
  For the remainder of this discussion, define κ=L/ℓ
  and
  
    _k-1=√(κ)/√(κ)+1_k-1+1/√(κ+1)_k-1,

  as in eq:wxy0.
  Assume further (these assumption will be explained later) that
  
    ‖_k-1-_k-1‖≤‖‖/L,

  and
  
    (_k-1-_k-1)^T=0.

  Let
  
    ℱ = {f:^n→:    , 
       , 
       ∇ f(_k-1)=}.

  Define
  
    q(,)=sup_f∈ℱ(2/ℓ)(f()-min{f})+‖-{f}‖^2/(2/ℓ)(f(_k-1)-min{f})+‖_k-1-{f}‖^2.

  Here min{f} and {f} are short-hand for
  min{f():∈^n} and {f():∈^n} respectively.
  Define ℳ_k=_k-1+{,_k-1-_k-1}.
  Let _k,_k be the vectors computed by the AG
  algorithm for this data (which lie in ℳ_k).   Then
  
    q(_k,_k)≤ 1-κ^-1/2,

  and
  
    min{q(,):,∈ℳ_k}≥ 1-κ^-1/2-O(κ^-1).


Note that eq:agbound1 and eq:agbound2 imply that
the choice of new iterate (_k,_k) made by AG is optimal up to a lower
order remainder term given the partial information used by AG.  In
contrast, IA is optimal (separately) for both terms in the numerator
of eq:qdef for the specific f and with no remainder term.


The proof of eq:agbound1 appears in <cit.>
(see (21), (22) and (43) in that paper, which use different notation
for AG).  Therefore, this proof establishes only eq:agbound2,
which involves constructing
a certain f∈ℱ to attain this bound.
It suffices to prove the result for the n=3
case, since we can extend f:^3→ to higher dimensions by adding
terms ℓ x_4^2+⋯+ℓ x_n^2 and appending 0's to _k-1,
_k-1 and .

The theorem is invariant under rigid motions of ^3, so we
can place _k-1 at an arbitrary point.  In addition,
we can rotate the two
vectors  and _k-1-_k-1 to arbitrary positions as long
as their lengths and their orthogonality are preserved. Starting with
, let g_0 denote ‖‖. 
Rotate  to

    =γ(
    [ -Lκ^-3/4(=-ℓκ^1/4);                  0;           ℓ/1-κ^-1;                    ]),

where γ is chosen to preserve the length of , in other words,

    γ=(g_0/ℓ) · (κ^1/2+(1-κ^-1)^-2)^-1/2.

Translate _k-1 as follows:

    _k-1=
    γ(
    [  -κ^-3/4;        0; 1/1-κ^-1 ]).

Finally, we rotate _k-1-_k-1 as follows:

    _k-1-_k-1=
    γ(
    [ κ^-5/4v^2(1-κ^-1)/(1-κ^-1)^2+κ^-1/2;                             vκ^-3/4; v^2κ^-1-κ^-3/2v^2/(1-κ^-1)^2+κ^-1/2 ])

where v is a scalar parameter that controls the length
of _k-1-_k-1.  In more detail, observe that
‖_k-1-_k-1‖=γκ^-3/4|v|+O(κ^-1),
while ‖‖/L=γκ^-3/4+O(κ^-1).  Therefore,
to assure eq:assum5, we restrict |v|≤ 1.
It is also straightforward
to check that eq:assum6 is satisfied.
With these two definitions in hand, we can now write:

    _k-1=_k-1+(_k-1-_k-1)=
    γ(
    [   -κ^-3/4+κ^-5/4v^2(1-κ^-1)/(1-κ^-1)^2+κ^-1/2;                                       vκ^-3/4; 1/1-κ^-1 +v^2κ^-1-κ^-3/2v^2/(1-κ^-1)^2+κ^-1/2 ]),

and

    _k-1 = _k-1-√(κ)(_k-1-_k-1)=
    γ(
    [   -κ^-3/4-κ^-3/4v^2(1-κ^-1)/(1-κ^-1)^2+κ^-1/2;                                      -vκ^-1/4; 1/1-κ^-1 -v^2κ^-1/2+κ^-1v^2/(1-κ^-1)^2+κ^-1/2 ]).


Next, we define:

    f()=Lx_1^2+√(Lℓ)x_2^2+ℓ x_3^2/2.

It is straightforward to verify that ∇ f is L-Lipschitz,
that f is ℓ-strongly convex, and that
∇ f(_k-1)=; thus f∈ℱ.  Also,
it is obvious that {f}= and min{f}=0.

We evaluate the two terms in the denominator of eq:qdef:

    f(_k-1) - min{f}   = f(_k-1) 
       =γ^2/2(L(κ^-3/4+O(κ^-5/4))^2+
    √(Lℓ)(vκ^-3/4)^2+
    ℓ(1+O(κ^-1))^2) 
       = γ^2ℓ/2(1+κ^-1/2+O(κ^-1)).

Also,

    ‖_k-1-{f}‖^2    = ‖_k-1‖^2 
       =γ^2(O(κ^-3/4)^2+
      v^2κ^-1/2+
      (1-v^2κ^-1/2+O(κ^-1))^2) 
       =γ^2(1-v^2κ^-1/2+O(κ^-1)).

Thus, the denominator of eq:qdef simplifies to
γ^2(2+(1-v^2)κ^-1/2+O(κ^-1)).

Next, we need to solve two constrained
quadratic optimization problems to obtain a lower bound on
the numerator of eq:qdef.  The constraint is
∈ℳ_k for the first and
∈ℳ_k for the second.
Imposing the constraint is simpler if it is written
as an inhomogenous linear equation;
one checks that

    ℳ_k   =_k-1+{_k-1-_k-1,}
       ={∈^3:κ^-1/4x_1/1-κ^-1-vκ^-1/4x_2+x_3=γ},

by substituting _k-1, _k-1-_k-1 and  into the
left-hand side of given equation and confirming that the values are γ, 0, 0
respectively.
It is also straightforward to check using a Lagrange-multiplier
argument that for any positive A,B,C and any (a,b,c) (0,0,0),

    min{Ax_1^2+Bx_2^2+Cx_3^2:ax_1+bx_2+cx_3=γ}=γ^2/a^2/A+b^2/B+c^2/C.

Using this result for the first term of the numerator of eq:qdef
yields that for any ∈ℳ_k,

    f()-min{f}   ≥1/2·γ^2/κ^-1/2/L(1-κ^-1)^2 +
      v^2κ^-1/2/√(Ll) +
      1/ℓ
       = γ^2ℓ/2(1+O(κ^-1)).

As for the second term of the numerator, for any ∈ℳ_k,

    ‖-{f}‖^2    ≥γ^2/κ^-1/2/(1-κ^-1)^2 +
        v^2κ^-1/2 +
        1
       = γ^2(1-(1+v^2)κ^-1/2+O(κ^-1)).

Therefore, a lower bound on the numerator of eq:qdef
is γ^2(2-(1+v^2)κ^-1/2+O(κ^-1)).
Finally, for any ,∈ℳ_k,

    q(,)   ≥γ^2(2-(1+v^2)κ^-1/2+O(κ^-1))/γ^2(2+(1-v^2)κ^-1/2+O(κ^-1))
       =1-κ^-1/2+O(κ^-1).



We now turn to the assumptions of the preceding theorem.
An assumption like eq:assum5 is necessary because,
in the situation that ‖‖≪‖_k-1-_k-1‖,
strong convexity (see eq:rdef0)
implies that the true minimizer ^* is
close to _k-1-/ℓ, meaning that the update
to _k implicit in AG is suboptimal. Thus,
an assumption along the lines of
eq:assum5
is necessary to establish the optimality of AG.

Orthogonality assumption eq:assum6 appears to be unnecessary
and rather is a limitation of our construction, which uses
a quadratic function f.  Intuitively,
we need to construct a function
that varies more rapidly in one direction than another.  We
used a quadratic function, whose level curves have fixed orthogonal axes,
which creates a requirement of orthogonality in the two directions.
However, a more general convex function may have level curves
whose axes of elongation vary from one level curve to the next.

As in the concluding remarks
of Section <ref>,
the proof of the optimality of the algorithm combined with
consideration of the assumptions uncovers situations when the algorithm
may be making suboptimal choices.  In the case of AG, this occurs
on iterations when ‖∇ f(_k-1)‖ is unexpectedly small.




§ A HYBRID NONLINEAR CONJUGATE GRADIENT



In this section, we propose a hybrid nonlinear conjugate gradient algorithm with
a convergence guarantee for smooth, strongly convex functions which is
related to an algorithm from the PhD thesis of the first author
<cit.>.  Classical
nonlinear conjugate gradient (NCG) methods such as the Fletcher-Reeves and
Polak-Ribière methods are known to have poor worst-case performance
for this class of functions–worse even than steepest descent.  See
<cit.> for more information.  The method developed in
this section guarantees O(log(1/ϵ)√(L/ℓ)) convergence,
the best possible, and reduces to the optimal CG 
algorithm in the case of a quadratic function.

The algorithm proposed below uses classical nonlinear conjugate
gradient steps mixed with geometric descent steps.
The rationale for developing this algorithm is as follows.  Classical NCG, although
it has no global convergence bound even for strongly convex functions,
behaves well on “nearly quadratic” functions.  For typical objective
functions occurring in practice, nearly quadratic behavior is expected close
to the solution.  Therefore, a method that can switch between steps with
a guaranteed complexity and NCG steps has the possibility of outperforming
both methods.

A summary of the algorithm is as follows.  At the beginning of
iteration k, the algorithm has a quadruple
(_k-1,_k-1,_k-1,σ̃_k-1).  From this
quadruple, a step of nonlinear conjugate gradient can be
applied.  For the line search, the line-search function of α,
namely, f(_k-1+α_k-1), is approximated by a univariate
quadratic, whose quadratic coefficient is obtained by computing
_k-1^T∇^2 f(_k-1)_k-1 using reverse-mode automatic
differentiation.  This approximation is exact in the case that f
itself is a quadratic function, in which case the hybrid algorithm
reproduces the steps of linear conjugate gradient.

The hybrid algorithm then computes _k=_k-1+α_k-1 and
computes _k 
as in the GD algorithm.  It checks whether f has decreased and
whether σ̃_k^2≤ (1-√(ℓ/L))σ̃_k-1^2.
If so, the iteration is over, and the nonlinear CG step is accepted.
If not, then a GD step is taken instead.
The detailed specification of the algorithm is as follows.





    _0:=;  _0:=_0;  _0:=
       _-1:=;
    σ̃_0 := √(2)‖∇ f(_0)‖/ℓ
       k = 1,2,…, 
       _k-1:=∇ f(_k-1)
        (_k^ CG, _k) :=
     CGSTEP(_k-1,_k-1,_k-2,_k-1,k) 
        (_k,ξ̃^*_k) :=
     YCOMPUTE(_k-1,_k-1,_k-1,σ̃_k-1) 
       γ̂_k^ CG := 2(f(_k^ CG)-f(_k-1))/ℓ
       σ̃_k^ CG := 
    √((ξ̃^*_k)^2+
    γ̂_k^ CG)
       γ̂_k^ CG≤ 0 
    (σ̃_k^ CG)^2≤(1-√(ℓ/L))σ̃_k-1^2 
    
       _k:=_k-1^ CG
       σ̃_k:=σ̃_k^ CG
       
       _k-1:=_k-1-_k-1/L 
       _k:={f():∈{_k-1,_k}}
       γ̂_k := 2(f(_k)-f(_k-1))/ℓ
       σ̃_k := 
    √((ξ̃^*_k)^2+
    γ̂_k)
       _k := _k-_k-1



    k == 1 
       _k := -_k-1
       
        := _k-1-_k-2
       β_k := 
    1/^T_k-1(
     - 2_k-1‖‖^2/^T_k-1)^T_k-1
       _k:= β_k_k-1-_k-1
       
       α_k := -_k^T_k-1/_k^T∇^2 f(_k-1)_k
       _k := _k-1+α_k_k 
        (_k,_k)



    _k-1 := _k-1-_k-1/ℓ
       ρ̃_k-1 := ‖_k-1‖ /ℓ
       σ̃_k-1^2 ≤ 2ρ̃_k-1^2 
       δ_k-1 := ‖_k-1-_k-1‖
       δ_k-1 > ρ̃_k-1ρ̃_k-1>|ρ̃_k-1^2-σ̃_k-1^2|^1/2
       λ_k := 
    δ_k-1^2+ρ̃_k-1^2-σ̃_k-1^2/2δ_k-1^2
       ξ̃^*_k :=
    1/2√(2ρ̃_k-1^2+2σ̃_k-1^2-δ_k-1^2-
    (ρ̃_k-1^2-σ̃_k-1^2)^2/δ_k-1^2)
       
       λ_k := 1 
       ξ̃_k^*:= σ̃_k-1
       
       
       λ_k := 0 
       ξ̃_k^* := ρ̃_k-1
       
       _k := (1-λ_k)_k-1 + λ_k_k-1
        (_k,ξ̃_k^*)


Some remarks on this procedure are as follows.  The variable
γ̂_k in eq:hatgamma1
and eq:hatgamma2 stands for γ_k-γ_k-1, where 
γ_k is defined as in eq:gammakdef.
The line-search implicit in eq:linesearch is carried
out with a univariate Newton method.  Because we have not made
sufficient assumptions about f to guarantee convergence of 
Newton's method, the Newton method is safeguarded with a bisection method.
The univariate second derivative of f needed for the Newton method
can be computed using reverse-mode automatic differentiation in time
proportional to evaluate f() (refer to <cit.>).  This
univariate second derivative is also needed in eq:qstep.

The formula for β_k in eq:hzbeta is from the
CG-Descent algorithm of Hager and Zhang <cit.>.  As mentioned
earlier, the formula for α_k appearing in eq:qstep is
based on a univariate quadratic Taylor expansion of the line-search function
at _k in the direction _k.  This formula is exact if f itself
is quadratic, but in all other cases it is speculative.  However, if it
yields a poor answer, the overall algorithm is still robust because 
when the CG step gives a poor answer, the GD algorithm serves as a backup.

The main theorem about this method, which follows from the material
presented so far, is as follows.


Assuming exact line-search in eq:linesearch, 
the Hybrid NCG algorithm produces a sequence of iterates
(_k,_k,σ_k) satisfying eq:sigma-ub and
eq:sigma-decr.  Furthermore, if f() is a quadratic function,
then Hybrid NCG produces the same sequence of iterates as linear
conjugate gradient.



We now turn to three important numerical issues with this method.
The first issue to note is that function YCOMPUTE has an “if”
statement eq:if2 not present in eq:gd.rpdef–eq:sigmadef2
when the GD algorithm was described.  
In the case of the “exact” GD algorithm, the condition of the
if-statement is guaranteed to hold as established in Section <ref>.
However, because the line-search is only approximate, eq:gradorth does
not hold exactly, and therefore the condition of the eq:if2 may occasionally
fail.  In this case, we safeguard its failure by defining _k:=_k-1
and σ_k:=σ_k-1 (so that σ̃_k is updated only due
to the decrease in the objective), i.e., we keep the same containing sphere for
the optimizer as on the previous step.

The second numerical issue concerns the computation of δ_k-1 in 
eq:compdelta.  This formula is prone to roundoff error as the algorithm
converges because _k-1 and _k-1 both tend to ^* in the limit.
In our implementation,
we addressed this issue by maintaining a separate program variable storing
the vector _k-_k.  This vector is updated using a recurrent formula
that is straightforward to derive; the recurrence updates _k-_k using
vectors that also tend to .  Given an accurate representation
of _k-_k it is clear that 
eq:compdelta can be computed without significant roundoff issues.
A similar issue and similar workaround is applied to eq:newp.

The third numerical issue concerns the subtractions in
eq:hatgamma1
and eq:hatgamma2, which are also prone to roundoff error
as _k  converges.  These errors could upset
the computation of σ̃_k.
Our implementation addressed this using “computational divided
differences”; see, e.g., <cit.>.




§ COMPUTATIONAL RESULTS



We implemented four methods: Geometric Descent (GD), Accelerated Gradient (AG),
Hybrid Nonlinear Conjugate Gradient (HyNCG, described in the preceding section), 
and NCG.  NCG stands for nonlinear conjugate gradient using the Hager-Zhang
formula for β_k given by eq:hzbeta.  (The entirety of their
NCG method is called “CG-Descent”; however, we did not implement other aspects
of CG-Descent such as the line search.)
The line search used by
GD, HyNCG and NCG is based on Newton's method and
is safeguarded with a bisection.  The techniques to address numerical
issues described in the
preceding section were applied in GD, HyNCG and NCG. (The line
search and the numerical techniques are not needed for AG).

We applied these four methods to two problem classes:
approximate BPDN and hinge-loss halfspace classification.

BPDN (basis pursuit denoising) refers to the 
unconstrained convex optimization
problem:

    min‖ A-‖̱^2 +λ‖‖_1

in which λ > 0 and A∈^m× n 
has fewer rows than columns, so that the problem
is neither strongly convex nor smooth.  However, the following
approximation (called APBDN)
is both smooth and strongly convex on any bounded domain:

    min‖ A-‖̱^2 +λ∑_i=1^n√(x_i^2+δ)

where δ>0 is a fixed scalar.  It is easy to see that as
δ→ 0, the original problem is recovered.  As
δ→ 0, ℓ→ 0 and L→∞,
where ℓ,L are the moduli of strong, smooth convexity.

In our tests of ABPDN 
we took A to be a subset of √(n) rows
of the discrete-cosine transform matrix of size n× n, where
n is an even power of 2.  (This matrix and its transpose, although
dense, can be applied in O(nlog n) operations.)  The subset of
rows was selected to be those numbered by the first m=√(n) prime
integers
in order to get reproducible
pseudorandomness in the choices.  Similarly, in 
order to obtain a pseudorandom $̱, we selected∈̱^maccording to the formulab_i=sin(i^2).  The value ofλwas fixed at10^-3in all tests; the convergence
criterion was‖∇f(_k)‖≤10^-8.
 Finally,
we variedδ=10^-2,10^-3,10^-4and we tried bothn=65536andn=262144.

The second test-case is the hinge-loss (HL) function for half-space
identification taken from <cit.>, which is as follows:f()=H(∘̱(A))+λ‖‖^2/2, whereAis a givenm×nmatrix,$̱ is a given m-vector of ± 1,  `∘' denotes Hadamard product (i.e., the entrywise
product of two vectors), λ>0 is a regularization parameter, and
H()̌=∑_i=1^mh(v_i) where 

    h(v)={[     0.5-v,      v≤ 0,; (1-v)^2/2,   v∈[0,1],;         0,      v≥ 1. ].

Minimizing f(·) corresponds to finding a hyperplane determined
by  of the form U={∈^n:^T=0} such that rows i of A,
i=1,…,m, for which
b_i=1 lie on one side of U (i.e., A(i,:)>0)
while rows i of A for which
b_i=-1 lie on the opposite side (i.e., A(i,:)<0).  
The objective function penalizes misclassified
points as well as penalizing a large value of .

This function is smooth and strongly convex.  As λ→ 0, the
strong convexity parameter ℓ vanishes.

Unlike <cit.>, who test GD applied
to this function on data sets available
on the web, we have tested the four algorithms on synthetic data for the
purpose of better control over experimental conditions.  
In our tests, m=200000, n=447 (so that n≈√(m)), 
λ=3· 10^-1, 3· 10^-2, 3· 10^-3.
For each i, i=1,…, m, b(i)=± 1 chosen at random with
probability 0.5.  If b(i)=1, then A(i,:)=[1,…,1]/√(n)+_i^T,
where _i is a noise vector chosen as a spherical Gaussian with 
covariance matrix diag(σ^2,…,σ^2), where σ=0.4.
If b(i)=-1, then A(i,:)=-[1,…,1]/√(n)+_i^T.
For these tests, the convergence test was ‖∇ f(_k)‖≤ 10^-6.

The results of all tests are shown in Table <ref>.  
For NCG and GD, the numbers in this table are the number of inner iterations (line
search steps), which is the dominant cost in these algorithms.
In the case of HyNCG, we have reported the sum of the number of CG steps
(which do not require a line-search) plus the number of inner line-search
iterations.  For AG we have
reported the number of outer iterations.
The notation DNC indicates that the algorithm did achieve the
requisite tolerance after
10^5 outer iterations.


One sees from the table that HyNCG was superior in every test case, sometimes by
a wide margin.  An unexpected feature of the table, for which we currently do not have
an explanation, is that in the case of the HL suite of problems, the 
number of iterations was nearly invariant with respect to variation in λ,
except for AG, whose running time grows steadily with decreasing λ.

To conclude this section, we also consider two hybrid algorithm
that do not use the potential.  They are as follows:
compute a step of both GD and CG, and then select the step that 
decreases either ‖∇ f(_k)‖ (denoted HyNCG/gr) or
f(_k) (denoted HyNCG/f) by the greatest amount.  The results of
this experiment are presented in Table <ref>.



The table shows that the hybrid
based on the potential outperforms the other two methods, often
by a factor of 2 and sometimes
by a large factor.
The reason for this follows from
the discussion in Section <ref>.
Although the norm of the gradient can be used (and in fact, was used
for all tests in this section) as a termination criterion, it is not
helpful for measuring progress step by step.
The two methods HyNCG/gr and HyNCG/f must carry out two evaluations
per iteration to decide which step is preferable.
In contrast, the hybrid HYNCG based on the
potential can select the CG step without trying an alternative
provided the potential shows sufficient decrease.  



§ CONCLUSIONS

We have demonstrated that a single computable potential bounds the
convergence of three algorithms, conjugate gradient, accelerated gradient
and geometric descent.  We have also pointed out other connections between
the algorithms, namely, their relationship to an idealized algorithm and
their relationship to the Bubeck-Lee-Singh lemma.
The existence of this potential enables the formulation
of a hybrid method for convex optimization that duplicates the steps of conjugate
gradient in the case of conjugate gradient but nonetheless achieves the
optimal complexity for general smooth, strongly convex problems.  
Directions for future work include the following.


  * The hybrid algorithm requires prior knowledge of ℓ,L; it would be interesting
to develop an algorithm with the same guarantees that does not need prior knowledge
of them.  Note that linear conjugate gradient does not need any such prior knowledge
of the coefficient matrix A.

  * It would be interesting to establish a theoretical result about the improved
performance of the hybrid algorithm in the case of “nearly quadratic” functions.

  * Although accelerated gradient has been extended well beyond the realm
of unconstrained smooth, strongly convex functions, none of the other
algorithms has been.  It would be interesting to extend the conjugate
gradient ideas outside this space.  
Also interesting is the extension to constrained or composite convex
minimization.  See, for example, <cit.>.

plain

