PEPit: computer-assisted worst-case analyses of first-order optimization methods in Python

Introduction

Due to their low cost per iteration, first-order optimization methods became a major tool in the modern numerical optimization toolkit. Those methods are particularly well suited when targeting only low to medium accuracy solutions, and play a central role in many fields of applications that include machine learning and signal processing. Their simplicity further allows both occasional and expert users to use them. On the contrary, when it comes to their analyses (usually based on worst-case scenarios), they are mostly reserved to expert users. The main goal of this work is to allow a simpler and reproducible access to worst-case analyses for first-order methods.

PEPit is a python package enabling computer-assisted worst-case analysis of a large family of first-order optimization methods. After being provided with a first-order method and a standard problem class, the package reformulates the problem of performing a worst-case analysis as a semidefinite program (SDP). This technique is commonly referred to as performance estimation problems (PEPs) and was introduced by (Drori and Teboulle, 2014; Drori, 2014). The package uses PEPs as formalized by (Taylor et al., 2017c, a).

In short, performing a worst-case analysis of a first-order algorithm usually relies on four main ingredients: a first-order algorithm (to be analyzed), a class of problems (containing the assumptions on the function to be minimized), a performance measure (measuring the quality of the output of the algorithm under consideration; for convenience here we assume that the algorithm aims at minimizing this performance measure and our analysis aims at finding a worst-case guarantee on it), and an initial condition (measuring the quality of the initial iterate). Performing the worst-case analysis (i.e., computing worst-case scenarios) corresponds to maximizing the performance measure of the algorithm on the class of problems, under a constraint induced by the initial condition. It turns out that such optimization problems can often be solved using SDPs in the context of first-order methods.

PEPs provide a principled approach to worst-case analyses, but usually relies on potentially tedious semidefinite programming (SDP) modelling and coding steps. The PEPit package eases the access to the methodology by automatically handling the modelling part, thereby limiting the amount of time spent on this tedious task and the risk of introducing coding mistakes in the process. In short, this work allows users to (i) write their first-order algorithms nearly as they would have implemented them, and (ii) let PEPit (a) perform the modelling and coding steps, and (b) perform the worst-case analysis numerically using tools for semidefinite programming in python (Mosek, 2010; Diamond and Boyd, 2016; O’Donoghue et al., 2016).

As a result, the package enables users to easily obtain worst-case analyses for most of the standard first-order methods, class of problems, performance measures and initial conditions. This is useful to numerically verify existing convergence guarantees, as well as to ease the development of new analyses and methods. To this end, the toolbox contains tools for analyzing classical scenarios of the first-order literature: standard problem classes (such as convex functions, smooth convex functions, Lipschitz convex functions, etc.) and algorithmic operations (such as gradient, proximal, or linear optimization oracles, etc.). Finally, the package contains more than $50$ examples and is designed in an open fashion, allowing users to easily add new ingredients (such as their own problem classes, oracles, or algorithm as an example).

This paper is organized as follows. First, Section 2 exemplifies the PEP approach on a very simple example, namely computing a worst-case contraction factor for gradient descent, and shows how to code this example in PEPit. Section 3 provides a roadmap through the package, Section 4 provides three additional numerical examples (including a composite minimization problem and a stochastic one), and some concluding remarks and perspectives are drawn in Section 5.

Related works.

The PEPit package relies on performance estimation problems as formalized in Taylor et al. (2017a). It also contains some improvements and generalization to other problem and algorithmic classes such as monotone and nonexpansive operators (Ryu et al., 2020; Lieder, 2021), stochastic methods and verification of potential (or Lyapunov/energy) functions (Hu et al., 2018; Fazlyab et al., 2018; Taylor and Bach, 2019) as inspired by the related control-theoretic IQC framework (Lessard et al., 2016). The package also contains numerous examples; e.g., recent analyses and developments from (Kim and Fessler, 2016; Van Scoy et al., 2017; Gu and Yang, 2020; Kim and Fessler, 2021; Kim, 2021; Lieder, 2021; Gannot, 2021; Taylor and Drori, 2021; Abbaszadehpeivasti et al., 2021). The package can be seen as an extended open source python version of the matlab package PESTO (Taylor et al., 2017b).

PEPit on a simple example

In this section, we illustrate the use of the package for studying the worst-case properties of a standard scenario: gradient descent for minimizing a smooth strongly convex function. The goal of this elementary example is twofold. First, we want to provide the base mathematical steps enabling the use of semidefinite programming for performing worst-case analyses, together with a corresponding PEPit code. Second, we want to highlight the main ingredients that can be generalized to other problem setups (e.g., Theorem 1 below providing “interpolation conditions” for the class of smooth strongly convex functions), allowing to analyze more algorithms under different assumptions (which are listed in Section 3).

For this example, we consider the convex optimization problem

Our goal for the rest of this section is to show how to compute the smallest possible $\tau(\mu,L,\gamma)$ (often referred to as the “contraction factor”) such that

It is relatively straightforward to establish that the smallest possible $\tau(\mu,L,\gamma)$ for which (2) is valid can be computed as the worst-case value of $\|x_{1}-y_{1}\|_{2}^{2}$ when $\|x_{0}-y_{0}\|^{2}_{2}\leqslant 1$. That is, we compute $\tau(\mu,L,\gamma)$ as the optimal value to the following optimization problem:

As written in (3), this problem involves an infinite-dimensional variable $f$. Our first step towards formulating (3) as an SDP consists in reformulating it by sampling $f$ (i.e., evaluating its function value and gradient) at the two points where its gradient is evaluated:

Using Theorem 1, we can formulate the problem of computing $\tau(\mu,L,\gamma)$ as a (nonconvex) quadratic problem:

Relying on a standard trick from semidefinite programming, one can convexify this problem using a Gram representation of the variable (this is due to maximization over $d$). That is, we formulate the problem using a positive semidefinite matrix $G\succcurlyeq 0$ defined as

Using this change of variable, we arrive to

which can be solved numerically using standard tools, see, e.g., (Diamond and Boyd, 2016; Mosek, 2010). Using numerical and/or symbolical computations, one can then easily arrive to $\tau(\mu,L,\gamma)=\max\{(1-L\gamma)^{2},(1-\mu\gamma)^{2}\}$ and hence that

To understand what PEPit can do, it is crucial to understand what elements allowed to cast the worst-case analysis as such a semidefinite program (which is what we refer to as the “modelling” of the problem). In short, the SDP reformulation of the worst-case computation problem was made possible due to $4$ main ingredients (see, e.g., (Taylor et al., 2017a, Section 2.2)):

the algorithmic steps can be expressed linearly in terms of the iterates and gradient values (i.e., step-sizes do not depend on the function at hand),

the class of functions has “interpolation condition”Interpolation conditions characterize the existence of a function (that has particular function values and gradients at given points) in the considered class by a list of constraints on those gradients, points, and function values. Such interpolation theorems (see, e.g., Theorem 1) have been obtained in the literature for various problem classes, see, e.g., Taylor et al. (2017c, a); Ryu et al. (2020). that are linear in $G$ and $F$,

the performance measure is linear (or convex piecewise linear) in $G$ and $F$,

the initial condition is linear in $G$ and $F$.

Those ingredients allow using PEPs much beyond the simple setup of this section. That is, PEPs apply for performing worst-case analyses involving a variety of first-order oracles, initial conditions, performance measures, and problem classes (see Section 3).

2 Code

where $x_{n}$ and $y_{n}$ are computed from $n$ iterations of gradient descent with step-size $\gamma$ starting from respectively $x_{0}$ and $y_{0}$. As illustrated in the previous section for the case $n=1$, computing the smallest possible such $\tau(\mu,L,\gamma,n)$ is equivalent to computing the worst-case value of $\|x_{n}-y_{n}\|_{2}^{2}$ under the constraint that $\|x_{0}-y_{0}\|_{2}^{2}\leqslant 1$ (note that we naturally have that $\tau(\mu,L,\gamma,n)\leqslant(\tau(\mu,L,\gamma,1))^{n}$). This is what we do in the following lines using PEPit.

Before going into the example, we have to include the right python imports. For this example, it is necessary to perform two imports.

Initialization of PEPit.

First, we set the stage by initializing a PEP object. This object allows manipulating the forthcoming ingredients of the PEP, such as functions and iterates.

For the sake of the example, let us pick some simple values for the problem class and algorithmic parameters, for which we perform the worst-case analysis below.

Specifying the problem class.

Second, we specify our working assumptions on the function to be optimized, and instantiate a corresponding object. Here, the minimization problem at hand was of the form (1) with a smooth strongly convex function.

Algorithm initialization.

Third, we can instantiate the starting points for the two gradient methods that we will run, and specify an initial condition on those points. To this end, two starting points $x_{0}$ and $y_{0}$ are introduced, one for each trajectory, and a bound on the initial distance between those points is specified as $\|x_{0}-y_{0}\|^{2}\leqslant 1$.

Algorithm implementation.

In this fourth step, we specify the algorithm in a natural format. In this example, we simply use the iterates (which are PEPit objects) as if we had to implement gradient descent in practice using a simple loop.

Setting up a performance measure.

It is crucial for the worst-case analysis to specify the metric for which we wish to compute a worst-case performance. In this example, we wish to compute the worst-case value of $\|x_{n}-y_{n}\|^{2}$, which we specify as follows.

Solving the PEP.

The last natural stage in the process is to solve the corresponding PEP. This is done via the following line, which will ask PEPit to perform the modelling steps and to call an appropriate SDP solver (which should be installed beforehand) to perform the worst-case analysis.

Output.

Running these pieces of code (see PEPit/examples/ for the complete example) for some specific values of the parameters $n=1$, $L=1$, $\mu=.1$ and $\gamma=1$, one obtains the following output.

Note that the size of the SDP is larger than that of Section 2 ($4\times 4$ instead of $3\times 3$ in (7)) because the modelling step is done in a slightly more generic way, which might not be exploiting all specificities of the problem at hand. For more complete examples of worst-case analyses using PEPit, see Section 4.

It is also possible to run the code for different values of the parameters, as exemplified on Figure 1. This simple example allows to observe that numerical values obtained from PEPit match the worst-case guarantee (1(a)), and to optimize the step-size numerically in Figure 1(b).

PEPit: general overview and content

In this section, describe the various choices of (i) elementary oracles used in algorithms, (ii) problem or function classes, (iii) performance measures, and (iv) initial conditions, that are naturally handled by PEPit. PEPit also allows studying methods for monotone inclusions and fixed point problems, but we do not cover them in this summary. In the optimization setting, the minimization problem under consideration has the form

The base black-box optimization oracles/operations available in PEPit are the following:

linear optimization (or conditional) steps,

PEPit also allows for their slightly more general approximate/inexact and Bregman (or mirror) versions. Those oracles might be combined with a few other operations, such as exact line-search procedures.

Problem classes.

A few base classes of functions are readily available within the package:

convex functions within different classes of assumptions possibly involving bounded gradients (Lipschitz functions), bounded domains, smoothness, and/or strong convexity. Those assumptions might be combined when compatible.

Convex indicator functions, possibly with a bounded domain.

Beyond the pure optimization setting, PEPit also allows using operators within different classes of assumptions (namely: nonexpansive, Lipschitz, cocoercive, maximally monotone, and strongly monotone operators) for studying first-order methods for monotone inclusions and variational inequalities.

Performance measures and initial conditions.

An important degree of freedom of the package is to allow a large panel of performance measures and initial conditions. Essentially, everything that can be expressed linearly (or slightly beyond) in function values and quadratically in gradient/iterates (i.e., linear in the Gram representation of Section 2) might be considered. Following the notation of Section 2.2 (and denoting by $x_{\star}$ an optimal point of $F$), typical examples of initial conditions include:

$\|x_{0}-x_{\star}\|^{2}_{2}\leqslant 1$,

$\|\nabla F(x_{0})\|^{2}_{2}\leqslant 1$ (when $F$ is differentiable, otherwise similar criterion involving some subgradient of $F$ might be used),

any linear combination of the above (see, e.g., examples in the potential functions folder of the package).

Similarly, typical examples of performance measures include:

$\|\nabla F(x_{n})\|^{2}_{2}$ (when $F$ is differentiable),

linear combinations and minimum values of the above, e.g., $\min_{0\leqslant i\leqslant n}\|\nabla F(x_{i})\|^{2}_{2}$ (when $F$ is differentiable).

Examples.

PEPit contains about $50$ examples that can readily be used, instantiating the different sets of black-box oracles, problem classes, and initial condition/performance measures. Those examples can be found in the folder PEPit/examples/.

Contributing.

PEPit is designed for allowing users to easily contribute to add features to the package. Classes of functions (or operators) as well as black-box oracles can be implemented by following the contributing guidelines from the documentation. We also welcome any new example for analyzing a method/setting that is not already present in the toolbox.

A few additional numerical examples

The following section provides a few additional numerical worst-case analyses obtained through PEPit; namely an accelerated gradient method (Nesterov, 2003), an accelerated Douglas-Rachford splitting (Patrinos et al., 2014), and point-SAGA (Defazio, 2016) (a proximal method for finite-sum minimization).

For this example, we focus again on the problem of minimizing a $L$-smooth $\mu$-strongly convex function (problem (11) with $F\in\mathcal{F}_{\mu,L}$). We consider a classical accelerated gradient method with constant momentum (Nesterov, 1983, 2003). It can be described as follows for $t\in\{0,...,n-1\}$ with $y_{0}=x_{0}$:

with $\kappa=\frac{\mu}{L}$, $\alpha=\frac{1}{L}$ and $\beta=\frac{1-\sqrt{\kappa}}{1+\sqrt{\kappa}}$. We decide to compute the smallest possible $\tau(n,L,\mu)$ such that the guarantee

A comparison between the output of PEPit and (13) is presented in Figure 2(a), where we see that (13) could be slightly improved to better match the worst-case behavior of the algorithm. The corresponding code can be found in accelerated_gradient_strongly_convex.py from the PEPit/examples/unconstrained_convex_minimization/ directory.

2 Analysis of an accelerated Douglas-Rachford splitting

In this section, we provide a simple PEPit example for studying an accelerated Douglas-Rachford splitting method. This method was introduced in Patrinos et al. (2014) where a worst-case analysis is provided for quadratic minimization. We perform a worst-case analysis numerically for a slightly more general setting:

where $f_{1}$ is closed proper and convex, and $f_{2}$ is $\mu$-strongly convex and $L$-smooth. This section focuses on the following accelerated Douglas-Rachford splitting method, described in (Patrinos et al., 2014, Section 4):

When $f_{2}$ is a $L$-smooth $\mu$-strongly convex quadratic function, the following worst-case guarantee is provided by (Patrinos et al., 2014, Theorem 5):

when $\theta=\frac{1-\alpha L}{1+\alpha L}$ and $\alpha<\frac{1}{L}$. A numerical worst-case guarantee for the case $L=1$, $\mu=0.01$, $\alpha=0.9$, $\theta=\frac{1-\alpha L}{1+\alpha L}$ is provided on Figure 2(b) for a few different values of $N$, where we use (14) as a reference for comparison. For each of those values, PEPit computed a tight (up to numerical precision) worst-case value for which we are not aware of any proved analytical worst-case guarantee beyond the quadratic setting. We see that PEPit provides an improvement over this guarantee, even when the problem under consideration is not quadratic. The corresponding PEPit code for this example can be found in the file accelerated_douglas_rachford_splitting.py from the directory PEPit/examples/composite_convex_minimization/.

3 Analysis of point-SAGA

In this section, we use PEPit for studying point-SAGA (Defazio, 2016), a stochastic algorithm for finite sum minimization:

where $\gamma=\frac{\sqrt{(n-1)^{2}+4n\frac{L}{\mu}}}{2Ln}-\frac{\left(1-\frac{1}{n}\right)}{2L}$ is the step-size. In this example, we use a Lyapunov (or potential / energy) function $V(x)=\frac{1}{L\mu}\frac{1}{n}\sum_{i\leqslant n}\|\nabla f_{i}(x)-\nabla f_{i}(x_{\star})\|^{2}_{2}+\|x-x_{\star}\|^{2}_{2},$ and compute the smallest $\tau(n,L,\mu)$ such that the guarantee

We compare (15) to PEPit’s tight (up to numerical precision) output in Figure 2(c). We see that the worst-case guarantee (15) can be slightly improved, although pretty accurate, particularly for large values of the condition number. The corresponding PEPit code of this example can be found in the file point_saga.py from the directory PEPit/examples/stochastic_convex_minimization/.

Conclusion

The PEPit package, briefly described in this paper, aims at providing a simplified access to worst-case analyses of first-order optimization methods in python. For doing that, it implements the performance estimation approach while allowing to avoid the possibly heavy semidefinite programming modeling steps. The first version of the package already contains about $50$ examples of first-order methods that can be analyzed through this framework. Those examples allow either reproducing or tightening, numerically, known worst-case analyses, or to provide new ones depending on the particular method and problem class at hand.

Overall, we believe that this package allows to quickly (in)validate proofs (step towards reproducible theory) which should help both the development and the review process in optimization. We also argue that this is a nice pedagogical tool for learning algorithms together with their worst-case properties just by playing with them. Possible extensions under consideration for future versions include an option for searching for Lyapunov (or potential/energy) functions, as well as a numerical proof assistant, and to incorporate recent extensions of PEP/IQCs to distributed and decentralized optimization (Sundararajan et al., 2020; Colla and Hendrickx, 2021).

The work of B. Goujaud and A. Dieuleveut is partially supported by ANR-19-CHIA-0002-01/chaire SCAI, and Hi!Paris. A. Taylor acknowledges support from the European Research Council (grant SEQUOIA 724063). This work was partly funded by the French government under management of Agence Nationale de la Recherche as part of the “Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

References