Direct Synthesis of Iterative Algorithms With Bounds on Achievable Worst-Case Convergence Rate

TL;DR

This paper introduces a novel method for directly synthesizing iterative algorithms with guaranteed worst-case convergence bounds, surpassing previous approaches by accommodating algorithms with arbitrary finite memory.

Contribution

It proposes a robust control synthesis approach to directly design algorithms with performance guarantees, extending beyond fixed-structure algorithms.

Findings

Stronger convergence bounds than previous methods.

Applicable to algorithms with arbitrary finite memory.

Demonstrates effectiveness through theoretical analysis.

Abstract

Iterative first-order methods such as gradient descent and its variants are widely used for solving optimization and machine learning problems. There has been recent interest in analytic or numerically efficient methods for computing worst-case performance bounds for such algorithms, for example over the class of strongly convex loss functions. A popular approach is to assume the algorithm has a fixed size (fixed dimension, or memory) and that its structure is parameterized by one or two hyperparameters, for example a learning rate and a momentum parameter. Then, a Lyapunov function is sought to certify robust stability and subsequent optimization can be performed to find optimal hyperparameter tunings. In the present work, we instead fix the constraints that characterize the loss function and apply techniques from robust control synthesis to directly search over algorithms. This approach yields stronger results than those previously available, since the bounds produced hold over algorithms with an arbitrary, but finite, amount of memory rather than just holding for algorithms with a prescribed structure.

Figures

Equations

$$x_{k+1}=x_k-\alpha \nabla f(x_k)+\beta (x_k-x_{k-1}).$$

$$\begin{gathered} {\left[\begin{array}{c}\nabla f(x)-\nabla f(y)\\ x-y\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}-2mL& L+m\\ L+m& -2\end{array}\right]\left[\begin{array}{c}\nabla f(x)-\nabla f(y)\\ x-y\end{array}\right]\ge 0 \\ \text{for all }x,y\in {\mathbb{R}}^n. \end{gathered}$$

$$\Vert x_{k+1}-x^{\star }\Vert \le \sqrt{1-\frac{2\alpha mL}{L+m}}\Vert x_k-x^{\star }\Vert .$$

$$\Psi +P^{\mathsf{T}}{\Theta }^{\mathsf{T}}Q+Q^{\mathsf{T}}\Theta P\prec 0.$$

$$W_P^{\mathsf{T}}\Psi W_P\prec 0\text{and}{W}_Q^{\mathsf{T}}\Psi W_Q\prec 0.$$

$$\left[\begin{array}{cc}X& X_2\\ X_2^{\mathsf{T}}& X_3\end{array}\right]\succ 0\text{and}{\left[\begin{array}{cc}X& X_2\\ X_2^{\mathsf{T}}& X_3\end{array}\right]}^{-1}=\left[\begin{array}{cc}Y& ?\\ ?& ?\end{array}\right]$$

$$X-Y^{-1}⪰0\text{if and only if}\left[\begin{array}{cc}X& I\\ I& Y\end{array}\right]⪰0.$$

$${\xi }_{k+1}$$

$$y_k$$

$${\zeta }_{k+1}$$

$$z_k$$

$$\sum _{k=0}^T{\rho }^{-2k}{z}_k^{\mathsf{T}}M{z}_k\ge 0\text{for }T=0,1,\dots$$

$$\left[\begin{array}[]{c|c}A&B\\ \hline\cr\rule{0.0pt}{11.19443pt}C&D\end{array}\right]\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle:$}}}{\vbox{\hbox{$\textstyle:$}}}{\vbox{\hbox{$\scriptstyle:$}}}{\vbox{\hbox{$\scriptscriptstyle:$}}}{=}}\left[\begin{array}[]{cc|c}A_{G}&0&B_{G}\\ B_{\Psi 1}C_{G}&A_{\Psi}&B_{\Psi 2}+B_{\Psi 1}D_{G}\\ \hline\cr\rule{0.0pt}{11.19443pt}D_{\Psi 1}C_{G}&C_{\Psi}&D_{\Psi 2}+D_{\Psi 1}D_{G}\\ \end{array}\right].$$

$$\begin{gathered} {\left[\begin{array}{cc}A& B\end{array}\right]}^{\mathsf{T}}P\left[\begin{array}{cc}A& B\end{array}\right]-{\rho }^2{\left[\begin{array}{cc}I& 0\end{array}\right]}^{\mathsf{T}}P\left[\begin{array}{cc}I& 0\end{array}\right] \\ +{\left[\begin{array}{cc}C& D\end{array}\right]}^{\mathsf{T}}M\left[\begin{array}{cc}C& D\end{array}\right]⪯0, \end{gathered}$$

$$w_{k+1}$$

$$q_{k+1}$$

$$y_k$$

$${\left[\begin{array}{c}{y}_k\\ u_k\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\left[\begin{array}{c}{y}_k\\ u_k\end{array}\right]\ge 0.$$

$$G(z)\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle:$}}}{\vbox{\hbox{$\textstyle:$}}}{\vbox{\hbox{$\scriptstyle:$}}}{\vbox{\hbox{$\scriptscriptstyle:$}}}{=}}\left[\begin{array}[]{c|c}A_{G}&B_{G}\\ \hline\cr\rule{0.0pt}{11.19443pt}C_{G}&D_{G}\end{array}\right]=\left[\begin{array}[]{cc|c}A_{K}&B_{K}&0\\ \tfrac{L+m}{2}C_{K}&1\!+\!\tfrac{L+m}{2}D_{K}&\tfrac{L-m}{2}\\ \hline\cr C_{K}&D_{K}&0\end{array}\right]$$

$${\left[\begin{array}{c}\nabla f(x)\\ x-x^{\star }\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}-2mL& L+m\\ L+m& -2\end{array}\right]\left[\begin{array}{c}\nabla f(x)\\ x-x^{\star }\end{array}\right]\ge 0\text{for all }x\in {\mathbb{R}}^n.$$

$$\left[\begin{array}[]{c|c}A&B\\ \hline\cr\rule{0.0pt}{11.19443pt}C&D\end{array}\right]\mathrel{\mathchoice{\vbox{\hbox{$\displaystyle:$}}}{\vbox{\hbox{$\textstyle:$}}}{\vbox{\hbox{$\scriptstyle:$}}}{\vbox{\hbox{$\scriptscriptstyle:$}}}{=}}\left[\begin{array}[]{c|c}A_{G}&B_{G}\\ \hline\cr\rule{0.0pt}{11.19443pt}\left[\begin{smallmatrix}C_{G}\\ 0\end{smallmatrix}\right]&\left[\begin{smallmatrix}D_{G}\\ 1\end{smallmatrix}\right]\end{array}\right].$$

$$\begin{gathered} {\left[\begin{array}{cc}{A}_G& B_G\end{array}\right]}^{\mathsf{T}}P\left[\begin{array}{cc}{A}_G& B_G\end{array}\right]-{\rho }^2{\left[\begin{array}{cc}I& 0\end{array}\right]}^{\mathsf{T}}P\left[\begin{array}{cc}I& 0\end{array}\right] \\ +{\left[\begin{array}{cc}{C}_G& D_G\\ 0& 1\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\left[\begin{array}{cc}{C}_G& D_G\\ 0& 1\end{array}\right]⪯0. \end{gathered}$$

$${\left[\begin{array}{cc}{A}_G& B_G\\ C_G& D_G\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}P& 0\\ 0& 1\end{array}\right]\left[\begin{array}{cc}{A}_G& B_G\\ C_G& D_G\end{array}\right]-\left[\begin{array}{cc}{\rho }^{2}P& 0\\ 0& 1\end{array}\right]⪯0.$$

$$\left[\begin{array}{cccc}{\rho }^{2}P& 0& A_G^{\mathsf{T}}& C_G^{\mathsf{T}}\\ 0& 1& B_G^{\mathsf{T}}& D_G^{\mathsf{T}}\\ A_G& B_G& P^{-1}& 0\\ C_G& D_G& 0& 1\end{array}\right]⪰0,\text{and}P\succ 0.$$

$$\left[\begin{array}[]{c;{2pt/2pt}c}A_{G}&B_G\\ \hdashline[2pt/2pt]C_{G}&D_G\end{array}\right]=\left[\begin{array}[]{cc;{2pt/2pt}c}0&0&0\\ 0&1&\tfrac{L-m}{2}\\ \hdashline[2pt/2pt]0&0&0\end{array}\right]\\ +\left[\begin{array}[]{cc}I&0\\ 0&\tfrac{L+m}{2}\\ \hdashline[2pt/2pt]0&1\end{array}\right]\underbrace{\begin{bmatrix}A_{K}&B_{K}\\ C_{K}&D_{K}\end{bmatrix}}_{K}\left[\begin{array}[]{cc;{2pt/2pt}c}I&0&0\\ 0&1&0\end{array}\right]$$

$$\begin{gathered} \underset{\Psi }{\underbrace{\left[\begin{array}{cccc}{\rho }^{2}P& \left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right]& \left[\begin{array}{c}0\\ 0\end{array}\right]\\ \left[\begin{array}{cc}0& 0\end{array}\right]& 1& \left[\begin{array}{cc}0& \frac{L-m}{2}\end{array}\right]& 0\\ \left[\begin{array}{cc}0& 0\\ 0& 1\end{array}\right]& \left[\begin{array}{c}0\\ \frac{L-m}{2}\end{array}\right]& P^{-1}& \left[\begin{array}{c}0\\ 0\end{array}\right]\\ \left[\begin{array}{cc}0& 0\end{array}\right]& 0& \left[\begin{array}{cc}0& 0\end{array}\right]& 1\end{array}\right]}} \\ +\mathrm{sym}\left[\begin{array}{cc}0& 0\\ 0& 0\\ \left[\begin{array}{c}I\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ \frac{L+m}{2}\end{array}\right]\\ 0& 1\end{array}\right]K\left[\begin{array}{cccc}\left[\begin{array}{cc}I& 0\end{array}\right]& 0& 0& 0\\ \left[\begin{array}{cc}0& 1\end{array}\right]& 0& 0& 0\end{array}\right]⪰0, \end{gathered}$$

$$\left[\begin{array}{cccc}I& 0& \left[\begin{array}{cc}0& 0\end{array}\right]& 0\\ 0& 1& \left[\begin{array}{cc}0& 0\end{array}\right]& 0\\ 0& 0& \left[\begin{array}{cc}0& 1\end{array}\right]& -\frac{L+m}{2}\end{array}\right]\Psi {\left[\begin{array}{cccc}I& 0& \left[\begin{array}{cc}0& 0\end{array}\right]& 0\\ 0& 1& \left[\begin{array}{cc}0& 0\end{array}\right]& 0\\ 0& 0& \left[\begin{array}{cc}0& 1\end{array}\right]& -\frac{L+m}{2}\end{array}\right]}^{\mathsf{T}}$$

$${\left[\begin{array}{ccc}\left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\end{array}\right]\\ 1& 0& 0\\ 0& I& 0\\ 0& 0& 1\end{array}\right]}^{\mathsf{T}}\Psi \left[\begin{array}{ccc}\left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 0\end{array}\right]\\ 1& 0& 0\\ 0& I& 0\\ 0& 0& 1\end{array}\right]$$

$$\left[\begin{array}{ccc}{\rho }^{2}P& \left[\begin{array}{c}0\\ 0\end{array}\right]& \left[\begin{array}{c}0\\ 1\end{array}\right]\\ \left[\begin{array}{cc}0& 0\end{array}\right]& 1& \frac{L-m}{2}\\ \left[\begin{array}{cc}0& 1\end{array}\right]& \frac{L-m}{2}& {\left[\begin{array}{c}0\\ 1\end{array}\right]}^{\mathsf{T}}{P}^{-1}\left[\begin{array}{c}0\\ 1\end{array}\right]+(\frac{L+m}{2}{)}^2\end{array}\right]$$

$$\left[\begin{array}{ccc}1& \left[\begin{array}{cc}0& \frac{L-m}{2}\end{array}\right]& 0\\ \left[\begin{array}{c}0\\ \frac{L-m}{2}\end{array}\right]& P^{-1}& \left[\begin{array}{c}0\\ 0\end{array}\right]\\ 0& \left[\begin{array}{cc}0& 0\end{array}\right]& 1\end{array}\right]$$

$$r+(\frac{L+m}{2}{)}^2-\left({\rho }^{-2}r+(\frac{L-m}{2}{)}^2\right)$$

$$mL\ge ({\rho }^{-2}-1)r$$

$${\left[\begin{array}{c}x-y\\ \varphi (x)-\varphi (y)\end{array}\right]}^{\mathsf{T}}\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\left[\begin{array}{c}x-y\\ \varphi (x)-\varphi (y)\end{array}\right]\ge 0.$$

Full text

11footnotetext: L. Lessard is with the Department of Electrical and Computer Engineering and the Wisconsin Institute for Discovery, University of Wisconsin–Madison, USA. laurent.lessard@wisc.edu22footnotetext: P. Seiler is with the Department of Electrical Engineering and Computer Science, University of Michigan, USA. pseiler@umich.edu

1 Introduction

First-order methods are a widely used class of iterative algorithms for solving optimization problems of the form $\min_{x\in\mathbb{R}^{n}}f(x)$, where $f$ is continuously differentiable and we only have access to first-order (gradient) measurements of $f$. An example of a first-order method is the heavy ball method [14], which uses iterations of the form

[TABLE]

Here, $\alpha$ and $\beta$ are parameters to be tuned. When $\beta=0$, we recover the well-known gradient descent algorithm. Algorithms such as (1) are called accelerated because when they are appropriately tuned, they can achieve faster convergence than by using gradient descent alone.

To quantify convergence, we typically make assumptions about $f$, and then bound the worst-case convergence rate of the algorithm over all functions that satisfy the assumptions. For example, suppose $f$ is strongly convex and has Lipschitz-continuous gradients, which means there exists $0<m\leq L$ such that [11, Thm. 2.1.12]

[TABLE]

If x^{\star}\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\operatorname*{\arg\min}_{x\in\mathbb{R}^{n}}f(x), then gradient descent with $0\leq\alpha\leq\tfrac{2}{L+m}$ satisfies [11, §2.1.5]

[TABLE]

In other words, the sequence of iterates $\{x_{0},x_{1},\dots\}$ converges to the global optimizer $x^{\star}$ at a geometric rate. The fastest rate bound is $\tfrac{L-m}{L+m}$, achieved when $\alpha=\tfrac{2}{L+m}$.

In an attempt to generalize these types of results and better understand the behavior of iterative algorithms and how they should be tuned, recent works have focused on a dynamical systems interpretation of optimization algorithms. In other words, the iterate sequence $\{x_{k}\}$ corresponds to the state, and the index $k$ serves the role of time. Viewed in this light, convergence of an algorithm is equivalent to the stability of the associated equilibrium point of the dynamical system.

A recently proposed approach uses robust control to analyze worst-case convergence rates [4, 10, 8, 7]. Here, the algorithm itself plays the role of the controller, and the gradient mapping $\nabla\!f$ is the uncertain plant. This viewpoint has also been extended to analyze operator-splitting methods [12] and distributed optimization algorithms [16].

Main contributions.

In the aforementioned robust control approach, the results obtained are general in the sense that they hold across all choices of parameters as in (3), but they are fundamentally constrained by the structure of the algorithm; the particular parametric form of the iterates such as (1).

In the present work, we apply tools from robust control synthesis to derive performance bounds that are independent of parametric forms such as (1) and hold over broad classes of algorithms. These results adapt the approach in [6] for the purpose of synthesizing algorithms that optimize the convergence rate.

Specifically, we consider the class of functions with sector-bounded gradients (a weaker version of (2)). We show that, for this function class, gradient descent with stepsize $\tfrac{2}{L+m}$ achieves the fastest possible worst-case convergence rate not only among all tunings of gradient descent, but among any algorithm representable as a linear time-invariant (LTI) system with finite state. This includes algorithms of the form (1) but where $x_{k+1}$ depends linearly on $\{x_{k},x_{k-1},\dots,x_{k-\ell}\}$ for some fixed $\ell$.

Our analysis also extends to the case where $\nabla\!f$ is characterized by a given and fixed integral quadratic constraint (IQC). In this case, we prove that when using the so-called off-by-one IQC introduced in [10], which is satisfied by the class of functions (2), the optimal rate over any LTI algorithm is achieved by a particular accelerated method known as the triple momentum method [17].

The remainder of the paper is organized as follows. In Section 2, we introduce the machinery from robust control synthesis that will be used throughout the rest of the paper. In Sections 4 and 5, we specialize the analysis to the sector-bounded and IQC cases, deriving the optimal LTI algorithms in both instances. We conclude in Section 7 with some closing remarks.

2 Preliminaries

Notation.

We use $\succ$ and $\succeq$ to denote matrix inequality in the definite and semi-definite sense, respectively. $X^{\mathsf{T}}$ denotes the matrix transpose. We use the short-hand notation: \operatorname{sym}X\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}X+X^{\mathsf{T}} with the convention that products are expanded as $\operatorname{sym}XY=XY+Y^{\mathsf{T}}X^{\mathsf{T}}$.

Linear matrix inequalities.

Most of the algebraic manipulations in the sequel involve linear matrix inequalities (LMI), which have now become a widespread tool in optimization and control. We refer the reader to [2] and references therein. We will make use of three fundamental results involving LMIs. The first is a basic property of the Schur complement, which can convexify certain semidefinite programs by converting them into an LMI.

Proposition 1**.**

Let $A,B,C$ be matrices of compatible dimensions such that we can define the block matrix: X\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\left(\begin{smallmatrix}A&B\\ B^{\mathsf{T}}&C\end{smallmatrix}\right) and suppose $X=X^{\mathsf{T}}$ and $C\succ 0$. Then, $X\succeq 0$ if and only if $A-BC^{-1}B^{\mathsf{T}}\succeq 0$. The result also holds when $\succeq$ is replaced by $\succ$.

The second fundamental result is called the matrix elimination lemma, by Gahinet and Apkarian [6]. This result applies to matrix inequalities with affine dependence on a matrix of decision variables and is useful for eliminating controller parameters to obtain results independent of controller state dimension.

Lemma 2** (​​[6]).**

Given a symmetric matrix $\Psi\in\mathbb{R}^{n\times n}$ and two matrices $P$, $Q$ of column dimension $n$, consider the problem of finding some matrix $\Theta$ of compatible dimensions such that

[TABLE]

Denote by $W_{P}$, $W_{Q}$ any matrices whose columns form bases for the null spaces of $P$ and $Q$ respectively. Then there exists $\Theta$ satisfying (4) if and only if

[TABLE]

The third fundamental result we will use is called the matrix completion lemma, by Packard et al. [13, Lem. 7.5]. This result describes how a pair of matrices can be augmented to satisfy an inverse relationship.

Lemma 3** (​​[13]).**

Suppose $X,Y\in\mathbb{R}^{n\times n}$, with $X=X^{\mathsf{T}}\succ 0$ and $Y=Y^{\mathsf{T}}\succ 0$. Let $m$ be a positive integer. Then there exist matrices $X_{2}\in\mathbb{R}^{n\times m}$, $X_{3}\in\mathbb{R}^{m\times m}$ such that $X_{3}=X_{3}^{\mathsf{T}}$, and

[TABLE]

if and only if $X-Y^{-1}\succeq 0$, and $\operatorname{rank}(X-Y^{-1})\leq m$.

Remark 4**.**

Applying Proposition 1, the condition in Lemma 3 can be expressed as an LMI in $X$ and $Y$:

[TABLE]

Also note that the condition $rank(X-Y^{-1})\leq m$ is trivially satisfied for any $m\geq n$.

Exponential stability with IQCs.

Finally, we will make use of an LMI for bounding the worst-case exponential convergence rate of a linear dynamical system subject to integral quadratic constraints (IQCs). Similar LMIs have appeared in [15, 9, 1].

Lemma 5**.**

Consider the discrete-time LTI system $G$, which maps $u\mapsto y$ and obeys the dynamics:

[TABLE]

Fix $\rho\in(0,1]$, and consider the auxiliary system $\Psi_{\rho}$, which maps $(y,u)\mapsto z$ and obeys the dynamics:

[TABLE]

Suppose $(y,u)$ satisfies the $\rho$-Integral Quadratic Constraint ($\rho$-IQC) defined by $(\Psi_{\rho},M)$. That is,

[TABLE]

Define state-space matrices for the map $u\mapsto z$:

[TABLE]

If there exists some $P\succ 0$ such that

[TABLE]

then the state of $G$ converges to zero at an exponential rate of $\rho$. Specifically, $\lim_{k\to\infty}\rho^{-k}\xi_{k}=0$ for all initial states $(\xi_{0},\zeta_{0})$.

3 Algorithms as control systems

Following the line of research developed in [10, 7], we formulate the problem of iterative algorithm analysis as a control problem. Specifically, we will consider algorithms set up as in Figure 1.

The iterative algorithm must contain a pure integrator, i.e. it must take the form $K(z)\frac{1}{z-1}$ where $K(z)$ is an LTI system that represents the algorithm. Assume $K(z)$ has a state-space realization $(A_{K},B_{K},C_{K},D_{K})$. Let $w\in\mathbb{R}$ and $q\in\mathbb{R}^{n_{K}}$ be the state of integrator and $K(z)$, respectively. The order of $K(z)$, denoted $n_{K}$, is unspecified at this point, i.e. the algorithm may have a finite but arbitrary amount of memory. A realization for $K(z)\frac{1}{z-1}$ from input $e_{k}$ to output $y_{k}$ is given by:

[TABLE]

This representation is very general. For example, the heavy-ball method is represented by $K(z)=\frac{-\alpha z}{(z-\beta)}$, which corresponds to $A_{K}=C_{K}=\beta$ and $B_{K}=D_{K}=-\alpha$. Meanwhile, the gradient descent algorithm corresponds to the static gain $K(z)=-\alpha$.

Dimensionality reduction.

Throughout this paper, we assume for convenience that $K(z)$ is a single-input, single-output (SISO) system. For algorithms such as heavy-ball (1) where $\alpha$ and $\beta$ are scalars, this assumption is not restrictive. Indeed, the transfer function for the heavy-ball algorithm takes the form $I_{d}\otimes K(z)\tfrac{1}{z-1}$, where $x_{k}\in\mathbb{R}^{d}$ in (1). The IQCs we will use to describe $\nabla\!f$ and hence the LMIs (8) will also have a block-Kronecker form. These facts allow us to restrict our attention to SISO $K(z)$ without any loss of generality. For further discussion on this type of dimensionality reduction, see [10, §4.2].

Loop transformation.

The case of interest in this paper is where the function $f$ is strongly convex and has Lipschitz gradients (2). This implies, in the language of robust control, that $\nabla\!f$ is sector-bounded and slope-restricted in the interval $(m,L)$. In order to streamline the rest of the paper, we will perform a loop transformation to shift this interval to $(-1,1)$. The transformation is shown in Figure 2.

The loop-shifting transformation normalizes the nonlinearity $\phi$ to lie in the sector $(-1,1)$. As a result, $u_{k}=\phi(y_{k})$ satisfies the pointwise-in-time constraint:

[TABLE]

The input to $K(z)\frac{1}{z-1}$ is transformed by the following algebraic equation: $e_{k}=\tfrac{L-m}{2}u_{k}+\tfrac{L+m}{2}y_{k}$. A realization for the transformed system $G(z)$ is obtained by combining this relation with (9). Eliminating $e_{k}$ yields the following state-space realization with state $(q_{k},w_{k})$:

[TABLE]

4 Sector-bounded nonlinearities

In this section we treat the analysis problem for the special case where $\nabla\!f$ is sector-bounded. This assumption corresponds to relaxing (2) to hold only between an arbitrary point and the optimizer of $f$. In other words, (2) reduces to:

[TABLE]

From the IQC standpoint, this means that we assume the transformed nonlinearity $\phi$ only satisfies the pointwise-in-time constraint (10).

Our goal is to find an algorithm $(A_{K},B_{K},C_{K},D_{K})$ such that the convergence rate $\rho$ is minimized. The sector IQC (10) has no dynamics, so let $z=\left(\begin{smallmatrix}y\\ u\end{smallmatrix}\right)$. In this case, $\Psi=I_{2}$ and the state matrices for the map $u\mapsto z$ in Lemma 5 simplify to:

[TABLE]

By Lemma 5, iterates converge with rate $\rho\in(0,1]$ if there exists $P\succ 0$ such that:

[TABLE]

We can rearrange the LMI above to the form:

[TABLE]

The design corresponds to synthesizing the algorithm $(A_{K},B_{K},C_{K},D_{K})$, Lyapunov matrix $P\succ 0$, and convergence rate bound $\rho\in(0,1]$. The matrix inequality (13) is nonlinear due to products of these variables. Applying Proposition 1, we can write the equivalent condition:

[TABLE]

From (11), we can write $(A_{G},B_{G},C_{G},D_{G})$ as an affine function of the controller parameters. In the equation below, dashed lines indicate block structure (not state-space notation).

[TABLE]

Substituting (15) into (14), we obtain:

[TABLE]

where $\operatorname{sym}$ is the shorthand notation we defined in Section 2 and $K$ is defined in (15). By Lemma 2, the LMI (16) is feasible if and only if a pair of conditions hold. In this case, the conditions are:

[TABLE]

where $\Psi$ is defined in (16). Carrying out the matrix multiplications and simplifying yields:

[TABLE]

Applying Proposition 1 to the first $2\times 2$ block of (17a) we obtain the following equivalent condition:

[TABLE]

where r\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\left[\begin{smallmatrix}0\\ 1\end{smallmatrix}\right]^{\mathsf{T}}\!P^{-1}\!\left[\begin{smallmatrix}0\\ 1\end{smallmatrix}\right]>0. This simplifies to:

[TABLE]

Next, Equation (17b) is feasible if and only if the first $2\times 2$ block is feasible. Applying Proposition 1 to the upper left corner of this $2\times 2$ block yields: $r\geq(\tfrac{L-m}{2})^{2}$. Combining with (18) yields $4mL\geq(\rho^{-2}-1)(L-m)^{2}$. Solving for $\rho$ yields: $\rho\geq\frac{L-m}{L+m}$, which is precisely the optimal gradient rate bound described in Section 1. Indeed, if we also set $P^{-1}=\begin{bmatrix}\varepsilon&0\\ 0&(\frac{L-m}{2})^{2}\end{bmatrix}$ for any $\varepsilon>0$, Equation (17) is satisfied. Setting $D_{K}=-\frac{2}{m+L}$ (and $A_{K}$ empty), which is the gradient method with optimal stepsize, renders (16) feasible too.

Discussion.

It is known that for gradient descent, choosing $\alpha=\tfrac{2}{m+L}$ optimizes the rate bound (3) for strongly convex functions with Lipschitz continuous gradients. Similar results have shown that the same upper bound (and optimal stepsize) holds if we constrain $f$ to be quadratic or if we relax $f$ to merely have sector-bounded gradients [10]. It is also straightforward to show that the rate bound of $\tfrac{L-m}{L+m}$ is tight, because it is attained when using this algorithm to optimize the scalar quadratic function f(x)\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}\tfrac{m}{2}(x-x^{\star})^{2}.

This bound assumes the algorithm in question is gradient descent. In this section, we showed that for the class of functions with sector-bounded gradients, gradient descent with $\alpha=\tfrac{2}{L+m}$ is in fact optimal over all algorithms that can be expressed as LTI systems (of any fixed order) in feedback with $\nabla\!f$, as in Figure 2. Put another way, when using sector-bounded gradients, we cannot improve our upper bound by using more complex algorithms.

5 Sector-bounded and slope restricted

We now extend our approach from Section 4 to the case where $\nabla\!f$ is both sector-bounded and slope-restricted. Specifically, we assume that $f$ is strongly convex with Lipschitz-continuous gradients so that $\nabla\!f$ satisfies the constraints in (2). The transformed nonlinearity $\phi$ thus satisfies the normalized constraint for all $x,y\in\mathbb{R}^{n}$:

[TABLE]

Such functions satisfy a general class of $\rho$-IQCs defined by Zames–Falb multipliers [3, 1, 5]. For now, we will assume a general IQC of the form (6), and we further break apart $z$ into its two components, thereby resulting in joint plant-IQC dynamics:

[TABLE]

Substituting $(A_{G},B_{G},C_{G},D_{G})$ from (11), obtain:

[TABLE]

This closed loop includes the algorithm (integrator and controller) and the IQC states. The closed-loop state matrices are denoted $\left(A,B,\begin{bmatrix}C_{1}\\ C_{2}\end{bmatrix},\begin{bmatrix}D_{1}\\ D_{2}\end{bmatrix}\right)$. In this case the LMI in Lemma 5 looks different because the IQC term is quadratic in $(z^{1},z^{2})$ instead of $(y,u)$. We obtain:

[TABLE]

As before, we can use Proposition 1 to convert this to a single LMI:

[TABLE]

Equation (20) will be affine in the algorithm parameters $(A_{K},B_{K},C_{K},D_{K})$ and everything will carry through as before if $C_{2}$ and $D_{2}$ do not depend on the algorithm parameters. If we examine the formulas of (19b), this amounts to requiring for example that $D_{\Psi}^{21}=0$ and $D_{\Psi}^{22}=1$. So the question becomes: is it possible to find a factorization $\Psi$ of the IQC with this property?

The general Zames–Falb $\rho$-IQC for a nonlinearity that is sector-bounded and slope-restricted on $(-1,1)$ is:

[TABLE]

Where $\hat{h}(z)$ is a discrete-time transfer function with impulse response $h_{k}$ satisfying:

[TABLE]

If we assume for now that the Zames–Falb $\rho$-IQC is causal, then a possible factorization is:

[TABLE]

Assume $\hat{h}$ is a finite impulse response filter, i.e. $\hat{h}(z)=h_{1}z^{-1}+h_{2}z^{-2}+\cdots+h_{k}z^{-k}$. Consider the realization:

[TABLE]

We see that $D_{\Psi}^{21}=0$ and $D_{\Psi}^{22}=1$, as desired. Therefore, Equation (20) holds, but the $\hat{h}$ coefficients only appear in $C_{1}$ and $C_{2}$, and the algorithm parameters only appear in an affine fashion and in the terms $A,B,C_{1},D_{1}$. The lower-left corner of (20) can be written in terms of the algorithm parameters by substituting (19). The result is (again, dashed lines indicate block structure)

[TABLE]

Let columns of $S^{\perp}$ and $T^{\perp}$ be bases for the left and right nullspaces of $S$ and $T$, respectively. These are given by:

[TABLE]

where $\begin{bmatrix}B_{1}^{\perp}&B_{2}^{\perp}\end{bmatrix}\begin{bmatrix}B_{\Psi}^{1}\\ 1\end{bmatrix}=0$. Substituting (21) into (20) and applying Lemma 2, the LMI condition (20) is feasible if and only if:

[TABLE]

Expanding and factoring the (1,1) block to expose $P$:

[TABLE]

Apply Proposition 1 to the (1,1) block of the first LMI above to obtain an equivalent LMI in Q\mathrel{\mathchoice{\vbox{\hbox{\displaystyle:}}}{\vbox{\hbox{\textstyle:}}}{\vbox{\hbox{\scriptstyle:}}}{\vbox{\hbox{\scriptscriptstyle:}}}{=}}P^{-1}. Next apply Proposition 1 to the (2,2) block of the second LMI above to obtain an equivalent LMI in $P$. This yields the following two conditions:

[TABLE]

[TABLE]

These LMIs are block $4\times 4$ and the dimensions are $(n_{K}+1+n_{\Psi}+1)$. To proceed further, we need to break $P$ and $Q$ into their $3\times 3$ blocks. These blocks have dimensions $(n_{K}+1+n_{\Psi})$. Substitute for $(S^{\perp},T^{\perp})$ and note that $C_{2}$ simplifies to: $C_{2}=\begin{bmatrix}D_{\Psi}^{21}C_{K}&D_{\Psi}^{21}D_{K}&C_{\Psi}^{2}\end{bmatrix}=\begin{bmatrix}0&0&C_{\Psi}^{2}\end{bmatrix}$. Upon further simplification, we find the first block-rows and block-columns of $P$ and $Q$ (which had dimension $n_{K}$ cancel out. The resulting LMIs are:

[TABLE]

The rank constraint is trivially satisfied if $n_{K}\geq 2$ and hence can be neglected. In this case (22) provides a set of three LMI conditions in $P$ and $Q$. In particular, $n_{K}=2$ corresponds to a third-order algorithm (two states in $K$ plus the additional integrator).

Discussion.

For each fixed value of $\rho$, if (22) is feasible for some $P$ and $Q$, then it is possible to synthesize a controller $(A_{K},B_{K},C_{K},D_{K})$ for which the rate of convergence $\rho$ is certified over the class of sector-bounded and slope-restricted functions.

Unlike previous analyses where the algorithm is fixed and one seeks to certify a certain convergence rate by searching over IQCs of a given class [4, 10, 8, 12, 16], here we fix the particular Zames-Falb IQC and instead certify a rate of convergence over all possible algorithms. This synthesis is performed numerically in the next section for the class of off-by-one IQCs.

6 Off-by-one IQC certification

A special instance of the general class of Zames–Falb multipliers discussed in Section 5 is the weighted off-by-one $\rho$-IQC [10]. For a given $\rho\in(0,1]$ and any $h_{1}\in[0,\rho^{2}]$, the function $\phi$ satisfies the $\rho$-IQC defined by:

[TABLE]

The filter $\Psi$ is parameterized by $\rho$ and the free variable $h_{1}$. It has one timestep of memory and the corresponding $\rho$-IQC captures the slope conditions that must hold between data $(y_{k},u_{k})=(y_{k},\phi(y_{k}))$ at one timestep and the previous timestep $(y_{k-1},u_{k-1})$. If $h_{1}=0$ then the constraint reduces to the sector IQC. In other words, this form of the weighted off-by-1 subsumes the sector bound used in Section 4. The more general Zames–Falb multipliers uses a filter $\Psi$ with additional memory to incorporate the slope constraints across multiple time steps.

Substituting the IQC (23) into the semidefinite program (22) and using $h_{1}=\rho^{2}$ as in [10], we obtain a convex program in $(P,Q)$ for each fixed value of $\rho$ (and for fixed values of $L$ and $m$). Performing a bisection search on $\rho$, we can obtain the minimum $\rho$ such that the LMI is feasible. We plot this result in Figure 3, where we can see that the numerical results agree perfectly with $\rho=1-\sqrt{m/L}$, which is the fastest known convergence rate achievable for strongly convex functions with Lipschitz-continuous gradients, and it is achieved by the triple momentum method [17].

7 Concluding remarks

We demonstrated that tools from robust synthesis can be leveraged for algorithm design. Unlike prior works that fix the algorithm size (fixed memory) and seek a Lyapunov candidate that guarantees a certain convergence rate over a class of functions, we instead fix the IQC used to characterize the class of functions, and then search over all possible algorithms (with finite but arbitrary memory) for the fastest certifiable rate.

We show that if $\nabla\!f$ is sector-bounded, no algorithm of any size can outperform gradient descent. We then show that if $\nabla\!f$ satisfies the weighted off-by-one $\rho$-IQC [10], no algorithm of any (finite) memory can outperform the triple momentum method. In short, we show that eliminating the control variables from the associated LMIs can lead to novel algorithm performance bounds that are independent of algorithm dimension.