Exact Worst-case Performance of First-order Methods for Composite Convex Optimization

Introduction

Consider the composite convex minimization problem

We focus on black-box oracle-based algorithms that use first-order information to approximately solve (CM) and in particular on obtaining exact and global worst-case guarantees on their performances. That is, for a given algorithm, we simultaneously seek to obtain worst-case guarantees—for example, on objective function accuracy—and an instance of (CM) for which the algorithm behaves as such. In this work, we investigate fixed-step linear first-order methods (FSLFOM), which include among others fixed-step projected, proximal, conditional, and inexact (sub)gradient methods.

This work builds on the recent idea of performance estimation, first developed by Drori and Teboulle in and followed up on by Kim and Fessler and the authors . The approach was initially tailored for obtaining upper bounds on the worst-case behavior of fixed-step gradient methods for unconstrained minimization of a single smooth convex objective function. Motivated by subsequent results (see among others ) we extend the framework of performance estimation to the composite case involving a much broader class of algorithms and function classes (see section 1.4 for more details about previous works).

Our performance estimation framework relies on formulating the worst-case computation problem as a tractable semidefinite program (SDP), which can be tackled with standard solvers . It enjoys the following attractive features:

Any primal feasible solution to this SDP leads to a lower bound on the worst-case performance of the method under consideration, by exhibiting a particular instance of (CM).

Any dual feasible solution to this SDP corresponds to an upper bound on the worst-case performance of the method under consideration, which can be converted into an explicit proof based on a combination of valid inequalities.

Finally, we use the standard notation $e_{i}$ for the unit vector having a single $1$ as its $i$th component.

2 Performance estimation problems

In , we introduced a formal definition for the performance estimation problem in the case of a black-box first-order method for unconstrained minimization of a single convex function $F$. We now generalize the performance estimation framework for handling multiple components in the objective function.

Third, we consider a method $\mathcal{M}$ whose iterates can be computed by combining past and current oracle information about $F$. This means that after the method has performed $i-1$ steps, the next iterate $x_{i}$ should be computable as a solution to an equation of the form

Note that the only unknown in this equation is $x_{i}$, and that it thus provides an implicit definition for the next step. We will see later that this assumption on $\mathcal{M}$ includes a large number of existing methods for composite optimization.

Finally, we consider a real-valued performance criterion $\mathcal{P}$ for evaluating the efficiency of the method. In what follows sequel, we assume without loss of generality that the lower the value of $\mathcal{P}$, the better the corresponding method.

The worst-case performance of method $\mathcal{M}$ on (CM) is then the optimal value of the following optimization problem, with both functions $\left\{F^{(k)}\right\}_{k\in K}$ and iterates $\left\{x_{i}\right\}_{i\in I}$ as variables, which we call a performance estimation problem (PEP):

Note that (PEP) is inherently an infinite-dimensional optimization problem, as functions $F^{(k)}$ appear as variables. However, a crucial observation is that, due to the black-box assumption on the objective components, this problem can be cast completely equivalently in a finite-dimensional fashion. Indeed, introducing the outputs of the oracle calls as variables, namely, $O^{(k)}_{i}=\mathcal{O}_{F^{(k)}}(x_{i})$ for all iterates $i\in I$ and oracles $k\in K$, we observe that steps of method $\mathcal{M}$ can be still be computed using only information contained in variables $O^{(k)}_{i}$, so that we can reformulate (PEP) as

Note the central role played by the interpolation conditions $\mathcal{O}_{F^{(k)}}(x_{i})=O^{(k)}_{i}$ $\forall i\in I$ and $k\in K$, which enforce the existence of functions $F^{(k)}$ compatible with the output of the oracles. In the next subsection we describe situations for which this formulation is tractable.

3 First-order methods and first-order convex interpolation

In the remainder of this work, we restrict ourselves to first-order oracles and methods. We now investigate the concept of (first-order) convex interpolability, in order to make existence constraints from (PEP2) tractable—more precise requirements are detailed in section 2. From the assumptions, the existence constraint for function $F^{(k)}$

which leads us to introduce the following general definition.

We conclude that identifying explicit conditions for convex interpolability by a given class of functions will be the key to eliminate the infinite-dimensional functional variables from (PEP) and transform it into a tractable estimation problem.

First-order convex interpolation was originally developed in for classes of (possibly) $L$-smooth and (possibly) $\mu$-strongly convex functions. In section 3, we extend these results to classes of functions involving simultaneously strong convexity, smoothness, gradient boundedness, and domain boundedness (for different norms). Those extensions also allow us to consider interpolation by indicator or support functions, which may among others be used for problems involving constraints.

Also, note that the notion of first-order interpolability can be adapted for nonconvex functions as well. Replacing the concept of subdifferentiability by standard differentiability can be used to study the convergence of first-order algorithms in the cases where some components $F^{(k)}$ are not convex (see section 3.4).

4 Prior work

The performance estimation approach on the same smooth unconstrained minimization is further studied in , where convex interpolation allows the derivation of an exact convex reformulation of the problem, leading to tight worst-case estimates. The obtained semidefinite formulation also forms the basis for this work.

Another recent and closely related approach for studying performances of first-order methods consists in viewing optimization algorithms as dynamical systems and to use the related stability theory in order to numerically analyze them. This idea is proposed by Lessard, Recht, and Packard in and is attractive because it requires solving a single SDP to obtain a bound that is valid for all subsequent iterations. This technique is particularly efficient for problems involving strong convexity, for which tight linear convergence rates are often recovered. However, as they aim at finding global rates of convergence, they are naturally more conservative than the general performance estimation approach.

For more details on the general topic of convergence analysis of first-order methods, we refer to the seminal books of Nemirovsky and Yudin , Polyak , Nesterov , and the more recent book of Bertsekas . Concerning the development of accelerated methods, we specifically refer to the original work of Nesterov , and to the later extensions to minimize smoothed convex functions and composite functions .

5 Paper organization and main contributions

This work is divided into three main parts. First, section 2 is concerned with putting in place the performance estimation framework for large classes of first-order algorithms, objective functions, performance criteria, and initialization conditions. The main idea of this section is to require every element of the PEP to be linearly Gram-representable (defined in Section 2.2). This section contains multiple examples of standard settings for which the methodology applies—including those covering (sub)gradient methods (along with their projected and proximal counterparts) and conditional gradient methods (CGMs).

Section 3 focuses on providing convex interpolation conditions for different classes of convex functions commonly arising in practice. Those classes include convex functions, possibly with strong convexity, smoothness, bounded domain, and bounded (sub)gradient requirements. The subclasses of indicator and support functions are also explicitly handled. Those classes of functions can all be used directly in the performance estimation framework of section 2, since their corresponding interpolation conditions are linearly Gram-representable. This section ends with an extension of the convex interpolation results to cope with smooth nonconvex functions in a linearly Gram-representable way.

In section 4, we apply our approach to several concrete first-order algorithms. We obtain improvements on the analysis of several well-known methods, either analytically or numerically, including the proximal point algorithm and the CGM. We also use those results to provide an extension of the OGM proposed by Kim and Fessler that incorporates a projection or a proximal operator to tackle constrained and composite problems.

Performance estimation framework for first-order algorithms

We start this section by formulating (f-PEP) in terms of a Gram matrix. This leads to a tractable convex formulation for (f-PEP)—once appropriate assumptions are made on the classes of objective function components, methods, performance criteria, and initialization conditions. Those assumptions are motivated by practical applications, which we also provide in the following. The main point underlying those assumptions is to ensure that every element of the PEP can be formulated in a linear way in terms of both the entries of a Gram matrix and the function values at the iterates.

We also denote by $B^{-1}P_{N}$ the matrix

Our goal for the next subsections is to show that in a lot of situations, the performance estimation problem (f-PEP) can be expressed exactly as an SDP in the $F_{N}$ and $G_{N}$ variables:

with $S$ some index set related to the constraints, and elements $a_{i},b_{i},c,D_{i}$, and $C$ of appropriate dimensions for writing the constraints and objective function linearly in terms of the Gram matrix $G_{N}$ and of the objective function values $F_{N}$.

2 Tractable formulation of the performance estimation problem

In this section, we present our main result, stating that computing the exact worst-case performance of a method on a class of functions is tractable and can, in many cases, be formulated as (SDP-PEP). We start with the concept of Gram-representability for the different ingredients of the PEP.

A class of functions is Gram-representable (resp., linearly Gram-representable) if and only if its interpolation conditions (INT) can be formulated using a finite number of convex (resp., linear) constraints involving only the matrix $G_{N}$ and the function values $F_{N}$.

The functional classes of smooth strongly convex functions, smooth convex functions with bounded (sub)gradients, and strongly convex functions with bounded domain are linearly Gram-representable. In addition, the particular subclasses of support and indicator convex functions share this same advantageous property. The details and proofs of these results are postponed to section 3.

A performance measure is Gram-representable (resp., linearly Gram-representable) if and only if it can be expressed as a concave (resp., linear) function involving only the matrix $G_{N}$ and the function values $F_{N}$.

An initialization condition is Gram-representable (resp., linearly Gram-representable) if and only if it can be expressed using a finite number of convex (resp., linear) constraints involving only the matrix $G_{N}$ and the function values $F_{N}$.

A first-order method is Gram-representable (resp., linearly Gram-representable) if and only if the computation of its iterates, implicitly defined by an equation of type (EQi), can be expressed using a finite number of convex (resp., linear) constraints involving only the matrix $G_{N}$ and the function values $F_{N}$.

where the last condition is linear in the entries of $G_{N}$.

We can now state our main results concerning Gram-representable situations.

It directly follows from Remark 6 and from the definitions of (linear) Gram-representability for the class of functions, first-order methods, performance measures, optimality condition of a solution, and initialization conditions: any solution to the corresponding optimization problem can be transformed into a particular instance of (CM), and vice versa.

The optimal value of (PEP) increases with dimension $d$. When (PEP) with Gram-representable elements attains a finite optimal value, Proposition 2.6 implies the existence of a function with dimension at most $(n+1)(N+2)$ that achieves the worst-case value.

The assumption $d\geq(n+1)(N+2)$ is referred to as the large-scale assumption in what follows. In terms of PEP, this assumption allows us to discard the nonconvex rank constraint and lead to a tractable semidefinite programming problem, which can be solved to global optimality efficiently (see, e.g., ). Without that assumption, our PEP is a nonconvex rank-constrained SDP, equivalent to a quadratic programming problem that is NP-hard in general (e.g., it has MAX-CUT and other nonconvex quadratic programs as particular cases). Approaches to handle rank constraints exist (e.g., via augmented Lagrangian techniques , via manifold optimization , or via Newton-like methods ) but in general only guarantee convergence to stationary points. This is not useful in the case of (SDP-PEP), as this only provides lower bounds on the worst-case performance.

Under the large-scale assumption, we obtain dimension-free guarantees (i.e., valid for any dimension, and tight as soon as $d\geq(n+1)(N+2)$), as is commonly found in the literature about first-order methods. In addition, we note that the dimension bound $(n+1)(N+2)$ is in fact (very) conservative for most standard algorithms—that is, the bound in the large-scale assumption can typically be significantly reduced; see Corollary 2.13.

The worst-case results provided by the SDP from Proposition 2.6 provide a tight worst-case achievable for any operator $B$ and any dual pairing ${\left<.,.\right>}$.

3 Linearly Gram-representable first-order methods

This class of first-order methods contains as particular cases the class of FSLFOM, whose iterations are defined by a linear equation (with known constant coefficients) involving the iterates and the corresponding (sub)gradients.

Note the class of FSLFOM is exactly the class of methods whose iterations can be written in the form (using first-order optimality conditions and convexity of $F^{(k)}$):

Note that coefficients in $m_{i}$ corresponding to columns describing subsequent iterates ($Bx_{j}$ and $g^{(k)}_{j}$ $\forall j>i$ and $k\in K$) must naturally be equal to zero, as well as those of columns related to the optimal solution ($Bx_{*}$ and $g^{(k)}_{*}$ $\forall k\in K$). In addition, any FSLFOM is linearly Gram-representable using the following formulation:

which is clearly linear in terms of the Gram matrix $G_{N}$. This can also easily be extended to cope with more general classes of linearly Gram-representable first-order methods, such as the following:

where $c_{i}^{(\text{low})},b_{i}^{(\text{low})}$ and $c_{i}^{(\text{up})},b_{i}^{(\text{up})}$ are some fixed parameters. Those could, for example, be used to require a sufficient decrease condition (involving function values in $F_{N}$) or to consider methods that perform an inexact computation of the next iterate in (FSLFOM), such as in

where $\epsilon_{i}\geq 0$ is some tolerance on the accuracy of the computation of $x_{i}$ in (FSLFOM).

Before going into the details of the PEPs for our class of linear fixed-step methods over different classes of convex functions, let us give several examples of methods fitting into the model provided by (FSLFOM) and (Inexact FSLFOM).

Fixed-step subgradient and gradient algorithms. Minimizing a convex function $F$ using a fixed-step subgradient method is naturally described as $x_{i}=x_{i-1}-\alpha_{i}B^{-1}g_{i-1}$, with $\alpha_{i}$ some step size and $g_{i-1}\in\partial F(x_{i-1})$. The method clearly belongs to the class of FSLFOM, and its linear Gram matrix representation can be obtained using formulation (5).

Proximal methods and proximal gradient methods. Consider a composite objective $F^{(1)}+F^{(2)}$, where $F^{(2)}$ admits a computable proximal operator. Minimizing this objective with a fixed-step proximal gradient method is usually described as performing an explicit (sub)gradient step on $F^{(1)}$ followed by a (proximal) minimization step involving $F^{(2)}$:

Optimality conditions on this last term allow writing each iteration as

Conditional gradient methods. Consider an objective function $F^{(1)}$ to be minimized over a closed convex set $Q$, whose indicator function is $F^{(2)}$. CGMs for this problem also fit the FSLFOM model. Indeed, their iterations take the following form (given a starting point $z_{0}$):

with coefficient $\lambda_{i}\in$ chosen beforehand. Using first-order necessary and sufficient optimality conditions on the intermediate optimization problem, we obtain that $y_{i}$ can be defined by the following equation:

This algorithm can also clearly be written as an FSLFOM; one only needs to merge the two sequences of iterates $y_{i}$ and $z_{i}$ into a single sequence, defining for example the iterates using $x_{2i}=z_{i}$ and $x_{2i+1}=y_{i}$ for every $i=0,1,\ldots$

Other noise models can also easily be used in the framework, such as the one proposed by d’Aspremont . However, the inexact $(\delta,L)$-oracles developed by Devolder, Glineur and Nesterov do not seem to easily fit into the approach.This is due to the fact that no necessary and sufficient interpolation conditions for functions admitting such an inexact oracle are known—that is, standard conditions are only necessary to guarantee interpolability. Using necessary conditions that are not sufficient still allows obtaining upper bounds on the worst-case behavior, but those may not be tight.

An even broader class of methods can be obtained by combining some of the above examples and/or restricting the functions to specific classes. For example, alternate projection-type algorithms are special cases of proximal methods applied to sums of convex indicator functions and hence can be represented in the FSLFOM format.

4 Simplified performance estimation problems

Note that for standard algorithms such as the above examples of FSLFOM, the SDP resulting from Proposition 2.6 can typically be further simplified, leading to a reduction in its size.

One can remove from the original SDP formulation those $p$ unnecessary points corresponding to $p$ function values in variable $F_{N}$ and to $p$ rows/columns in the Gram matrix variable $G_{N}$. Furthermore, the $N$ equations defining the iterations allow us to further substitute $N$ variables, i.e., to remove $N$ columns from $P_{N}$ and hence $N$ rows/columns from the Gram matrix variable $G_{N}$. The dimension of $G_{N}$ can finally be decreased by one, using the fact that one of the $g_{*}^{(k)}$ may also be discarded, by substituting it using the optimality condition defining $x_{*}$.

Under the assumptions of Corollary 2.13, the large-scale assumption becomes $d\geq(n+1)(N+2)-N-p{-1}$. For example, when considering methods where only the output from a single oracle (among the $n$ possible $F^{(k)}$) is used at each iteration, we have that $p=(n-1)(N+1)$, which leads to $d\geq N+n+2$.

The original SDP from Proposition 2.6 may be challenging to solve in practice, because of its potentially large size on the one hand and because it may lack an interior on the other hand. We observe that the simplified PEP described above typically improves the situation for both issues, reducing the size of the problem and solving in a lot of cases the issue of a lack of interior points.

Convex interpolation

In this section, we study convex interpolation problems for different standard classes of convex functions. The underlying motivation is to obtain discrete characterizations of convex functions commonly arising in the context of convex optimization via first-order methods. More specifically, the classes of convex functions of interest for this section are all linearly Gram-representable (see Definition 2.2). Therefore, using those classes within the performance estimation framework will lead to tractable formulations providing tightness guarantees.

The main technical tools from this section are borrowed from convex analysis; we refer to the seminal works for details.

Consider a proper, closed and convex function $f$. The main characteristics of interest for us are the following, all commonly appearing in the context of first-order convex optimization:

Alternatively, domain and gradient boundedness can be specified in terms of diameters instead of radii.

As some characteristics are incompatible with each other (e.g., gradient boundedness is incompatible with strong convexity, domain boundedness is incompatible with smoothness), we define the following three classes of functions combining specific pairs of properties.

2 Interpolation conditions

The set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$ is $\mathcal{F}_{\mu,L}$-interpolable if and only if the following set of conditions holds for every pair of indices $i\in I$ and $j\in I$:

The set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$ is $\mathcal{S}_{D,\mu}$- ($D$-bounded, $\mu$-strongly convex) (resp., $\mathcal{S}^{\prime}_{D,\mu}$-) interpolable if and only if the following set of conditions holds for every pair of indices $i\in I$ and $j\in I$:

Observe that $f$ is $\mu$-strongly convex (convex domain, and maximum of $\mu$-strongly convex functions) and that it does interpolate the set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$. First, we have

using interpolation conditions. By noting that the maximum is bigger than taking individually the component $j$, we also have that

which allows us to conclude that $f(x_{j})=f_{j}$. To obtain that $g_{j}\in\partial f(x_{j})$, let us write

This interpolation result can be used immediately to develop interpolation conditions for the class of convex functions with bounded gradient, using the conjugate duality between smoothness and strong convexity on the one hand and gradient and domain boundedness on the other hand.

The set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$ is $\mathcal{C}_{M,L}$- ($L$-smooth with $M$-bounded subgradients) (resp., $\mathcal{C}^{\prime}_{M,L}$-) interpolable if and only if the following set of conditions holds for every pair of indices $i\in I$ and $j\in I$:

which are respectively equivalent to conditions (8) and (9).

3 Indicator and support functions

The use of projection (to deal with constraints) and regularization is so recurrent in optimization that we dedicate the next lines to interpolation procedures specifically tailored to deal with them.

In our setting, an indicator function is a closed convex function taking only values and $\infty$, for which it can be shown that the domain must be a closed convex set. As explained earlier, this class of functions is particularly interesting when considering projection operators in the context of performance estimation, as a proximal step over an indicator function is equivalent to a projection on its domain.

This corresponds to a particular case of the $\mathcal{S}_{D,\mu}$- (or $\mathcal{S}^{\prime}_{D,\mu}$-) interpolation problem with $\mu=0$. Note, however, that indicator function interpolation is not completely straightforward from $\mathcal{S}^{\prime}_{D,\mu}$-interpolation, as, for example, requiring the corresponding interpolation constraints in addition to $f_{i}=0$ would not a priori guarantee that the interpolated function from Theorem 3.5 would satisfy $f(x)=0$ on $\operatorname{dom}f$.

The set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$ is $\mathcal{I}_{D}$- (resp., $\mathcal{I}^{\prime}_{D}$-) interpolable, i.e., interpolable by a $D$-bounded indicator, if and only if the following inequalities hold for every pair of indices $i\in I$ and $j\in I$:

We start with the simpler case $D=\infty$, by considering the polyhedral set

with $a_{j}=g_{j}$ and $b_{j}={\left<g_{j},x_{j}\right>}$. The construction guarantees that $x_{i}\in Q$. Indeed, by condition (10) we have ${\left<g_{j},x_{i}\right>}\leq{\left<g_{j},x_{j}\right>},$ which is equivalent to ${\left<a_{j},x_{i}\right>}\leq b_{j}$ using the definitions of $a_{j}$ and $b_{j}$, and therefore guarantees that $x_{i}\in Q$.

Support functions are very commonly used in applications. In particular, all norms, which are used for regularization, are support functions (e.g., the $l_{1}$ norm is the support function of the unit ball for ${\left\lVert.\right\rVert}_{\infty}$).

The set $\left\{(x_{i},g_{i},f_{i})\right\}_{i\in I}$ is $\mathcal{I}^{*}_{M}$- (resp., $\mathcal{I}^{\prime*}_{M}$-) interpolable, i.e., interpolable by a support function with $M$-bounded subgradients, if and only if the following inequalities hold for every pair of indices $i\in I$ and $j\in I$:

4 Smooth nonconvex interpolation

In this short section, we derive interpolation conditions for smooth, not necessarily convex, functions. Those conditions are also linearly Gram-representable and can be used to obtain tight versions of (f-PEP) for nonconvex optimization.

where the equivalences are obtained by expressing $f$ and $\nabla f$ in terms of $h$ and $\nabla h$ (or reciprocally), which proves our statement.

From Lemma 3.13 and Theorem 3.4, it is now straightforward to establish the desired interpolation conditions.

Algorithm analysis

Consider a simple model with only one convex (possibly nonsmooth) term in the objective function,

Using an observation made in Section 2.3, we see that iterations can also be written in the form of an implicit method $x_{k+1}=x_{k}-\alpha_{k+1}B^{-1}g_{k+1}$, for some $g_{k+1}\in\partial F(x_{k+1})$, and hence belong to the class (FSLFOM).

For a recent overview and motivations concerning proximal algorithms, we refer the reader to the work of Combettes and PesquetThis work among others features a large list of known proximal operators. and to the review works of Bertsekas and Parikh and Boyd . For a historical point of view on those methods, we refer to the pioneer works of Moreau and Rockafellar and the analysis of Güler .

The standard convergence result for the proximal point algorithm is provided by Güler in [18, Theorem 2.1]:

We first prove that the bound is tight. For given $N$, $R$ and step sizes $\{\alpha_{k}\}_{1\leq k\leq N}$, we consider the $l_{1}$-shaped one-dimensional function

Indeed, note that for $x\neq 0$, we have $\nabla F(x)=\textrm{sign}(x)\frac{\sqrt{B}R}{2\sum_{k=1}^{N}\alpha_{k}}$. Hence,

The proof of the upper bound is based on considering a simplified formulation of (f-PEP) for the proximal point algorithm, computing its dual and exhibiting a feasible solution to that dual. Because it is a little longer it is relegated to Appendix A.

Let us consider another convergence measure based on the residual subgradient norm. Studying a PEP similar to the one above, we obtained strong numerical evidence for the following conjecture.

In particular, the choice $g_{N}=\frac{Bx_{N-1}-Bx_{N}}{\alpha_{N}}$ is a subgradient satisfying the inequality.

Observe that this bound cannot be improved, as it is attained on the (one-dimensional) $l_{1}$-shaped function $F(x)=\frac{\sqrt{B}R\lvert x\rvert}{\sum_{k=1}^{N}\alpha_{k}}$ started from $x_{0}=-R/\sqrt{B}$. The particular choice of subgradient suggested in the theorem corresponds to the subgradient appearing in the proximal operation when written as an implicit subgradient step.

This sort of convergence results in terms of the residual (sub)gradient norm is particularly interesting when considering dual methods. In that case, the dual residual gradient norm corresponds to the primal distance to feasibility (see, e.g., ).

2 Fast gradient methods

In this section, we consider the two-term composite objective function

In the following, we call the standard fast proximal gradient method FPGM1 (FISTA ) and introduce FPGM2, a variant with slightly better guarantees, and POGM, a novel proximal version of the optimized gradient method . FPGM2 and POGM illustrate how PEPs can be used in the development of new optimization algorithms; their study in this paper remains, however, entirely numerical.

The first variants of accelerated proximal methods we are considering use a standard proximal step after an explicit gradient step for generating the so-called primary sequence $\left\{y_{k}\right\}_{k}$.

In this algorithm, we refer to coefficients $\alpha_{k}$ as inertial parameters. We use two standard variants: $\alpha^{(a)}_{k}=\frac{k-1}{k+2}$—among others proposed in —and $\alpha^{(b)}_{k}=\frac{\theta_{k-1}-1}{\theta_{k}}$, with

and $\theta_{0}=1$— see . For both variants, the standard convergence result is (see, e.g., )

2.2 New fast proximal gradient methods (FPGM2)

Secondary sequences $\left\{x_{k}\right\}$ are usually converging slightly faster than primary sequences $\left\{y_{k}\right\}$ in the unconstrained case ($F^{(2)}=0$), as observed in . However, some issues may arise with the secondary sequences of FPGM1 when applied to constrained or proximal problems: iterates may in some cases become infeasible, or the objective may become unbounded (see Table 1). We therefore propose a new variant of a fast proximal gradient method called FPGM2, also with two different step size policies, that does not suffer from theses drawbacks. Part of the underlying motivation behind FPGM2 is also the ability to generalize it later to the optimized gradient method.

The design of FPGM2 is based on two ideas: on the one hand, it should be equivalent to the standard fast gradient method in the case of smooth unconstrained convex minimization, and on the other hand, it should not move after two consecutive iterates have reached the same optimal point for (11) (i.e., $x_{k-1}=x_{k}=x_{*}$ implies $x_{k+1}=x_{*}$).

In this algorithm, we use the coefficients $\gamma_{k}=\frac{\alpha_{k}+1}{L}$. Note that we introduced two intermediate sequences: on the one hand sequence $\left\{\gamma_{k}\right\}_{k}$, corresponding to the step sizes to be taken by the proximal steps, and on the other hand sequence $\left\{z_{k}\right\}_{k}$, which keeps track of the subgradient used in the proximal steps (note that $\frac{1}{\gamma_{k}}(z_{k}-x_{k})$ corresponds to the subgradient used in the proximal step from $z_{k}$ to $x_{k}$). Although FPGM2 may look more intricate than the classical FPGM1, it is in fact simpler, as it involves only one sequence on which both implicit (proximal) and explicit (gradient) steps are being taken. Indeed, explicit steps are taken using gradient values of $F^{(1)}$ at $x_{k}$, and subgradients used in the proximal steps are subgradients of $F^{(2)}$ also at $x_{k}$. This can also be seen by rewriting the iterations of FPGM2 using the secondary sequence $\left\{x_{k}\right\}_{k}$ only, in the following way:

Comparing the different variants of FPGM2 on Figure 1 (right plot) leads to the same conclusion as for FPGM1: inertial parameters $\alpha^{(b)}$ perform slightly better than $\alpha^{(a)}$.

All finite convergence results reported in the table actually correspond to specific worst-case functions that we could identify numerically, which means that they provide rigorous lower bounds. After solving the corresponding PEPs numerically (for $L=R=1$ and $1\leq N\leq 100$), we conjecture them to be equal to the exact worst-case guarantees.

We observe that the worst-case guarantees for FPGM2 are slightly better than for FPGM1. Guarantees for the unconstrained case are slightly better than those for the constrained and proximal cases, which are equal. Note that the secondary sequence of FPGM1 is not guaranteed to be feasible in the constrained case, and that the corresponding objective value may be unbounded in the proximal case (for any $N\geq 1$).

The worst-case functions identified numerically for the unconstrained case are Huber-shaped functions . In the constrained case, we identified one-dimensional linear optimization problems of the form $\min_{x\geq 0}{c}{x}$ as worst-cases, where $c$ is a constant defined by

where $\{h_{N,j}^{(1)}\}$ correspond to the step sizes used in FPGM according to the notation introduced in (FSLFOM), under the particular choice of $t_{N,N}=1$, and $t_{N,j}=0$ for $0\leq j\leq N-1$). Finally, for the proximal case, our worst-case has function $F^{(1)}(x)=cx$ with the same $c$ as above, and function $F^{(2)}(x)$ may be chosen equal to zero for $x\geq 0$ and to $sx$ for $x<0$, for any negative value of the slope $s<0$.

3 A proximal optimized gradient method

In this section, we consider again the nonsmooth composite convex minimization problem (11). In particular, we investigate the possibility of obtaining an optimized method for this setting (i.e., a method whose worst-case performance is the best possible).

Our proposal consists in extending the optimized gradient method (OGM) developed by Kim and Fessler in , which was originally tailored for smooth unconstrained minimization ($F^{(2)}=0$). In the unconstrained smooth minimization setting, this first-order method was recently shown in to have the best achievable worst-case guarantee for the criterion $F_{N}-F_{*}$.

The new method we propose, called POGM, has been obtained by combining ideas obtained from the original OGM and the nonstandard placement of the proximal operator used for speeding up the convergence of fast proximal gradient methods (FPGM2). It was designed using the same two principles as FPGM2 (see Remark 4.4): on the one hand, it is equivalent to OGM when applied to smooth unconstrained convex minimization problems, and on the other hand, it remains at an optimal point when it reaches one.

In this algorithm, we use the sequence $\gamma_{k}=\frac{1}{L}\frac{2\theta_{k-1}+\theta_{k}-1}{\theta_{k}}$ and the inertial coefficients proposed in :

Simply trying to generalize OGM using the standard proximal step on the primary sequence $\left\{y_{i}\right\}$ (as for FPGM1) does not lead to a converging algorithm. We obtained numerical evidence, i.e., worst-case functions showing that the worst-case bound for this candidate algorithm does not decrease after each iteration (in other words, its worst-case rate is not converging to zero). Therefore we have to introduce the same idea used in FPGM2 concerning the place of the proximal operator.

We compare POGM to FPGM with inertial coefficients $\alpha^{(b)}_{k}$ in Figure 2. We obtain worst-case performances about twice better for POGM when compared to both FPGM1 and FPGM2 between $1$ and $100$ iterations. Also, we observe that the bound for POGM (equivalent to OGM when $F^{(2)}=0$) is approximately $12\%$ worse than that for OGM in the worst-case.

Of course, POGM suffers from the drawback of requiring the knowledge of the number of iterations in advance (because the rule to compute the last coefficient $\theta_{N}$ differs from the rule to compute all the previous ones). This practical disadvantage is not easily solved: if the last $\theta_{N}$ is updated with the same rule as all the previous coefficients, performance is degraded by a nonnegligible factor, rendering it even slower than FPGM (note that this is already the case for smooth unconstrained minimization ).

4 A conditional gradient method

Consider the constrained smooth convex optimization problem

The standard global convergence guarantee for this method (see e.g., [20, Theorem 1]) is

which we compare with the exact bound provided by PEP in Figure 3(a) (see section 2.3, which shows that CGM fits into the (FSLFOM) format). The numerical guarantees we obtained by solving the PEP for up to a hundred iterations are between two and three times better than the standard guarantee.

5 Alternate projection and Dykstra methods

In this section, we numerically investigate the difference between the worst-case behaviors of the standard alternate projection method (APM) for finding a point in the intersection of two convex sets and the Dykstra method (DAPM) for finding the closest point in the intersection of two convex sets. APM is a particular instance of subgradient-type descentIt can be shown that $\frac{x-\Pi_{Q_{k}}(x)}{||x-\Pi_{Q_{k}}(x)||}$ is a subgradient of the function $f(x)$ (at $x$ such that $f(x)=||x-\Pi_{Q_{k}}(x)||$). Therefore, in the case of two sets $Q_{1},Q_{2}$, and assuming that $x$ is feasible for one of the two sets (say, $Q_{1}$), a projection onto the other one corresponds to a subgradient step on $f$ with step size $||x-\Pi_{Q_{2}}(x)||$. Hence, APM is an instance of a subgradient method for $k>1$ (when $x_{k}$ is feasible for one of the two sets). applied to the problem

whose objective function is convex and nonsmooth (with Lipschitz constant $M=1$). Therefore, its expected global convergence rate is $\mathcal{O}(\frac{1}{\sqrt{N}})$ (see [15, Theorem A.1]). We compare below the convergence of both APM and DAPM with the standard lower bound for subgradient schemes $\frac{MR}{\sqrt{N+1}}$ as a reference.

Conclusion

In this work, we presented a performance estimation approach to analyze first-order algorithms for composite optimization problems. The results of were largely extended to handle both larger classes of (composite) objective functions and larger classes of first-order algorithms (also in a more general setting for handling pairs of conjugate norms).

Our contribution was essentially threefold: first, we developed specific interpolation conditions for different classes of convex and nonconvex functions; then, we exploited those interpolation conditions to formulate the exact worst-case problem for fixed-step linear first-order methods, and finally we applied that methodology to provide tight analyses for different first-order methods. Among others, we presented a new analytical guarantee for the proximal point algorithm that is twice better than previously known and improved the standard worst-case guarantee for the conditional gradient method by more than a factor of two. On the way, we also proposed an extension of the optimized gradient method proposed by Kim and Fessler that incorporates a projection or a proximal operator.

As further research, we believe this methodology should be applied to refine analyses of methods fitting in the context of fixed-step linear first-order methods, and possibly extended to handle dynamic step size rules. To this end, a possibility is to explore convex relaxations of the resulting possibly nonconvex performance estimations problems. As an example, we believe it would be interesting to analyze algorithms involving line-search, such as backtracking or Armijo–Wolfe procedures (a first step in that direction is taken in , which study the worst-case behavior of steepest descent with exact line-search). Moreover, it seems to us that the performance estimation approach could be used to refine the analyses of randomized coordinate descent-type algorithms . Performance estimation problems also opened the door for looking toward optimized methods, as proposed by Kim and Fessler for unconstrained smooth convex minimization.

Finally, algorithmic analyses using performance estimation problems are intrinsically limited by our ability to solve semidefinite problems, both numerically (when the number of iterations is large) or analytically (to obtain results valid for any number of iterations). Therefore, any idea leading to (convex) programs that are easier to solve while maintaining reasonable guarantees would be very advantageous.

Software. An easy-to-use MATLAB implementation of the approach is available at https://github.com/AdrienTaylor/Performance-Estimation-Toolbox.

References

Appendix A Proof of upper bound in Theorem 4.1

In order to express the corresponding PEP in the simplest form, we heavily rely on some straightforward simplifications of (SDP-PEP) (see Corollary 2.13 and Remark 2.14). Let us denote by $P_{N}$ the matrix containing the information harvested after $N$ iterations: $P_{N}=[g_{1}\ g_{2}\ \ldots\ g_{N}\ Bx_{0}]$ (we use the notation $g_{i}$ for subgradients $g_{i}\in\partial F(x_{i})$), and by $G_{N}$ its corresponding Gram matrix (see section 2.1). Also, we introduce the step size vectors $m_{k}$ that express each iterate $x_{k}$ in terms of $x_{0}$ and the subgradients $\{g_{i}\}_{1\leq i\leq N}$, that is $x_{k}=P_{N}m_{k}\ (k=0,\ldots,N).$ Using the standard notation $e_{i}$ for the unit vector having a single $1$ as its $i$th component, this results in the following explicit expressions for $m_{k}$: $m_{k}=e_{N+1}-\sum_{i=1}^{k}\alpha_{i}e_{i},$ along with $m_{0}=e_{N+1}$ and $m_{*}=0$ (where we assumed without loss of generality that $x_{*}=0$).

In order to perform the worst-case analysis for PPA, we now formulate the performance estimation problem (f-PEP) as the following SDP, the simplified version of (SDP-PEP) where the $x_{k}$’s ($k=1,\ldots,N$) have been substituted using the equation defining the iterates $x_{k}=x_{k-1}-\alpha_{k}B^{-1}g_{k}$:

with matrices $2A_{ij}=e_{j}(m_{i}-m_{j})^{\top\!}+(m_{i}-m_{j})e_{j}^{\top\!}$ (where $e_{*}=0$) coming from the nonsmooth convex interpolation inequalities (see condition (7)). In order to obtain an analytical upper bound for PPA, we consider the Lagrangian dual to (PPA-PEP), which is given by the following:

(where the constraint corresponding to $f_{*}$ can be discarded since it is clear that letting $f_{*}=0$ does not change the optimal solution of (PPA-PEP)). Note that the set of equality constraints can be assimilated to a set of flow constraints on a complete directed graph. That is, considering a graph where the optimum and each iterate correspond to nodes, each nonnegative $\lambda_{ij}$ corresponds to the flow on the edge going from node $j$ to node $i$ (we choose this direction by convention). This flow constraint imposes that the outgoing flow equals the ingoing flow for every node, except at the node for final iterate $N$, where the outgoing flow should be equal to $1$, and at the optimum node, where the incoming flow should be equal to $1$. We show that the following choice is a feasible point of the dual (PPA-dPEP).

and $\lambda_{ij}=0$ otherwise. First, we clearly have $\lambda_{ij}\geq 0$ and some basic computations allow us to verify that the equality constraints from (PPA-dPEP) are satisfied:

It remains to show that the corresponding dual matrix $S$ is positive semidefinite.

In order to reduce the number of indices to be used, we will note $\lambda_{i}=\lambda_{i,i+1}$ and $\mu_{i}=\lambda_{*,i}$. Then, using the equality constraints, we arrive at the following dual matrix:

Using the values of $\mu_{i}$, $\lambda_{i}$, and $\tau$ along with elementary computations allows to verify that $\forall i\in\left\{1,\ldots,N\right\}$,