arXiv is now an independent nonprofit! Learn more
License: CC BY 4.0
arXiv:2508.00775v2 [eess.SY] 04 Jun 2026

Learning to Optimize with Guarantees: A Complete Characterization of Linearly Convergent Algorithms

Andrea Martin andrmar@kth.se    Ian R. Manchester ian.manchester@sydney.edu.au    Luca Furieri luca.furieri@eng.ox.ac.uk School of Electrical Engineering and Computer Science, and Digital Futures, KTH Royal Institute of Technology, Sweden Australian Centre for Robotics and School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Australia Department of Engineering Sciences, University of Oxford, United Kingdom
Abstract

The design of many classical optimization algorithms is driven by the certification of linear convergence rates over classes of optimization problems. In this paper, we consider the problem of improving the average-case performance of an algorithm over a specific distribution of problem instances. While this task can be tackled by embedding trainable components into the algorithm updates, a key challenge is to preserve worst-case guarantees across the entire problem class. For classes of composite optimization problems, we show that all linearly convergent algorithms can be parametrized in terms of a baseline linearly convergent algorithm, and a set of trainable, exponentially-decaying modifications to its update rule; crucially, this parametrization excludes all—and only—the algorithms that do not converge linearly. Our results apply to improving the average-case performance of classical algorithms such as gradient descent for nonconvex, gradient-dominated functions; Nesterov’s accelerated method for smooth, strongly convex functions; and projected gradient methods for optimization over polyhedral feasible sets. We illustrate how our characterization can be used for learning to optimize with linear convergence and feasibility guarantees. Numerical results showcase benefits over classical optimizers when solving ill-conditioned systems of linear equations and running a model predictive control scheme on a linear dynamical system.

keywords:
Learning to optimize; parametric optimization; fixed-point iterations; neural network control; numerical algorithms; application of nonlinear analysis and design.
thanks: This paper was not presented at any IFAC meeting. Corresponding author: A. Martin. This work was supported by Digital Futures and the Swiss National Science Foundation (SNSF) through the Ambizione grant PZ00P2_208951.

, ,

1 Introduction

Guarantees of fast convergence are crucial whenever optimization must be executed under tight computational budgets, as in large-scale machine learning (ML) or real-time model predictive control (MPC). Worst-case linear convergence guarantees have been developed for iterative optimization algorithms over several classes of objective functions, whose structure—e.g., strong convexity, smoothness, or gradient dominance—can be exploited by first-order schemes such as standard gradient descent and Nesterov’s accelerated method [22]. A growing body of work leverages the analogy between worst-case convergence rates and robust-control techniques such as integral quadratic constraints (IQCs), leading to characterizations of accelerated algorithms with provably optimal rates across families of convex functions [15, 27, 29].

Worst-case rate guarantees establish a baseline performance in terms of the number of iterations required to achieve a certain level of precision. However, the performance of an algorithm in a specific application does not depend solely on its worst-case convergence rate. First, there exist fundamental trade-offs between the speed of convergence and the robustness of algorithmic behavior; see, for instance, the speed/covariance trade-off for accelerated methods in strongly convex optimization analyzed in [21]. This raises the question of how to appropriately define algorithm performance. A second challenge is that scenarios encountered in applications rarely span the entirety of the space of problem instances for which the worst-case guarantee is tight, resulting in overly conservative average-case performance. This introduces another trade-off: how to tailor an algorithm performance to specific instances without compromising the original uniform guarantees. A prime example of such a situation is MPC [20], where the optimization problems to be solved online often share the same objective and system dynamics constraints, and only differ in the initial state. In such cases, a solver tailored to this sub-family of problems could converge in significantly fewer iterations than a general-purpose algorithm tuned for optimal worst-case convergence rates.

The learning to optimize (L2O) literature addresses the challenge of adopting user-defined performance metrics beyond mere convergence rates and designs algorithms that are tailored to such metrics using ML. For instance, [4] proposes an algorithm performance metric that balances convergence speed with solution precision, and accordingly designs neural network (NN) update rules. However, general-purpose NN update rules come with no guarantees. Convergence with learned updates has been addressed through possibly conservative safeguarding mechanisms [11], or by limiting the use of ML components to learning initializations [25] or to hyperparameter tuning [12] of classical algorithms. These approaches demonstrate performance exceeding that of state-of-the-art classical algorithms upon training and inherit their convergence guarantees—yet the reliance on classical update rules inherently limits the expressivity of these learning-based design methods.

Beyond the optimal tuning of classical algorithms, another line of research seeks to use ML to design entirely new convergent update rules, aiming to discover application-specific shortcuts unknown to classical update rules. This has been achieved by taking simple gradient descent as a baseline and enhancing it through learned optimal deviations from such gradient-based updates. The work [19] characterizes the class of all and only those deviation functions that ensure convergence to stationary points in nonconvex, unconstrained smooth optimization, enabling learned optimization for user-defined performance metrics and outperforming finely tuned Adam [14] in NN training. The work [5] uses deep learning to train deviations from gradient descent and saturates these updates with the norm of measured gradients, ensuring convergence for composite convex optimization. The numerical studies of [19, 5] empirically demonstrate that convergence rates superior to those of classical algorithms can be achieved through training on gradient descent perturbations.

Our main goal is to address a question that has remained open in the literature on L2O. Given any algorithm with guaranteed linear convergence over a class of optimization problems, how can we improve its average-case performance over a subclass of problems of interest without sacrificing its guarantees over the entire class? A theoretical study of these trade-offs is important towards making learned optimization a standard and reliable component of linearly convergent algorithm design.

The goal of characterizing all convergent optimization algorithms for a class of problems is analogous to the canonical control-theoretic problem of parameterizing all stabilizing controllers for a given system, initially studied for linear systems [30, 3] and more recently for nonlinear and in particular neural-network systems (see, e.g. [28, 8, 18]). Structurally, our approach is analogous in that it augments an existing convergent algorithm with a decaying perturbation, whereas the Youla parameterization augments an existing stabilizing controller with a stable perturbation. Connection can also be drawn between our work and the recently-studied input-to-state stability properties of gradient algorithms applied to optimization problems satisfying the Polyak–Łojasiewicz (PL) condition [6]. In this paper we study a wider class of optimization algorithms and provide a characterization of perturbations achieving precise bounds on convergence rates.

Contributions: Given any existing optimization algorithm that achieves linear convergence to a set of fixed points at a specified rate—henceforth the baseline algorithm—our main contributions are as follows. First, we characterize conditions on the baseline algorithm under which adding exponentially decaying perturbations preserves the same linear convergence rate, up to a higher-order polynomial term. These conditions identify fundamental trade-offs between the frequency of perturbations and their worst-case impact on the linear convergence rate. Second, we establish a completeness result for linearly convergent optimization: every update rule that converges linearly at a given rate can be written as the sum of the baseline algorithm and a suitably designed exponentially decaying augmentation of its update rule that preserves linear convergence. Finally, we instantiate our results for augmenting state-of-the-art linearly convergent algorithms, notably 1) Nesterov’s accelerated gradient (NAG) [22] method for strongly convex and smooth optimization, 2) gradient descent for classes of nonconvex Polyak–Łojasiewicz (PL) functions, and 3) projected gradient descent for convex optimization with polyhedral feasible sets. Numerical examples showcase the potential of learned linearly convergent optimization in augmenting the performance of classical algorithms under tight iteration budgets.

Notation: The set of all sequences 𝐱=(x0,x1,x2,)\mathbf{x}=(x_{0},x_{1},x_{2},\ldots) where xtnx_{t}\in\mathbb{R}^{n} for all tt\in\mathbb{N} is denoted as n\ell^{n}. For a function g:nmg:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}, we write g(𝐱)=(g(x0),g(x1),)mg\left(\mathbf{x}\right)=(g(x_{0}),g(x_{1}),\ldots)\in\ell^{m}. For mm\in\mathbb{N}, we use 𝒫m(x)\mathcal{P}_{m}(x) to denote the set of positive and monotonically non-decreasing polynomials of degree at most mm in the variable xx with p(1)1p(1)\geq 1. We define the set of fixed points of an operator π\pi, assumed non-empty, as 𝙵𝚒𝚡π\mathtt{Fix}_{\pi}. For mm\in\mathbb{N} and γ(0,1)\gamma\in(0,1), we denote by exp(m,γ)\ell_{exp}(m,\gamma) the class of signals 𝐱\mathbf{x} for which there exists a polynomial pm(t)𝒫m(t)p_{m}(t)\in\mathcal{P}_{m}(t) such that |xt|pm(t)γt|x_{t}|\leq p_{m}(t)\gamma^{t}, where |||\cdot| denotes any vector norm. We write 𝕀S\mathbb{I}_{S} for the indicator function of a set SS (0 if xSx\in S, ++\infty otherwise). The distance from xx to a set SS is defined as dist(x,S)=infyS|xy|\mathrm{dist}(x,S)=\inf_{y\in S}|x-y|. For u,vnu,v\in\mathbb{R}^{n}, uvu\leq v indicates the element‐wise inequality.

2 Problem formulation

We consider composite optimization problems of the form

minxdf(x)+g(x),\displaystyle~\min_{x\in\mathbb{R}^{d}}\quad f(x)+g(x)\,, (1)

where xdx\in\mathbb{R}^{d} is the decision variable, f:df:\mathbb{R}^{d}\to\mathbb{R} is proper and β\beta-smooth, and g:d{+}g:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\} is convex, proper, and lower semi-continuous, but potentially nonsmooth. We let F(x)=f(x)+g(x)F(x)=f(x)+g(x) for brevity, and we assume that the set of optimizers 𝒳=argminxdF(x)\mathcal{X}^{\star}=\operatorname*{arg\,min}_{x\in\mathbb{R}^{d}}~F(x) is non-empty. In particular, we note that (1) subsumes constrained optimization problems of the form

minxdf0(x)\displaystyle~\min_{x\in\mathbb{R}^{d}}\quad f_{0}(x) (2a)
subjecttofi(x)0,i[1,M],\displaystyle\operatorname*{subject~to}\quad f_{i}(x)\leq 0\,,\quad\forall i\in[1,M]\,, (2b)

where f0:df_{0}:\mathbb{R}^{d}\to\mathbb{R} is β\beta-smooth, and each function fi:df_{i}:\mathbb{R}^{d}\to\mathbb{R} with i[1,M]i\in[1,M] defines a non-empty convex feasibility set 𝒳id\mathcal{X}_{i}\subseteq\mathbb{R}^{d}. In fact, one can rewrite (2) as an instance of (1) letting f(x)=f0(x)f(x)=f_{0}(x) and g(x)=maxi[1,M]𝕀𝒳i(x)g(x)=\max_{i\in[1,M]}~\mathbb{I}_{\mathcal{X}_{i}}(x).

A standard method to solve problem (1) is to analytically construct iterations of the form:

ξt+1=π(F,ξt),xt=ϕ(F,ξt),ξ0n,t,\xi_{t+1}=\pi(F,\xi_{t})\,,~x_{t}=\phi(F,\xi_{t})\,,~\xi_{0}\in\mathbb{R}^{n}\,,~t\in\mathbb{N}\,, (3)

where ξtn\xi_{t}\in\mathbb{R}^{n} is the state variable, the decision xtdx_{t}\in\mathbb{R}^{d} is the output variable, and the operator π\pi is designed so that its set of fixed points 𝙵𝚒𝚡π\mathtt{Fix}_{\pi}—the set of points ξ\xi^{\star} such that π(F,ξ)=ξ\pi(F,\xi^{\star})=\xi^{\star}—is related to 𝒳\mathcal{X}^{\star} through ϕ\phi; that is, a point x𝒳x^{\star}\in\mathcal{X}^{\star} can be reconstructed from a point ξ𝙵𝚒𝚡π\xi^{\star}\in\mathtt{Fix}_{\pi} as per x=ϕ(F,ξ)x^{\star}=\phi(F,\xi^{\star}).

A key metric for the performance of algorithms (3) when applied to a class of problems FF\in\mathcal{F} is how fast they converge to 𝙵𝚒𝚡π\mathtt{Fix}_{\pi}. Classical optimization algorithms often come with convergence guarantees that hold for the worst-case instance of FF\in\mathcal{F}. However, optimal control methods such as MPC require efficiently finding solutions to the instances of (1) that are encountered during deployment, where the objective FF is drawn from a specific distribution 𝒟\mathcal{D}_{\mathcal{F}} over the class \mathcal{F}. Motivated as such, in this work we investigate the following question.

Given a set of problem instances F𝒟F\sim\mathcal{D}_{\mathcal{F}} and a baseline algorithm π\pi to solve (1), how can we improve its average-case performance over 𝒟\mathcal{D}_{\mathcal{F}}, while retaining worst-case convergence guarantees over the entire class \mathcal{F}?

In particular, this paper characterizes algorithms 𝝂\bm{\nu} that achieve linear convergence to 𝙵𝚒𝚡π\mathtt{Fix}_{\pi} for classes of functions \mathcal{F}. We will showcase how to leverage such characterization for learning-based algorithm design in Section 4.

Definition 1.

An algorithm ξt+1=νt(F,ξt:0)\xi_{t+1}=\nu_{t}(F,\xi_{t:0}) is said to be linearly convergent to 𝙵𝚒𝚡π\mathtt{Fix}_{\pi} for \mathcal{F} with rate γ(0,1)\gamma\in(0,1) if there exists a polynomial p(t)𝒫m(t)p(t)\in\mathcal{P}_{m}(t) such that for any FF\in\mathcal{F}

dist(ξt,𝙵𝚒𝚡π)p(t)γtdist(ξ0,𝙵𝚒𝚡π),ξ0n,\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi})\leq p(t)\gamma^{t}\operatorname{dist}(\xi_{0},\mathtt{Fix}_{\pi})\,,~\forall\xi_{0}\in\mathbb{R}^{n}\,, (4)

at all times, where dist(,)\operatorname{dist}(\cdot,\cdot) is a distance function. In this case, we write that 𝛎𝚙𝙴𝚡𝚙π(m,γ\bm{\nu}\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma). When the focus is not on the polynomial order and only on the exponential convergence rate, we write 𝛎𝚙𝙴𝚡𝚙^π(γ)\bm{\nu}\in\widehat{\mathtt{pExp}}^{\pi}_{\mathcal{F}}(\gamma). Additionally, if (4) holds with a polynomial such that p(1)=1p(1)=1, then we say that 𝛎\bm{\nu} is monotonically linearly convergent to 𝙵𝚒𝚡π\mathtt{Fix}_{\pi} since the distance to the set of fixed points shrinks monotonically with the iterations, and we write 𝛎𝙴𝚡𝚙π(γ)\bm{\nu}\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma).111Note that if (4) holds with a polynomial p(t)p(t) such that p(1)=1p(1)=1, then dist(ξt,𝙵𝚒𝚡π)γdist(ξt1,𝙵𝚒𝚡π)γtdist(ξ0,𝙵𝚒𝚡π)\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi})\leq\gamma\operatorname{dist}(\xi_{t-1},\mathtt{Fix}_{\pi})\leq\dots\leq\gamma^{t}\operatorname{dist}(\xi_{0},\mathtt{Fix}_{\pi}). Hence, (4) holds with the constant polynomial p(t)=1p(t)=1 and dist(ξt,𝙵𝚒𝚡π)\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi}) shrinks monotonically.

To enhance the average-case performance of a baseline algorithm π\pi on problem instances F𝒟F\sim\mathcal{D}_{\mathcal{F}}, we aim to design augmentations vtnv_{t}\in\mathbb{R}^{n} to its update rule that do not jeopardize its convergence guarantees. Specifically, we will consider augmented update rules defined by the iterations:

ξt+1\displaystyle\xi_{t+1} =νt(F,ξt:0)=π(F,ξt)+vt(F,ξt:0),ξ0n,\displaystyle=\nu_{t}(F,\xi_{t:0})=\pi(F,\xi_{t})+v_{t}(F,\xi_{t:0}),~\xi_{0}\in\mathbb{R}^{n}, (5)

and we will establish conditions ensuring that 𝝂\bm{\nu} is linearly convergent as per (4).

3 A complete characterization of linearly convergent optimization algorithms

In this section, we establish our main results on how introducing an augmentation term vtv_{t} affects the worst-case linear convergence guarantees of a given baseline optimization algorithm π\pi. We first abstract away from the specific form of π\pi and the class of functions \mathcal{F} it is designed to optimize; we only assume that π\pi is a linearly convergent fixed-point algorithm as per Definition 1. In Section 3.1, we present corollaries that reveal several classes of problems (1) and corresponding baseline algorithms π\pi that are compatible with our framework.

Consider a baseline algorithm π𝚙𝙴𝚡𝚙^π(γ)\pi\in\widehat{\mathtt{pExp}}^{\pi}_{\mathcal{F}}(\gamma) achieving linear convergence as per Definition 1. The property (4) implies that the signal dist(ξt,𝙵𝚒𝚡π)\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi}) decays exponentially up to a polynomial factor p(t)p(t), for any initial condition ξ0n\xi_{0}\in\mathbb{R}^{n} and any objective FF\in\mathcal{F}. We now characterize to what extent injecting exponentially decaying signals vtv_{t} in the iterates of (5) can deteriorate the convergence guarantee of π\pi.

Theorem 1.

Consider the recursion (5) and assume that π𝚙𝙴𝚡𝚙π(m,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma). Choose any NN\in\mathbb{N} such that

ρ=p(N)γN<1,\rho=p(N)\gamma^{N}<1\,,

and any auxiliary signal 𝐰exp(m,ρ)\mathbf{w}\in\ell_{exp}(m,\rho). For every tt\in\mathbb{N}, construct the augmentation signal vtv_{t} in (5) as follows:

vt={wt+1N1 if t+1modN=0,0otherwise.v_{t}=\begin{cases}w_{\frac{t+1}{N}-1}\quad\text{ if }t+1\operatorname{mod}N=0\,,\\ 0\qquad\qquad\text{otherwise}\,.\end{cases} (6)

Then, the iterates of (5) satisfy:

dist(𝝃,𝙵𝚒𝚡π)exp(m+1,p(N)Nγ).\operatorname{dist}(\bm{\xi},\mathtt{Fix}_{\pi})\in\ell_{exp}(m+1,\sqrt[N]{p(N)}\gamma)\,. (7)

We report the proof of Theorem 1 in Appendix A. Theorem 1 establishes a trade-off between how often we inject an exponentially decaying perturbation—as measured by NN\in\mathbb{N}—and the degradation of the convergence rate. In particular, if π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma), then p(1)=1p(1)=1 by definition and we observe that the asymptotic rate γ\gamma does not change even with N=1N=1, as only the order of the polynomial factor in (4) is affected. This means that one can apply perturbation 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) at every time instant as per (5) without deteriorating the asymptotic convergence rate. For the general case where π𝚙𝙴𝚡𝚙π(m,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma), instead, the convergence rate increases at most to the value p(N)Nγ(γ,1)\sqrt[N]{p(N)}\gamma\in(\gamma,1). As expected, when NN tends to infinity, we recover the original rate of the baseline algorithm because limNp(N)N=1\lim_{N\to\infty}\sqrt[N]{p(N)}=1; this corresponds to the case vt=0v_{t}=0 for all tt.

It remains crucial to understand how large is the class of linearly convergent algorithms ξt=νt(F,ξt:0)\xi_{t}=\nu_{t}(F,\xi_{t:0}) that can be achieved by perturbing a baseline algorithm π\pi with an exponentially decaying learned update vtv_{t} as per (5). Our next result establishes conditions on the baseline algorithm π\pi that guarantee that any algorithm in 𝚙𝙴𝚡𝚙π(m,γ\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma) can be represented—provided that such target algorithm satisfies the following regularity condition.

Definition 2.

Given a linearly convergent algorithm 𝛎𝚙𝙴𝚡𝚙π(m,γ\bm{\nu}\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma) with iterates 𝛏\bm{\xi}, define the sequence of updates 𝐮\mathbf{u} satisfying

ut=ξt+1ξt.u_{t}=\xi_{t+1}-\xi_{t}.

We say that 𝛎\bm{\nu} is regular if this sequence of updates vanishes with the same exponential rate, i.e. 𝐮exp(m,γ)\mathbf{u}\in\ell_{exp}(m,\gamma)

In other words, the definition above excludes pathological cases of linearly convergent algorithms that can cycle indefinitely among the points in 𝙵𝚒𝚡π\mathtt{Fix}_{\pi} even when dist(ξt,𝙵𝚒𝚡π)=0\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi})=0. As an example, consider the case where the objective FF is constant, the baseline algorithm is ξt+1=π(F,ξt)=ξt\xi_{t+1}=\pi(F,\xi_{t})=\xi_{t} with ξ00\xi_{0}\neq 0, and 𝙵𝚒𝚡π=n\mathtt{Fix}_{\pi}=\mathbb{R}^{n}. In this case, the algorithm ξt+1=νt(F,ξt:0)=ξt\xi_{t+1}=\nu_{t}(F,\xi_{t:0})=-\xi_{t} converges to 𝙵𝚒𝚡π\mathtt{Fix}_{\pi} yet is not regular according to Definition 2. In fact, by definition of 𝝂\bm{\nu}, we have ut=2ξ00u_{t}=-2\xi_{0}\neq 0 when tt is even and ut=2ξ00u_{t}=2\xi_{0}\neq 0 when tt is odd, hence the resulting 𝐮\mathbf{u} is not exponentially decaying.

We are ready to present our completeness result.

Theorem 2.

Let π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma) be such that π(F,ξ)\pi(F,\xi) is Lipschitz continuous in ξ\xi. Consider the augmented algorithm 𝛎\bm{\nu} with iterates ξt\xi_{t} defined as per (5), and any target algorithm χt+1=σt(F,χt:0)\chi_{t+1}=\sigma_{t}(F,\chi_{t:0}) such that 𝛔𝚙𝙴𝚡𝚙π(m,γ\bm{\sigma}\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma). If 𝛔\bm{\sigma} is regular, there exists a sequence 𝐯(F,𝛏)exp(m,γ)\mathbf{v}(F,\bm{\xi})\in\ell_{exp}(m,\gamma) such that the iterations of 𝛎\bm{\nu} initialized with ξ0=χ0\xi_{0}=\chi_{0} are equivalent to those of 𝛔\bm{\sigma}. Additionally, the augmented algorithm 𝛎\bm{\nu} belongs to 𝚙𝙴𝚡𝚙^π(γ)\widehat{\mathtt{pExp}}^{\pi}_{\mathcal{F}}(\gamma) for any 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) with mm\in\mathbb{N}.

A few remarks are in order. First, the completeness property of Theorem 2 is key in the context of automating the design of augmented algorithms, e.g., by learning from sampled problems, as it implies that (5) encompasses all regular linearly convergent algorithms 𝝈𝚙𝙴𝚡𝚙π(m,γ\bm{\sigma}\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma). Second, when we learn an augmentation term 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) by searching over the entire space of exponentially decaying updates, it is desirable that the baseline algorithm satisfies the stronger condition π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma). This is because, if π𝚙𝙴𝚡𝚙π(m,γ)𝙴𝚡𝚙π(γ)\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma~)\setminus\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma), Theorem 1 only guarantees linear convergence—with the worsened rate p(N)Nγ\sqrt[N]{p(N)}\,\gamma—for those 𝐯\mathbf{v} chosen exactly as in (6). In other words, without π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma), most perturbations in exp(m,γ)\ell_{exp}(m,\gamma) may not preserve linear convergence. Third, the assumption that π(F,ξ)\pi(F,\xi) is Lipschitz continuous in ξ\xi is mild; we will show in the following that this condition is satisfied for baseline algorithms widely used for convex and composite optimization.

The rest of this section is dedicated to establishing how Theorem 1 and Theorem 2 can be used to augment existing solvers for composite optimization problems in the form (1) drawn from specific classes \mathcal{F}.

3.1 Results for smooth optimization

We first consider the case (1) where g(x)=0g(x)=0 for all xdx\in\mathbb{R}^{d}, leaving us with the task of minimizing a β\beta-smooth function F(x)=f(x)F(x)=f(x). Our first result focuses on classes of possibly nonconvex functions for which standard gradient descent achieves monotonic linear convergence.

Corollary 1.

Let RSIβ,μ\mathcal{F}_{RSI}^{\beta,\mu} be the class of β\beta-smooth functions satisfying the restricted secant inequality (RSI) with constant μ>0\mu>0, that is, those for which it holds

F(x)(xx)μ2dist(x,𝒳)2,xd,\nabla F(x)^{\top}(x-x^{\star})\geq\frac{\mu}{2}\operatorname{dist}(x,\mathcal{X}^{\star})^{2}\,,~\forall x\in\mathbb{R}^{d}\,, (8)

for any xargmin𝒳dist(x,𝒳)x^{\star}\in\operatorname*{arg\,min}_{\mathcal{X}^{\star}}\operatorname{dist}(x,\mathcal{X}^{\star}). Let π\pi be the gradient descent update rule π(F,ξt)=ξtηF(ξt)\pi(F,\xi_{t})=\xi_{t}-\eta\nabla F(\xi_{t}) with ξt=xt\xi_{t}=x_{t}, η=μβ2\eta=\frac{\mu}{\beta^{2}}, and γ=1μ2β2\gamma=\sqrt{1-\frac{\mu^{2}}{\beta^{2}}}. Then, any regular algorithm 𝛔𝚙𝙴𝚡𝚙RSIβ,μπ(m,γ)\bm{\sigma}\in\mathtt{pExp}_{\mathcal{F}_{RSI}^{\beta,\mu}}^{\pi}(m,\gamma) can be written as

xt+1=νt(F,xt:0)=xtηF(xt)+vt(F,xt:0),x_{t+1}=\nu_{t}(F,x_{t:0})=x_{t}-\eta\nabla F(x_{t})+v_{t}(F,x_{t:0})\,, (9)

with 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma). Vice versa, for any 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma), the algorithm (9) is such that 𝛎𝚙𝙴𝚡𝚙^RSIβ,μπ(γ)\bm{\nu}\in\widehat{\mathtt{pExp}}_{\mathcal{F}_{RSI}^{\beta,\mu}}^{\pi}(\gamma).

Proof .

By Theorem 2.1 of [17], it holds that (4) holds for the gradient descent algorithm ξt+1=π(F,ξt)=ξtηF(ξt)\xi_{t+1}=\pi(F,\xi_{t})=\xi_{t}-\eta\nabla F(\xi_{t}) with η=μβ2\eta=\frac{\mu}{\beta^{2}} with γ=1μ2β2(0,1)\gamma=\sqrt{1-\frac{\mu^{2}}{\beta^{2}}}\in(0,1). Further, we have that π(F,ξt)\pi(F,\xi_{t}) is Lipschitz continuous in ξt\xi_{t} since |π(F,x)π(F,y)|=|xyηF(x)+ηF(y)|(1+ηβ)|xy||\pi(F,x)-\pi(F,y)|=|x-y-\eta\nabla F(x)+\eta\nabla F(y)|\leq(1+\eta\beta)|x-y|, where the last inequality follows from the β\beta-smoothness of FRSIβ,μF\in\mathcal{F}_{RSI}^{\beta,\mu}. The result then follows by applying Theorem 2.

The result of Corollary 1 enables learning over the class of all the regular linearly convergent algorithms in 𝚙𝙴𝚡𝚙RSIβ,μπ(m,γ)\mathtt{pExp}_{\mathcal{F}_{RSI}^{\beta,\mu}}^{\pi}(m,\gamma), while ensuring that the augmented algorithm (9) never leaves the class 𝚙𝙴𝚡𝚙^RSIβ,μπ(γ)\widehat{\mathtt{pExp}}_{\mathcal{F}_{RSI}^{\beta,\mu}}^{\pi}(\gamma), irrespectively of how “badly” the augmentation term 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) may be chosen.

A few comments regarding the generality of the class of functions RSIβ,μ\mathcal{F}_{RSI}^{\beta,\mu} are in order. First, RSI\mathcal{F}_{RSI} encompasses certain nonconvex functions, as highlighted in [31]. Second, it holds that SCβ,μcPLβ,μRSIβ,μ\mathcal{F}_{SC}^{\beta,\mu}\subset\mathcal{F}_{cPL}^{\beta,\mu}\subset\mathcal{F}_{RSI}^{\beta,\mu} , where SCβ,μ\mathcal{F}_{SC}^{\beta,\mu} is the set of β\beta-smooth and μ\mu-strongly convex functions complying with

F(y)F(x)+F(x)(yx)+μ2|yx|2,F(y)\geq F(x)+\nabla F(x)^{\top}(y-x)+\frac{\mu}{2}|y-x|^{2}\,, (10)

for some μ>0\mu>0, and cPLβ,μ\mathcal{F}_{cPL}^{\beta,\mu} is the set of all the β\beta-smooth and convex functions that comply with the Polyak–Łojasiewicz (PL) inequality

F(x)minxdF(x)12μ|F(x)|2,F(x)-\min_{x\in\mathbb{R}^{d}}F(x)\leq\frac{1}{2\mu}|\nabla F(x)|^{2}\,, (11)

for some μ>0\mu>0.

Remark 1.

It is well known that RSIβ,μPLβ,μ24β\mathcal{F}_{RSI}^{\beta,\mu}\subseteq\mathcal{F}_{PL}^{\beta,\frac{\mu^{2}}{4\beta}}, where PLβ,μ\mathcal{F}_{PL}^{\beta,\mu} is the set of all possibly nonconvex functions satisfying (11), see [13]. For functions in PLβ,μ\mathcal{F}_{PL}^{\beta,\mu}, the gradient descent rule π(F,x)=1βF(x)\pi(F,x)=-\frac{1}{\beta}\nabla F(x) achieves linear convergence in the function value as per

F(xt)F(1μβ)t(F(x0)F),F(x_{t})-F^{\star}\leq\left(1-\frac{\mu}{\beta}\right)^{t}(F(x_{0})-F^{\star})\,,

where F=minxdF(x)F^{\star}=\min_{x\in\mathbb{R}^{d}}~F(x). However, π\pi induces a monotonically linearly convergent sequence of iterates only if the restricted secant inequality (8) also holds, see [17].

Corollary 1 ensures a complete parametrization of regular linearly convergent algorithms with the same rate γ\gamma as gradient descent for all functions in RSIβ,μ\mathcal{F}_{RSI}^{\beta,\mu}. For the special case of strongly convex functions FSCβ,μF\in\mathcal{F}_{SC}^{\beta,\mu}, one typically wants to augment ad-hoc algorithms tailored to SCβ,μ\mathcal{F}_{SC}^{\beta,\mu} such as NAG method [22], the heavy-ball method [23, 10], or optimal-rate algorithms such as those characterized in [15, 29].

Motivated as such, we show compatibility of the proposed framework with the augmentation of accelerated algorithms for objectives FSCβ,μF\in\mathcal{F}_{SC}^{\beta,\mu}.

Corollary 2.

Consider the NAG algorithm

π(F,ξt)=[1+αα10]ξt+[η0]F([1+αα]ξt),\pi(F,\xi_{t})=\begin{bmatrix}1+\alpha&-\alpha\\ 1&0\end{bmatrix}\xi_{t}+\begin{bmatrix}-\eta\\ 0\end{bmatrix}\nabla\!F\left(\begin{bmatrix}1+\alpha&-\alpha\end{bmatrix}\xi_{t}\right)\,, (12)

where ξt=[xtxt1]\xi_{t}=\begin{bmatrix}x_{t}^{\top}&x_{t-1}^{\top}\end{bmatrix}^{\top} and α0\alpha\geq 0 is the momentum coefficient. Let η=1β\eta=\frac{1}{\beta} and α=κ1κ+1\alpha=\frac{\sqrt{\kappa}-1}{\sqrt{\kappa}+1}, where κ=βμ1\kappa=\frac{\beta}{\mu}\geq 1 is the condition number. Choose any target rate degradation factor τ(1,1γ)\tau\in(1,\frac{1}{\gamma}), where γ=11κ\gamma=\sqrt{1-\frac{1}{\sqrt{\kappa}}}. Then, for any NN\in\mathbb{N} such that p(N)<τNp(N)<\tau^{N} and 𝐯\mathbf{v} constructed as per (6) using any 𝐰exp(m,τγ)\mathbf{w}\in\ell_{exp}(m,\tau\gamma), the augmented algorithm 𝛎(F,ξt:0)\bm{\nu}(F,\xi_{t:0}) defined by ξt+1=π(F,ξt)+vt,\xi_{t+1}=\pi(F,\xi_{t})+v_{t}\,, is such that 𝛎𝚙𝙴𝚡𝚙^SCβ,μπ(τγ)\bm{\nu}\in\widehat{\mathtt{pExp}}_{\mathcal{F}_{SC}^{\beta,\mu}}^{\pi}(\tau\gamma).

Proof .

It is well known that the NAG algorithm (12) applied to the class SCβ,μ\mathcal{F}_{SC}^{\beta,\mu} with the parameters α\alpha and η\eta as above is such that π𝚙𝙴𝚡𝚙π(0,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(0,\gamma), see [22, 15]. Since p(N)N<τ<1γ\sqrt[N]{p(N)}<\tau<\frac{1}{\gamma}, we have that p(N)γN<1p(N)\gamma^{N}<1 and Theorem 1 applies.

While Corollary 2 focuses on the case where NAG is used as the baseline algorithm π\pi in (5), we remark that the results extend analogously to any baseline algorithm π𝚙𝙴𝚡𝚙SCβ,μπ(m,γ)\pi\in\mathtt{pExp}_{\mathcal{F}_{SC}^{\beta,\mu}}^{\pi}(m,\gamma) such as those with optimal convergence rates designed using IQCs as per [15, 29]. As also discussed after Theorem 1, we note that enhancing accelerated algorithms, which are not monotonic in general, involves a trade-off between keeping the worst-case degradation rate τ\tau as small as possible and the frequency at which we can apply a learned update. Last, we remark that one can always impose a target τ(1,1γ)\tau\in(1,\frac{1}{\gamma}). Indeed, a large enough NN\in\mathbb{N} such that p(N)<τNp(N)<\tau^{N} always exists since, for NN large enough, the exponential term dominates over the polynomial one.

3.2 Results for nonsmooth optimization

We now turn our attention to the case (1) where the objective F(x)=f(x)+g(x)F(x)=f(x)+g(x) is nonsmooth. Our first result focuses on the class cPL,μ\mathcal{F}_{cPL}^{\!~\infty,\mu} of potentially nonsmooth proper, lower semi-continuous, convex functions that comply with the following inequality

F(x)minxdF(x)12μdist(0,F(x))2,F(x)-\min_{x\in\mathbb{R}^{d}}F(x)\leq\frac{1}{2\mu}\operatorname{dist}(0,\partial F(x))^{2}\,, (13)

where F(x)\partial F(x) is the convex subdifferential of FF at xx, defined as

F(x)={sd:F(y)F(x)+s(yx),yd}.\partial F(x)=\left\{s\in\mathbb{R}^{d}:F(y)\geq F(x)+s^{\top}(y-x)\,,~\forall y\in\mathbb{R}^{d}\right\}\,.

In particular, note that (13) corresponds to (11) when FF is differentiable.

Corollary 3.

Consider the class of functions FcPL,μF\in\mathcal{F}_{cPL}^{\!~\infty,\mu}. Let π\pi be the proximal point algorithm performing the iterations

xt+1=proxFc(xt)=minxdF(x)+12c|xxt|2,x_{t+1}=\operatorname{prox}_{F}^{c}(x_{t})=\min_{x\in\mathbb{R}^{d}}~F(x)+\frac{1}{2c}|x-x_{t}|^{2}\,, (14)

where c>0c>0. Let γ=min{11+cμ,11+c2βμ}(0,1)\gamma=\min\left\{\frac{1}{\sqrt{1+c\mu}},\frac{1}{\sqrt{1+\frac{c^{2}}{\beta\mu}}}\right\}\in(0,1). Then, any regular algorithm 𝛔𝚙𝙴𝚡𝚙cPL,μπ(m,γ)\bm{\sigma}\in\mathtt{pExp}_{\mathcal{F}_{cPL}^{\infty,\mu}}^{\pi}(m,\gamma) can be written as

xt+1=νt(F,xt:0)=proxFc(xt)+vt(F,xt:0),x_{t+1}=\nu_{t}(F,x_{t:0})=\operatorname{prox}_{F}^{c}(x_{t})+v_{t}(F,x_{t:0})\,, (15)

with 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma). Vice versa, for any 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma), the algorithm (15) is such that 𝛎𝚙𝙴𝚡𝚙^cPL,μπ(γ)\bm{\nu}\in\widehat{\mathtt{pExp}}_{\mathcal{F}_{cPL}^{\infty,\mu}}^{\pi}(\gamma).

Proof .

We begin by showing the baseline algorithm π(F,xt)\pi(F,x_{t}) in (14) is 11-Lipschitz continuous in xtx_{t}. Let u=π(F,x)u=\pi(F,x) and v=π(F,y)v=\pi(F,y). Since uu and vv are obtained by solving the minimization (14), it holds that 0F(u)+1c(ux)0\in\partial F(u)+\frac{1}{c}(u-x) and 0F(v)+1c(vy)0\in\partial F(v)+\frac{1}{c}(v-y) by the first-order optimality conditions. In other words, there exist subgradients guF(u)g_{u}\in\partial F(u) and gvF(v)g_{v}\in\partial F(v) such that xu=cgux-u=cg_{u} and yv=cguy-v=cg_{u}. Since FF is closed, convex, and proper, F\partial F is maximally monotone [24] and (gugv)(uv)0(g_{u}-g_{v})^{\top}(u-v)\geq 0. Observe that

(gugv)(uv)\displaystyle(g_{u}-g_{v})^{\top}(u-v) =1c(xu(yv))(uv)\displaystyle=\frac{1}{c}(x-u-(y-v))^{\top}(u-v)
=1c(xy)(uv)1c|uv|20.\displaystyle=\frac{1}{c}(x-y)^{\top}(u-v)-\frac{1}{c}|u-v|^{2}\geq 0\,.

Therefore, we have (xy)(uv)|uv|2(x-y)^{\top}(u-v)\geq|u-v|^{2}. On the other hand, by the Cauchy-Schwarz inequality, |xy||uv|(xy)(uv)|x-y||u-v|\geq(x-y)^{\top}(u-v). Combining the two last inequalities, we have |uv|2|xy||uv||u-v|^{2}\leq|x-y||u-v|, from which 11-Lipschitz continuity follows. Similarly to Corollary 1, the result then follows by combining our Theorem 2 with the linear convergence result of the proximal point method (14) when applied to functions FcPL,μF\in\mathcal{F}_{cPL}^{\infty,\mu} from [17, Theorem 4.2] and the definitions of error bound and quadratic growth from [13].

The result of Corollary 3 holds for any objective FcPL,μF\in\mathcal{F}_{cPL}^{\infty,\mu}. In particular, cPL,μ\mathcal{F}_{cPL}^{\infty,\mu} encompasses the class of optimization problems (1), where fSCβ,μf\in\mathcal{F}_{SC}^{\beta,\mu} and gCg\in\mathcal{F}_{C}^{\infty}, that is, g(x)g(x) is nonsmooth and convex, for which ad-hoc algorithms have been developed to exploit the structure underlying these composite problems. Our next result focuses on the case where g(x)g(x) represents the indicator function of a set of linear constraints (2b) to address constrained optimization problems of the form (2) with f0SCβ,μf_{0}\in\mathcal{F}_{SC}^{\beta,\mu}.

Corollary 4.

Consider the constrained optimization problem (2) with f0SCβ,μf_{0}\in\mathcal{F}_{SC}^{\beta,\mu} and fi(x)=Aixbif_{i}(x)=A_{i}x-b_{i} for all i[1,M]i\in[1,M] and define the feasible set 𝒳={xd:fi(x)0,i[1,M]}\mathcal{X}=\{x\in\mathbb{R}^{d}:f_{i}(x)\leq 0,~\forall i\in[1,M]\}. Let g(x)=𝕀𝒳(x)g(x)=\mathbb{I}_{\mathcal{X}}(x) and define comp\mathcal{F}_{\text{comp}} as the set of all such functions F(x)=f0(x)+g(x)F(x)=f_{0}(x)+g(x). Let π\pi be the proximal gradient descent method performing the iterations

xt+1\displaystyle x_{t+1} =minxdg(x)+12|x(xtηf(xt))|2\displaystyle=\min_{x\in\mathbb{R}^{d}}~g(x)+\frac{1}{2}|x-(x_{t}-\eta\nabla f(x_{t}))|^{2} (16)
=proxg(xtηf(xt))=proj𝒳(xtηf(xt)),\displaystyle=\operatorname{prox}_{g}(x_{t}-\eta\nabla f(x_{t}))=\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t}))\,,

where η(0,1β]\eta\in(0,\frac{1}{\beta}]. Consider any regular algorithm χt+1=σt(F,χt:0)\chi_{t+1}=\sigma_{t}(F,\chi_{t:0}) with feasible iterates χt𝒳\chi_{t}\in\mathcal{X} and such that 𝛔𝚙𝙴𝚡𝚙compπ(m,γ)\bm{\sigma}\in\mathtt{pExp}_{\mathcal{F}_{\text{comp}}}^{\pi}(m,\gamma), with γ=1ημ\gamma=1-\eta\mu. Then, there exists 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) with

AivtbiAiproj𝒳(xtηf(xt)),i[1,M],A_{i}v_{t}\leq b_{i}-A_{i}\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t}))\,,~\forall i\in[1,M]\,, (17)

at all times such that the augmented algorithm

xt+1\displaystyle x_{t+1} =νt(F,xt:0)\displaystyle=\nu_{t}(F,x_{t:0}) (18)
=proj𝒳(xtηf(xt))+vt(F,xt:0),\displaystyle=\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t}))+v_{t}(F,x_{t:0})\,,

is equivalent to 𝛔\bm{\sigma}. Vice versa, for any 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) such that (17) holds at all times and for all i[1,M]i\in[1,M], the algorithm (18) is such that the iterates xt𝒳x_{t}\in\mathcal{X} and 𝛎𝚙𝙴𝚡𝚙^compπ(γ)\bm{\nu}\in\widehat{\mathtt{pExp}}_{\mathcal{F}_{\text{comp}}}^{\pi}(\gamma).

Proof .

The baseline algorithm π\pi defined in (16) is such that π𝙴𝚡𝚙compπ(1ημ)\pi\in\mathtt{Exp}_{\mathcal{F}_{\text{comp}}}^{\pi}(1-\eta\mu) as shown in [9, Theorem 11.5]. Now, consider any regular algorithm χt+1=σt(F,χt:0)\chi_{t+1}=\sigma_{t}(F,\chi_{t:0}) with feasible iterates χt𝒳\chi_{t}\in\mathcal{X} and such that 𝛔𝚙𝙴𝚡𝚙compπ(m,γ)\bm{\sigma}\in\mathtt{pExp}_{\mathcal{F}_{\text{comp}}}^{\pi}(m,\gamma), with γ=1ημ\gamma=1-\eta\mu. Its iterates are equivalent to those of (18) by choosing

vt=π(F,χt)+σt(F,χt:0),x0=χ0.v_{t}=-\pi(F,\chi_{t})+\sigma_{t}(F,\chi_{t:0})\,,\quad x_{0}=\chi_{0}\,. (19)

Next, we verify that AivtbiAiproj𝒳(xtηf(xt))A_{i}v_{t}\leq b_{i}-A_{i}\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t})) for all i[1,M]i\in[1,M]. We have

Aivt\displaystyle A_{i}v_{t} =Aiπ(F,χt)+Aiσt(F,χt:0)\displaystyle=-A_{i}\pi(F,\chi_{t})+A_{i}\sigma_{t}(F,\chi_{t:0})
=Aiproj𝒳(χtηf(χt))+Aiχt+1.\displaystyle=-A_{i}\operatorname{proj}_{\mathcal{X}}\left(\chi_{t}-\eta\nabla f(\chi_{t})\right)+A_{i}\chi_{t+1}\,. (20)

By definition, Aiχt+1biA_{i}\chi_{t+1}\leq b_{i} because χt+1\chi_{t+1} lies in 𝒳\mathcal{X} and the claim follows by direct substitution in (20) since the sequence χt\chi_{t} is equal to xtx_{t}. Next, we verify that π(F,x)\pi(F,x) is Lipschitz in xx. Since the projection onto an affine subspace is 11-Lipschitz, it holds that

|π(F,x)π(F,y)|\displaystyle|\pi(F,x)-\pi(F,y)| |xy+ηf(y)ηf(x)|\displaystyle\leq|x-y+\eta\nabla f(y)-\eta\nabla f(x)|
(1+ηβ)|xy|.\displaystyle\leq(1+\eta\beta)|x-y|\,.

Last, analogously to the proof of Theorem 2, it holds that 𝐯exp(m,1ημ)\mathbf{v}\in\ell_{exp}(m,1-\eta\mu) because 𝛔\bm{\sigma} is regular. Vice versa, if 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma) is such that at all times AivtbiAiproj𝒳(xtηf(xt))A_{i}v_{t}\leq b_{i}-A_{i}\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t})) for all i[1,M]i\in[1,M], the iterates of (18) are feasible because

Aixt+1=Aiproj𝒳(xtηf(xt))+Aivtbi,\displaystyle A_{i}x_{t+1}=A_{i}\operatorname{proj}_{\mathcal{X}}(x_{t}-\eta\nabla f(x_{t}))+A_{i}v_{t}\leq b_{i}\,, (21)

and 𝛎𝚙𝙴𝚡𝚙^compπ(1ημ)\bm{\nu}\in\widehat{\mathtt{pExp}}_{\mathcal{F}_{\text{comp}}}^{\pi}(1-\eta\mu) by Theorem 1.

Leveraging the composite structure of (2), Corollary 4 addresses the requirement of ensuring feasibility of all iterates of (18) in optimization problems with polytopic constraints. In fact, while Corollary 3 guarantees convergence rates of the augmented algorithm 𝝂\bm{\nu}, feasibility of iterates xtx_{t} of (15) may be lost for arbitrary choices of 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma).

4 Learning to optimize with linear convergence guarantees

In this section, we illustrate a natural application of our characterization of linearly convergent algorithms to design new update rules with improved average‑case performance relative to a baseline linearly convergent algorithm while preserving its worst‑case convergence rate. We begin by formulating a meta‑optimization problem over the augmentation sequence vt(F,ξt:0)v_{t}(F,\xi_{t:0}), where we aim at minimizing the expected cost incurred by the augmented algorithm under a given data distribution while enforcing a target worst‑case convergence rate. Beyond the special case in which the algorithmic costs are quadratic in the iterates, the augmentation sequences are linear functions, and the objective functions FF\in\mathcal{F} are quadratic, the meta‑optimization is nonconvex. This motivates parameterizing vt(F,ξt:0)v_{t}(F,\xi_{t:0}) as an exponentially decaying NN update and learning its parameters via empirical algorithmic cost minimization. We conclude by discussing the significant speed-ups achieved by learned augmentations through examples drawn from linear regression and linear quadratic MPC.222The Python source code that reproduces our experiments is publicly available at https://github.com/andrea-martin/LinearlyConvergentL2O.

4.1 Augmenting average-case performance of a baseline algorithm through NN updates

Let 𝒟\mathcal{D}_{\mathcal{F}} denote a distribution of objective functions FF\in\mathcal{F} and let π𝚙𝙴𝚡𝚙π(m,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma) be a known linearly convergent baseline algorithm for \mathcal{F}. The problem of designing an augmentation signal vtv_{t} to enhance the average-case performance of π\pi relative to 𝒟\mathcal{D}_{\mathcal{F}} is expressed as

min𝒟θ\displaystyle~\min_{\mathcal{D}_{\theta}} 𝔼θ𝒟θ𝔼F𝒟ξ0𝒟0[𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(F,𝝃)]\displaystyle\mathbb{E}_{\theta\sim\mathcal{D}_{\theta}}~\mathbb{E}_{\begin{subarray}{c}F\sim\mathcal{D}_{\mathcal{F}}\\ \xi_{0}\sim\mathcal{D}_{0}\end{subarray}}\left[\mathtt{AlgoCost}(F,\bm{\xi})\right] (22a)
subjectto\displaystyle\operatorname*{subject~to} ξt+1=ξt+π(F,ξt)+vt(F,ξt:0,θ),\displaystyle\quad\xi_{t+1}=\xi_{t}+\pi(F,\xi_{t})+v_{t}(F,\xi_{t:0},\theta),\quad (22b)
𝝃exp(m¯,γ¯),F,\displaystyle\quad\bm{\xi}\in\ell_{exp}(\overline{m},\overline{\gamma}),~\forall F\in\mathcal{F}\,, (22c)

where 𝒟θ\mathcal{D}_{\theta} and 𝒟0\mathcal{D}_{0} represent distributions over the parameters θq\theta\in\mathbb{R}^{q} and the initial conditions ξ0Ξ0n\xi_{0}\in\Xi_{0}\subseteq\mathbb{R}^{n}, respectively, m¯m\overline{m}\geq m and γ¯[γ,1)\overline{\gamma}\in[\gamma,1) in (22c) specify a target convergence rate over all functions FF\in\mathcal{F}, and 𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝()\mathtt{AlgoCost}(\cdot) measures the cost incurred by the algorithm (22b) in optimizing a function FF\in\mathcal{F} when starting from an initial condition ξ0Ξ0\xi_{0}\in\Xi_{0}. We refer to [4, 16] for commonly used algorithm performance metrics.

To cast (22) as a finite-dimensional learning problem, we perform a number of standard approximations. First, we parametrize classes of distributions 𝒟θ(ψ)\mathcal{D}_{\theta}(\psi), where ψqψ\psi\in\mathbb{R}^{q_{\psi}} is a new set of trainable parameters, that is, we let θ𝒟θ(ψ)\theta\sim\mathcal{D}_{\theta}(\psi). For instance, one can parametrize the mean and variance of Gaussian distributions. Second, we approximate the expectations in (22a) with their sample-average counterparts

1NθNFN0i=1Nθj=1NFk=1N0𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(Fj,𝝃(θi,ξ0k)),\frac{1}{N_{\theta}N_{F}N_{0}}\sum_{i=1}^{N_{\theta}}~\sum_{j=1}^{N_{F}}~\sum_{k=1}^{N_{0}}~\mathtt{AlgoCost}(F_{j},\bm{\xi}(\theta_{i},{\xi_{0}}_{k}))\,, (23)

on the sets ^={F1,,FNF}\hat{\mathcal{F}}=\{F_{1},\dots,F_{N_{F}}\}, Ξ^0={Ξ01,,Ξ0N0}\hat{\Xi}_{0}=\{\Xi_{0_{1}},\dots,\Xi_{0_{N_{0}}}\}, and Θ^={θ1,,θNθ}\hat{\Theta}=\{\theta_{1},\dots,\theta_{N_{\theta}}\} with elements independently sampled from the distributions 𝒟F\mathcal{D}_{F}, 𝒟0\mathcal{D}_{0}, and 𝒟θ\mathcal{D}_{\theta}, respectively. These two approximations allow us to evaluate the empirical algorithm performance explicitly as per (23), and backpropagate through the parameters ψ\psi.

Remark 2.

The by-design convergence guarantee (22c) ensures that the trajectories 𝛏\bm{\xi} never diverge, and therefore, commonly chosen algorithmic costs [4] admit a uniform upperbound AlgoCost()C\operatorname{AlgoCost}(\cdot)\leq C over any bounded sets of parameters (θ,F,ξ0)(\theta,F,\xi_{0}), without wrapping the values between intervals defined a-priori [26]. Therefore, the use of stochastic algorithm augmentations enables a direct application of PAC-Bayes generalization bounds [2] as follows. Given a prior distribution of augmentation parameters θ𝒟θ(ψ0)\theta\sim\mathcal{D}_{\theta}(\psi_{0}) with compact support, and a trained posterior distribution of augmentation parameters θ𝒟θ(ψ)\theta\sim\mathcal{D}_{\theta}(\psi^{\star}) with compact support, one has by [2, Theorem 2.1] that, for any λ>0\lambda>0 and δ(0,1)\delta\in(0,1),

(\displaystyle\mathbb{P}\Bigg( 𝔼θ𝒟θ(ψ)[R(θ)]𝔼θ𝒟θ(ψ)[r(θ)]+λC28(NF+N0)\displaystyle\mathbb{E}_{\theta\sim\mathcal{D}_{\theta}(\psi^{\star})}[R(\theta)]\leq\mathbb{E}_{\theta\sim\mathcal{D}_{\theta}(\psi^{\star})}[r(\theta)]+\frac{\lambda C^{2}}{8\,(N_{F}+N_{0})}
+KL(𝒟θ(ψ)𝒟θ(ψ0))+log(1δ)λ)1δ,\displaystyle~+\!\frac{\mathrm{KL}\!\left(\mathcal{D}_{\theta}(\psi^{\star})\,\middle\|\,\mathcal{D}_{\theta}(\psi_{0})\right)+\log\!\left(\tfrac{1}{\delta}\right)}{\lambda}\Bigg)\geq 1-\delta\,, (24)

where R(θ)=𝔼F𝒟ξ0𝒟0[𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(F,𝛏)],R(\theta)\;=\;\mathbb{E}_{\begin{subarray}{c}F\sim\mathcal{D}_{\mathcal{F}}\\ \xi_{0}\sim\mathcal{D}_{0}\end{subarray}}\big[\,\mathtt{AlgoCost}(F,\bm{\xi})\,\big], and r(θ)=1NFN0j=1NFk=1N0𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(Fj,𝛏(θ,ξ0k))r(\theta)\;=\;\frac{1}{N_{F}N_{0}}\sum_{j=1}^{N_{F}}\sum_{k=1}^{N_{0}}\mathtt{AlgoCost}\!\big(F_{j},\bm{\xi}(\theta,{\xi_{0}}_{k})\big) is its empiric value, and KL()\mathrm{KL}(\cdot\!\parallel\!\cdot) denotes the Kullback–Leibler divergence. The bound holds for every set of samples {Fj}j=1NF\{F_{j}\}_{j=1}^{N_{F}} and {ξ0k}k=1N0\{{\xi_{0}}_{k}\}_{k=1}^{N_{0}} obtained by sampling from 𝒟\mathcal{D}_{\mathcal{F}} and 𝒟0\mathcal{D}_{0}, respectively (i.e., for any realization of these draws). The inequality (2) enables quantifying how much the performance of the algorithm—trained over a set of problems ^\hat{\mathcal{F}} and Ξ^0\hat{\Xi}_{0}—generalizes to future instances (F,ξ0)(𝒟,𝒟0)(F,\xi_{0})\sim(\mathcal{D}_{\mathcal{F}},\mathcal{D}_{0}). We refer the interested reader to more advanced PAC-Bayes bounds [2] and their specialization to learned optimization [26] for allowing optimal tuning of the quantities in (2) as part of the learning procedure.

Last, we parametrize exponentially decaying augmentations vt(F,ξt:0,θ)v_{t}(F,\xi_{t:0},\theta) such that (22c) holds for any θq\theta\in\mathbb{R}^{q}. To achieve this, we write vtv_{t} as the product of a trainable exponentially decaying sequence of scalars and a NN operator encoding the direction of vtv_{t}. Specifically, we let θ=(α,φ)\theta=(\alpha,\varphi), where α+m¯+1\alpha\in\mathbb{R}^{\overline{m}+1}_{+} are the coefficients of a polynomial and φr\varphi\in\mathbb{R}^{r} are the weights of a NN, and let

vt=(s=0m¯αsts)γ¯ttanh(Dt(F(ξt:0),F(ξt:0),ξt:0,φ)),\displaystyle v_{t}=\left(\sum_{s=0}^{\overline{m}}\alpha_{s}t^{s}\right)\overline{\gamma}^{t}\operatorname{tanh}(D_{t}(F(\xi_{t:0}),\nabla F(\xi_{t:0}),\xi_{t:0},\varphi))\,, (25)

where DtD_{t} is chosen as the output of a long short-term memory (LSTM) network with 2 layers. For ease of implementation, in the examples below we consider deterministic algorithm updates by letting θ\theta be drawn from a Dirac distribution with trainable mean. When dealing with polyhedral constraints as per Corollary 4, we use Agmon’s iterative method [1] to enforce (17) after generating vtv_{t} through (25).

4.2 Augmenting NAG for ill-conditioned systems of linear equations

In our first example, we consider the problem of solving systems of linear equations of the form Ax=bAx=b, where A0A\succ 0 and bb are sampled from a joint distribution 𝒟A,b\mathcal{D}_{A,b}. Although each instance of this problem admits the analytical solution x=A1bx^{\star}=A^{-1}b, directly computing the matrix inverse becomes numerically unstable when the condition number κ(A)\kappa(A) of AA grows very large. To address this, we instead approximate xx^{\star} using iterative methods by solving the equivalent quadratic program:

minxd|Axb|2=minxdxAAx2bAx+bb.\min_{x\in\mathbb{R}^{d}}~|Ax-b|^{2}=\min_{x\in\mathbb{R}^{d}}~x^{\top}A^{\top}Ax-2b^{\top}Ax+b^{\top}b\,. (26)

For our experiment, we assume that A=A^+δAA=\hat{A}+\delta_{A}, where A^\hat{A} corresponds to the matrix 𝚋𝚌𝚜𝚜𝚝𝚔𝟶𝟸\mathtt{bcsstk02} in the dataset [7] and is such that λmin(A^A^)17.75\lambda_{\operatorname{min}}(\hat{A}^{\top}\hat{A})\approx 17.75, κ(A^A^)1.17×108\kappa(\hat{A}^{\top}\hat{A})\approx 1.17\times 10^{8}, and δA\delta_{A} is a matrix with entries drawn from a standard Gaussian distribution. Similarly, we assume that b=b^+δbb=\hat{b}+\delta_{b}, where b^\hat{b} is a vector with entries b^i=0.5\hat{b}_{i}=0.5 and δb\delta_{b} is a vector with entries drawn from a zero-mean Gaussian random variable with standard deviation 0.20.2.

We first solve the optimization problem (26) using standard gradient descent (GD) and NAG methods, with step-size and momentum chosen according to the optimal tuning for quadratic functions given in [15, Proposition 1]. Using the latter method as a baseline optimizer, we then train a two-layer LSTM to learn an augmentation term vtv_{t} improving the resulting empirical performance (23) over a dataset train\mathcal{F}_{\text{train}} of 10241024 realizations of AA and bb. Specifically, we pick the algorithmic cost function in (22) as 𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(A,b,𝐱)=t=0T|Axtb|2\mathtt{AlgoCost}(A,b,\mathbf{x})=\sum_{t=0}^{T}~|Ax_{t}-b|^{2} for T=10000T=10000, and perform meta-optimization using Adam with a learning rate of 10310^{-3} for 100100 epochs.

We then construct a test dataset by sampling 256256 independent realizations of AA and bb and compare the average-case performance of our learned optimizer against standard methods in solving (1) in Figure 1. As expected, the introduction of a momentum term enables NAG to converge significantly faster than standard gradient descent. Nevertheless, both methods still require a large number of iterations to reach solutions with high accuracy. Remarkably, we observe that our linearly convergent L2O method learns to initially follow the direction of the positive gradient, therefore initially increasing the cost function rather than decreasing it. This behavior does not pertain to classical optimizers and effectively accelerates the accumulation of momentum in the early stages. As demonstrated by Figure 1, this learned behavior results in improved transient performance, without affecting the asymptotic convergence rate.

Refer to caption
Figure 1: Comparison between the average-case performance of classical and learned optimizers in solving the linear regression problem (26); shaded areas and solid lines denote standard deviations and mean values, respectively.

Finally, to assess the impact of our linear convergence guarantees, we repeat the training of a 2-layer LSTM optimizer without enforcing an exponential decay of the augmentations vtv_{t} in (5). We observe that, without clipping the magnitude of vtv_{t} as per (25), the training procedure becomes highly sensitive to the initial parameter choice θ0\theta_{0}. In fact, since the matrix AAA^{\top}A is ill-conditioned, the algorithmic cost function 𝙰𝚕𝚐𝚘𝙲𝚘𝚜𝚝(A,b,𝐱)\mathtt{AlgoCost}(A,b,\mathbf{x}) dramatically increases if the algorithm updates are poorly chosen, causing numerical instability problems when performing backpropagation to update the parameters θ\theta. As a result, the lack of formal guarantees makes training algorithms to perform a large number of optimization steps numerically challenging. In this example, we thus consider a number of optimization steps of T=30T=30. To make things worse, Figure 2 showcase how learned optimizers may diverge when unrolled for more optimization steps than they were during training—even when training and test problems are drawn from the same distribution. These observations highlight that restricting our search space to linearly convergent algorithms yield not only theoretical but also practical computational advantages.

Refer to caption
Figure 2: Comparison between the average-case performance of classical and learned optimizers in solving the linear regression problem (26); without convergence guarantees, learned optimizers may diverge when unrolled for more steps than they were during training (vertical dashed line).

4.3 Augmenting projected gradient descent for MPC

In our second example, we consider a model predictive control setting where only a limited number of optimization steps can be executed in real-time. In particular, we study a discrete-time linear dynamical system described by the state space equation

xt+1=Axt+But+wt=[1101]xt+[01]ut+wt,x_{t+1}=Ax_{t}+Bu_{t}+w_{t}=\begin{bmatrix}1&1\\ 0&1\end{bmatrix}x_{t}+\begin{bmatrix}0\\ 1\end{bmatrix}u_{t}+w_{t}\,, (27)

where xt2x_{t}\in\mathbb{R}^{2} is the system state, utu_{t}\in\mathbb{R} is the control input, and wt2w_{t}\in\mathbb{R}^{2} represents a zero-mean Gaussian process noise term.

The goal at each time step is to solve the following finite-horizon linear quadratic control problem:

minu0,,uT1\displaystyle~\min_{u_{0},\dots,u_{T-1}}~ t=0T1xtQxt+utRut+xTQTxT\displaystyle\sum_{t=0}^{T-1}x_{t}^{\top}Qx_{t}+u_{t}^{\top}Ru_{t}+x_{T}^{\top}Q_{T}x_{T} (28a)
subjectto\displaystyle\operatorname*{subject~to}~ xt+1=Axt+But,\displaystyle x_{t+1}=Ax_{t}+Bu_{t}\,, (28b)
ut𝒰t,xt𝒳,x02;\displaystyle u_{t}\in\mathcal{U}_{t}\,,x_{t}\in\mathcal{X}\,,x_{0}\in\mathbb{R}^{2}\,; (28c)

for simplicity, we assume that the weighting matrices Q,RQ,R and QTQ_{T} are identity matrices of appropriate dimensions. We further set 𝒰t={u:u0.25}\mathcal{U}_{t}=\{u:||u||_{\infty}\leq 0.25\} and 𝒳t=2\mathcal{X}_{t}=\mathbb{R}^{2} to account for actuation constraints yet sidestep recursive feasibility issues for ease of exposition, as our focus lies in efficiently solving the underlying quadratic program.

By introducing the stacked notation 𝐮=[u0uT1]\mathbf{u}=[u_{0}^{\top}\dots u_{T-1}^{\top}]^{\top} and 𝐱=[x0xT]\mathbf{x}=[x_{0}^{\top}\dots x_{T}^{\top}]^{\top}, (28) can be equivalently rewritten as:

min𝐮𝐮(𝐆𝐐𝐆+𝐑)𝐮+2x0𝐅𝐐𝐆𝐮+x0𝐅𝐐𝐅x0\displaystyle~\min_{\mathbf{u}}~\mathbf{u}^{\top}(\mathbf{G}^{\top}\mathbf{Q}\mathbf{G}+\mathbf{R})\mathbf{u}+2x_{0}^{\top}\mathbf{F}^{\top}\mathbf{Q}\mathbf{G}\mathbf{u}+x_{0}^{\top}\mathbf{F}^{\top}\mathbf{Q}\mathbf{F}x_{0} (29a)
subjectto𝐱=𝐅x0+𝐆𝐮,𝐮𝓤,x02,\displaystyle\operatorname*{subject~to}~\mathbf{x}=\mathbf{F}x_{0}+\mathbf{G}\mathbf{u},\mathbf{u}\in\bm{\mathcal{U}},x_{0}\in\mathbb{R}^{2}\,, (29b)

where 𝐐=blkdiag(ITQ,QT)\mathbf{Q}=\operatorname{blkdiag}(I_{T}\otimes Q,Q_{T}), 𝐑=ITR\mathbf{R}=I_{T}\otimes R, and 𝐅\mathbf{F} and 𝐆\mathbf{G} are block matrices comprising the system Markov parameters to encode the system dynamics over a prediction horizon of length T=20T=20.

To solve the quadratic program above, we start from an initial guess 𝐮(0)\mathbf{u}^{(0)} equal to zero and employ the projected gradient descent (PGD) method (16) with step size η=1λmax(𝐆𝐐𝐆+𝐑)3.8105\eta=\frac{1}{\lambda_{\operatorname{max}}(\mathbf{G}^{\top}\mathbf{Q}\mathbf{G}+\mathbf{R})}\approx 3.8\cdot 10^{-5}. Using this as our baseline optimization algorithm, we then learn an augmentation term vtv_{t} parametrized as per (25) to minimize the total average cost (29a) over 100100 optimization steps 𝐮(1),,𝐮(100)\mathbf{u}^{(1)},\dots,\mathbf{u}^{(100)} when each component of the initial state x0x_{0} of the system (27) is drawn from a zero-mean Gaussian random variable with standard deviation of 0.50.5. Specifically, we perform meta-optimization using Adam with a learning rate of 51035\cdot 10^{-3} for 6565 epochs.

Refer to caption
Figure 3: Comparison between the average-case performance of classical and learned optimizers in solving the quadratic program (29); solid lines represent mean values, and shaded areas cover values up to the 90%90\% percentile.

As highlighted by Figure 3, while both algorithms converge in a small neighborhood of the solution of (29) after 100100 iterations, our learned optimizer showcase improved transient performance. Numerically, we verify that this is because during training vtv_{t} learns to promptly saturate the constraint ut𝒰tu_{t}\in\mathcal{U}_{t} when needed, leading to better performance in the testing phase. Remarkably, as shown in Figure 4, we observe a similar trend even when we close the loop between the linear dynamical system (27) and our learned optimizer. In particular, despite the repeated solution of (29) in a receding horizon fashion induces a new distribution over the state space, using our learned optimizer—trained on initial conditions x0x_{0} randomly drawn from a Gaussian distribution—enables a significant reduction of the closed-loop cost incurred over an horizon of length N=30N=30 when only a limited optimization steps can be performed in real-time.

Refer to caption
Figure 4: Comparison between the closed-loop cost incurred by the receding horizon control law that results from unrolling classical and learned optimizer to approximate the solution of the quadratic program (29); solid lines represent mean values, and shaded areas cover values up to the 90%90\% percentile.

5 Conclusion

We have shown that many classical linearly convergent algorithms—ranging from standard gradient descent to accelerated and proximal methods—can be enriched with an exponentially-decaying augmentation term whose role is to improve average‐case performance on a specified distribution of problems. Crucially, this can be achieved without sacrificing their provable worst‐case linear rates. Specifically, (i) these perturbations can be injected at any desired frequency below a threshold, trading off average-case improvement against worst-case rate degradation in a quantifiable way; and (ii) every regular linearly convergent method admits this form of decomposition. In practice, these results enable average-case improvement of baseline algorithms for composite optimization, by designing an exponentially decaying update function akin to designing a robustly stabilizing feedback policy for nonlinear control systems. Numerical examples confirm the potential for significant reduction in the number of optimization steps over general-purpose solvers.

While we have focused on application to learned optimization and average-case improvement, the characterization of all linearly convergent algorithms has independent interest. Accordingly, one direction of interest is to exploit the characterization to deriving update rules with optimal performance from the lens of typical control costs, drawing novel parallels between optimal control theory and accelerated optimization. Important venues of further investigation include analyzing the computational overhead arising from computing learned augmentations, guarantees for derivative-free optimization, design of nonexpansive fixed-point iterations, and applications to time-varying and online optimization.

References

  • [1] Shmuel Agmon. The relaxation method for linear inequalities. Canadian Journal of Mathematics, 6:382–392, 1954.
  • [2] Pierre Alquier. User-friendly introduction to pac-bayes bounds. arXiv preprint arXiv:2110.11216, 2021.
  • [3] Brian DO Anderson. From youla–kucera to identification, adaptive and nonlinear control. Automatica, 34(12):1485–1506, 1998.
  • [4] Marcin Andrychowicz, Misha Denil, Sergio Gomez, Matthew W Hoffman, David Pfau, Tom Schaul, Brendan Shillingford, and Nando De Freitas. Learning to learn by gradient descent by gradient descent. Advances in neural information processing systems, 29, 2016.
  • [5] Sebastian Banert, Jevgenija Rudzusika, Ozan Öktem, and Jonas Adler. Accelerated forward-backward optimization using deep learning. SIAM Journal on Optimization, 34(2):1236–1263, 2024.
  • [6] Leilei Cui, Zhong-Ping Jiang, Eduardo D Sontag, and Richard D Braatz. Perturbed gradient descent algorithms are small-disturbance input-to-state stable. arXiv preprint arXiv:2507.02131, 2025.
  • [7] Timothy A Davis and Yifan Hu. The university of florida sparse matrix collection. ACM Transactions on Mathematical Software (TOMS), 38(1):1–25, 2011.
  • [8] Clara Lucía Galimberti, Luca Furieri, and Giancarlo Ferrari-Trecate. Parametrizations of all stable closed-loop responses: From theory to neural network control design. Annual Reviews in Control, 60:101012, 2025.
  • [9] Guillaume Garrigos and Robert M Gower. Handbook of convergence theorems for (stochastic) gradient methods. arXiv preprint arXiv:2301.11235, 2023.
  • [10] Euhanna Ghadimi, Hamid Reza Feyzmahdavian, and Mikael Johansson. Global convergence of the heavy-ball method for convex optimization. In 2015 European control conference (ECC), pages 310–315. IEEE, 2015.
  • [11] Howard Heaton, Xiaohan Chen, Zhangyang Wang, and Wotao Yin. Safeguarded learned convex optimization. In Proceedings of the AAAI Conference on Artificial Intelligence, volume 37, pages 7848–7855, 2023.
  • [12] Jeffrey Ichnowski, Paras Jain, Bartolomeo Stellato, Goran Banjac, Michael Luo, Francesco Borrelli, Joseph E Gonzalez, Ion Stoica, and Ken Goldberg. Accelerating quadratic optimization with reinforcement learning. Advances in Neural Information Processing Systems, 34:21043–21055, 2021.
  • [13] Hamed Karimi, Julie Nutini, and Mark Schmidt. Linear convergence of gradient and proximal-gradient methods under the polyak-łojasiewicz condition. In Joint European conference on machine learning and knowledge discovery in databases, pages 795–811. Springer, 2016.
  • [14] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014.
  • [15] Laurent Lessard, Benjamin Recht, and Andrew Packard. Analysis and design of optimization algorithms via integral quadratic constraints. SIAM Journal on Optimization, 26(1):57–95, 2016.
  • [16] Ke Li and Jitendra Malik. Learning to optimize. In International Conference on Learning Representations, 2017.
  • [17] Feng-Yi Liao, Lijun Ding, and Yang Zheng. Error bounds, PL condition, and quadratic growth for weakly convex functions, and linear convergences of proximal point methods. In 6th Annual Learning for Dynamics & Control Conference, pages 993–1005. PMLR, 2024.
  • [18] Ian R Manchester, Ruigang Wang, and Nicholas H Barbara. Neural networks in the loop: Learning with stability and robustness guarantees. Annual Review of Control, Robotics, and Autonomous Systems, in press, 2026.
  • [19] Andrea Martin and Luca Furieri. Learning to optimize with convergence guarantees using nonlinear system theory. IEEE Control Systems Letters, 8:1355–1360, 2024.
  • [20] David Mayne. Nonlinear model predictive control: Challenges and opportunities. Nonlinear model predictive control, pages 23–44, 2000.
  • [21] Hesameddin Mohammadi, Meisam Razaviyayn, and Mihailo R Jovanović. Robustness of accelerated first-order algorithms for strongly convex optimization problems. IEEE Transactions on Automatic Control, 66(6):2480–2495, 2020.
  • [22] Yurii Evgen’evich Nesterov. A method of solving a convex programming problem with convergence rate O(1k2){O}\bigl(\frac{1}{k^{2}}\bigr). In Doklady Akademii Nauk, volume 269, pages 543–547. Russian Academy of Sciences, 1983.
  • [23] Boris T Polyak. Some methods of speeding up the convergence of iteration methods. Ussr computational mathematics and mathematical physics, 4(5):1–17, 1964.
  • [24] Ernest K Ryu and Stephen Boyd. Primer on monotone operator methods. Appl. comput. math, 15(1):3–43, 2016.
  • [25] Rajiv Sambharya, Georgina Hall, Brandon Amos, and Bartolomeo Stellato. Learning to warm-start fixed-point optimization algorithms. Journal of Machine Learning Research, 25(166):1–46, 2024.
  • [26] Rajiv Sambharya and Bartolomeo Stellato. Data-driven performance guarantees for classical and learned optimizers. arXiv preprint arXiv:2404.13831, 2024.
  • [27] Carsten Scherer and Christian Ebenbauer. Convex synthesis of accelerated gradient algorithms. SIAM Journal on Control and Optimization, 59(6):4615–4645, 2021.
  • [28] Arjan J Van der Schaft and AJ Van Der Schaft. L2-gain and passivity techniques in nonlinear control, volume 2. Springer, 2000.
  • [29] Bryan Van Scoy, Randy A Freeman, and Kevin M Lynch. The fastest known globally convergent first-order method for minimizing strongly convex functions. IEEE Control Systems Letters, 2(1):49–54, 2017.
  • [30] Dante Youla, Hamid Jabr, and Jr Bongiorno. Modern Wiener-Hopf design of optimal controllers–Part II: The multivariable case. IEEE Transactions on Automatic Control, 21(3):319–338, 1976.
  • [31] Hui Zhang and Wotao Yin. Gradient methods for convex minimization: better rates under weaker conditions. arXiv preprint arXiv:1303.4645, 2013.

Appendix A Proof of Theorem 1

We first prove the result by assuming that π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma). This is instrumental towards establishing the general result. Let δt=dist(ξt,𝙵𝚒𝚡π)\delta_{t}=\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi}) for compactness. By the algorithm definition (5) and the triangle inequality, we have that for every FF\in\mathcal{F}

δt\displaystyle\delta_{t} =dist(π(F,ξt1)+vt1,𝙵𝚒𝚡π)\displaystyle=\operatorname{dist}(\pi(F,\xi_{t-1})+v_{t-1},\mathtt{Fix}_{\pi})
=infc𝙵𝚒𝚡πdist(π(F,ξt1)+vt1,c)\displaystyle=\inf_{c\in\mathtt{Fix}_{\pi}}~\operatorname{dist}(\pi(F,\xi_{t-1})+v_{t-1},c)
infc𝙵𝚒𝚡πdist(π(F,ξt1),c)+|vt1|\displaystyle\leq\inf_{c\in\mathtt{Fix}_{\pi}}~\operatorname{dist}(\pi(F,\xi_{t-1}),c)+|v_{t-1}|
=dist(π(F,ξt1),𝙵𝚒𝚡π)+|vt1|.\displaystyle=\operatorname{dist}(\pi(F,\xi_{t-1}),\mathtt{Fix}_{\pi})+|v_{t-1}|\,.

Assuming that π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma), we have that (4) holds with p(t)=1p(t)=1. It follows that δtγδt1+|vt1|\delta_{t}\leq\gamma\delta_{t-1}+|v_{t-1}|. Iterating this inequality, we deduce that

δt\displaystyle\delta_{t} γtδ0+k=0t1γk|vt1k|\displaystyle\leq\gamma^{t}\delta_{0}+\sum_{k=0}^{t-1}\gamma^{k}|v_{t-1-k}|
γtδ0+k=0t1γkp(t1k)γt1k\displaystyle\leq\gamma^{t}\delta_{0}+\sum_{k=0}^{t-1}\gamma^{k}p(t-1-k)\gamma^{t-1-k}
γt(δ0+1γk=0t1p(k)),\displaystyle\leq\gamma^{t}\left(\delta_{0}+\frac{1}{\gamma}\sum_{k=0}^{t-1}p(k)\right)\,,

where we used the fact that 𝐯exp(m,γ)\mathbf{v}\in\ell_{exp}(m,\gamma). Let q(t)=k=0tp(k)q(t)=\sum_{k=0}^{t}p(k) and note that the right-hand side of the above can be written as γtr(t)\gamma^{t}r(t) where r(t)=δ0+1γq(t1)r(t)=\delta_{0}+\frac{1}{\gamma}q(t-1). We study q(t)q(t). By linearity of summation, q(t)q(t) can be equivalently rewritten as

k=0tj=0majkj=amk=0tkm++a1k=0tk+a0k=0t1,\sum_{k=0}^{t}\sum_{j=0}^{m}a_{j}k^{j}=a_{m}\sum_{k=0}^{t}k^{m}+\dots+a_{1}\sum_{k=0}^{t}k+a_{0}\sum_{k=0}^{t}1\,,

where aja_{j}\in\mathbb{R} with j{0,,m}j\in\{0,\dots,m\} is the jj-th coefficient of the polynomial p()p(\cdot). Faulhaber’s formula implies that q(t)q(t) is a polynomial of degree m+1m+1 in the variable tt with coefficient bm+1=amm+1b_{m+1}=\frac{a_{m}}{m+1}. Furthermore, q(t)q(t) is positive and monotonically non-decreasing by construction, that is, q(t)𝒫m+1(t)q(t)\in\mathcal{P}_{m+1}(t). Note that q(t)𝒫m+1(t)q(t)\in\mathcal{P}_{m+1}(t) implies r(t)𝒫m+1(t)r(t)\in\mathcal{P}_{m+1}(t). Hence, we conclude that δtr(t)γt\delta_{t}\leq r(t)\gamma^{t} for all tt\in\mathbb{N}, which proves the result for the case π𝙴𝚡𝚙π(γ)\pi\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\gamma).

We now turn our attention to the general case where π\pi is any linearly convergent algorithm in 𝚙𝙴𝚡𝚙π(m,γ\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma) certified by a polynomial p(t)𝒫m(t)p(t)\in\mathcal{P}_{m}(t). We distinguish two cases. In the case p(1)γ<1p(1)\gamma<1, we have that π𝚙𝙴𝚡𝚙π(m,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma) also lies in 𝙴𝚡𝚙π(p(1)γ)\mathtt{Exp}^{\pi}_{\mathcal{F}}(p(1)\gamma), and the analysis above applies with rate p(1)γp(1)\gamma in place of γ\gamma.

In the case p(1)γ1p(1)\gamma\geq 1, we necessarily have than N2N\geq 2. For any wtexp(m,γ)w_{t}\in\ell_{exp}(m,\gamma), consider the recursion

ζk+1=πN(F,ζk)+wt,\zeta_{k+1}=\pi^{N}(F,\zeta_{k})+w_{t}\,, (30)

where πN()\pi^{N}(\cdot) denotes the repeated application of π\pi defined as πN(F,ζk)=π(F,πN1(F,ζk))\pi^{N}(F,\zeta_{k})=\pi(F,\pi^{N-1}(F,\zeta_{k})) with π1(F,ζk)=π(F,ζk)\pi^{1}(F,\zeta_{k})=\pi(F,\zeta_{k}).

We first observe that, if ζ0=ξ0\zeta_{0}=\xi_{0} and 𝐯\mathbf{v} is constructed as per (6), then (30) is equivalent to (5) in the sense that ζk=ξNk\zeta_{k}=\xi_{Nk} for every kk\in\mathbb{N}. By construction, πN𝙴𝚡𝚙π(ρ)\pi^{N}\in\mathtt{Exp}^{\pi}_{\mathcal{F}}(\rho) and therefore complies with (4) with p(t)=1p(t)=1 and γ=ρ\gamma=\rho. Hence, as proven above, it holds that

dist(𝜻,𝙵𝚒𝚡πN)exp(m+1,ρ).\operatorname{dist}(\bm{\zeta},\mathtt{Fix}_{\pi^{N}})\in\ell_{exp}(m+1,\rho)\,.

We now argue that 𝙵𝚒𝚡πN=𝙵𝚒𝚡π\mathtt{Fix}_{\pi^{N}}=\mathtt{Fix}_{\pi}. Clearly, 𝙵𝚒𝚡πN𝙵𝚒𝚡π\mathtt{Fix}_{\pi^{N}}\supseteq\mathtt{Fix}_{\pi} since ξ=π(F,ξ)\xi^{\star}=\pi(F,\xi^{\star}) for every ξ𝙵𝚒𝚡π\xi^{\star}\in\mathtt{Fix}_{\pi} and thus ξ=πN(F,ξ)\xi^{\star}=\pi^{N}(F,\xi^{\star}). To show that 𝙵𝚒𝚡πN𝙵𝚒𝚡π\mathtt{Fix}_{\pi^{N}}\subseteq\mathtt{Fix}_{\pi}, assume there exists ζ𝙵𝚒𝚡πN\zeta^{\star}\in\mathtt{Fix}_{\pi^{N}} such that ζ𝙵𝚒𝚡π\zeta^{\star}\notin\mathtt{Fix}_{\pi}. Since π𝚙𝙴𝚡𝚙π(m,γ\pi\in\mathtt{pExp}^{\pi}_{\mathcal{F}}(m,\gamma), we have that limtdist(πt(F,ζ),𝙵𝚒𝚡π)=0\lim_{t\to\infty}~\operatorname{dist}(\pi^{t}(F,\zeta^{\star}),\mathtt{Fix}_{\pi})=0. At the same time, dist(πτN(F,ζ),𝙵𝚒𝚡π)>0\operatorname{dist}(\pi^{\tau N}(F,\zeta^{\star}),\mathtt{Fix}_{\pi})>0 for any τ\tau\in\mathbb{N} because πτN(F,ζ)=ζ𝙵𝚒𝚡π\pi^{\tau N}(F,\zeta^{\star})=\zeta^{\star}\notin\mathtt{Fix}_{\pi}. This is a contradiction, and thus 𝙵𝚒𝚡πN=𝙵𝚒𝚡π\mathtt{Fix}_{\pi^{N}}=\mathtt{Fix}_{\pi}. We conclude that

dist(𝜻,𝙵𝚒𝚡π)exp(m+1,ρ),\operatorname{dist}(\bm{\zeta},\mathtt{Fix}_{\pi})\in\ell_{exp}(m+1,\rho)\,,

and therefore there exists a polynomial q(k)𝒫m+1(k)q(k)\in\mathcal{P}_{m+1}(k) such that dist(ζk,𝙵𝚒𝚡π)q(k)ρk\operatorname{dist}(\zeta_{k},\mathtt{Fix}_{\pi})\leq q(k)\rho^{k} for all kk\in\mathbb{N}.

Next, we note that, for any s{1,,N1}s\in\{1,\dots,N-1\}

dist(ξNk+s,𝙵𝚒𝚡π)\displaystyle\operatorname{dist}(\xi_{Nk+s},\mathtt{Fix}_{\pi}) =dist(πs(F,ξNk),𝙵𝚒𝚡π)\displaystyle=\operatorname{dist}(\pi^{s}(F,\xi_{Nk}),\mathtt{Fix}_{\pi})
p(s)γsdist(ξNk,𝙵𝚒𝚡π)\displaystyle\leq p(s)\gamma^{s}\operatorname{dist}(\xi_{Nk},\mathtt{Fix}_{\pi})
p(s)γsq(k)ρk,\displaystyle\leq p(s)\gamma^{s}q(k)\rho^{k}\,,

where we used the fact that dist(ξNk,𝙵𝚒𝚡π)=dist(ζk,𝙵𝚒𝚡π)\operatorname{dist}(\xi_{Nk},\mathtt{Fix}_{\pi})=\operatorname{dist}(\zeta_{k},\mathtt{Fix}_{\pi}) for any kk\in\mathbb{N}. Letting t=Nk+st=Nk+s, and using the fact that p(),q()𝒫m+1p(\cdot),q(\cdot)\in\mathcal{P}_{m+1}, we obtain

dist(ξt,𝙵𝚒𝚡π)\displaystyle\operatorname{dist}(\xi_{t},\mathtt{Fix}_{\pi}) p(N1)γq(tsN)ρtsN\displaystyle\leq p(N-1)\gamma q\left(\left\lfloor\frac{t-s}{N}\right\rfloor\right)\rho^{\left\lfloor\frac{t-s}{N}\right\rfloor}
p(N1)γq(tN)ρtN+1N\displaystyle\leq p(N-1)\gamma q\left(\frac{t}{N}\right)\rho^{\left\lfloor\frac{t-N+1}{N}\right\rfloor}
p(N1)γρ21Nq(tN)r(t)𝒫m+1(t)(ρ1N)t.\displaystyle\leq\underbrace{\frac{p(N-1)\gamma}{\rho^{2-\frac{1}{N}}}q\left(\frac{t}{N}\right)}_{r(t)\in\mathcal{P}_{m+1}(t)}\left(\rho^{\frac{1}{N}}\right)^{t}\,.

Since ρ=p(N)γN\rho=p(N)\gamma^{N}, we have that ρ1N=p(N)Nγ\rho^{\frac{1}{N}}=\sqrt[N]{p(N)}\gamma. This concludes the proof.

Appendix B Proof of Theorem 2

Let vt(F,ξt:0)=π(F,χt)+σt(F,χt:0)v_{t}(F,\xi_{t:0})=-\pi(F,\chi_{t})+\sigma_{t}(F,\chi_{t:0}). We first show by induction that ξt=χt\xi_{t}=\chi_{t} at all times, starting from the base case ξ0=χ0\xi_{0}=\chi_{0}, which holds by construction. Assume now that ξt:0=χt:0\xi_{t:0}=\chi_{t:0}. We aim to prove that ξt+1=χt+1\xi_{t+1}=\chi_{t+1}. This holds because

ξt+1\displaystyle\xi_{t+1} =π(F,ξt)π(F,χt)+σt(F,χt:0)\displaystyle=\pi(F,\xi_{t})-\pi(F,\chi_{t})+\sigma_{t}(F,\chi_{t:0})
=σt(F,χt:0)=χt+1.\displaystyle=\sigma_{t}(F,\chi_{t:0})=\chi_{t+1}\,.

It remains to show that the sequence vt(F,ξt:0)=π(F,χt)+σt(χt:0)v_{t}(F,\xi_{t:0})=-\pi(F,\chi_{t})+\sigma_{t}(\chi_{t:0}) belongs to exp(m,γ)\ell_{exp}(m,\gamma). To prove this, we rewrite vtv_{t} as

vt=(π(F,χt)χt)+σt(F,χt:0)χt.\displaystyle v_{t}=-(\pi(F,\chi_{t})-\chi_{t})+\sigma_{t}(F,\chi_{t:0})-\chi_{t}\,. (31)

Since π(F,)\pi(F,\cdot) is Lipschitz continuous, letting χtp\chi_{t}^{p} be any element of argminχ𝙵𝚒𝚡π|χχt|2\operatorname*{arg\,min}_{\chi\in\mathtt{Fix}_{\pi}}|\chi-\chi_{t}|^{2}, there exists a constant Lπ+L_{\pi}\in\mathbb{R}_{+} such that

|π(F,χt)χt|\displaystyle|\pi(F,\chi_{t})-\chi_{t}| =|π(F,χt)χtp+χtpχt|\displaystyle=|\pi(F,\chi_{t})-\chi_{t}^{p}+\chi_{t}^{p}-\chi_{t}|
=|π(F,χt)π(F,χtp)+χtpχt|\displaystyle=|\pi(F,\chi_{t})-\pi(F,\chi_{t}^{p})+\chi_{t}^{p}-\chi_{t}|
(Lπ+1)|χtχtp|\displaystyle\leq(L_{\pi}+1)|\chi_{t}-\chi_{t}^{p}|
=(Lπ+1)dist(χt,𝙵𝚒𝚡π),\displaystyle=(L_{\pi}+1)\operatorname{dist}(\chi_{t},\mathtt{Fix}_{\pi})\,,

and hence (π(F,𝝌)𝝌)exp(m,γ)-(\pi(F,\bm{\chi})-\bm{\chi})\in\ell_{exp}(m,\gamma). We further have that 𝝈(F,𝝌)𝝌exp(m,γ)\bm{\sigma}(F,\bm{\chi})-\bm{\chi}\in\ell_{exp}(m,\gamma) by the regularity assumption on 𝝈\bm{\sigma} as per Definition 2. Since the sum of signals in exp(m,γ)\ell_{exp}(m,\gamma) belongs to exp(m,γ)\ell_{exp}(m,\gamma), we conclude the proof by inspection of (31).