Computer-aided Lyapunov analyses in optimization

Many convergence proofs for first-order methods in optimization follow a common pattern: construct a Lyapunov function, establish an inequality describing its decrease, and use this inequality to show that a performance measure converges to zero. Within many standard problem classes and prescribed families of Lyapunov functions, these tasks can often be formulated exactly as semidefinite programs (SDPs). Existing solvers can then be used to search numerically for feasible solutions. Convergence analysis itself thus becomes a tractable optimization problem.

Closely related is the performance estimation problem (PEP) methodology, introduced by Drori and Teboulle [DT14]. For many first-order methods and problem classes, the search for a worst-case example can be reformulated as an SDP, giving tight worst-case guarantees [THG17]. The standard PEP methodology is now a mature technology for worst-case analysis, supported by software packages such as PEPit [GMG+24].

However, standard PEP formulations analyze a fixed number of iterations, with the number of variables and constraints typically growing quadratically with this horizon. Deriving closed-form guarantees valid for all iterations can also be challenging. Automated Lyapunov analyses provide a complementary approach: restricting the search to prescribed families of Lyapunov analyses can give smaller SDPs whose number of variables and constraints is independent of the iteration horizon, or grows only linearly with it, facilitating concise, closed-form convergence proofs [TVSL18, TB19, UBTG25, UDGT+26]. A closely related framework uses integral quadratic constraints (IQCs), developed in robust control [MR97], to obtain convergence guarantees through linear matrix inequalities (LMIs) that can be checked by semidefinite programming [LRP16, MSE21]. Although restricting the Lyapunov search may introduce conservatism, tight worst-case guarantees can still be recovered in many settings.

In this post, we derive the corresponding SDPs explicitly for two examples: an inexact gradient method, for which the Lyapunov function is fixed in advance, and a gradient method with memory, for which we search for a Lyapunov function within a prescribed family. The aim is to explain the main ideas in a simple setting and provide a starting point for analyses of other methods and problem classes. For the first example, we also show how to extract a closed-form proof of a tight convergence guarantee from numerical SDP solutions. In the second example, we show how allowing a more flexible class of Lyapunov functions can give tighter convergence guarantees. This post assumes familiarity with smooth and strongly convex optimization, basic matrix algebra, and Lagrangian duality. On a first reading, the expandable derivations and proofs can be skipped.

1. Problem setting

For simplicity, both examples use the familiar setting of smooth and strongly convex optimization, although the tools developed here apply more broadly. We consider the unconstrained problem

where \(f:\mathbb{R}^d\to\mathbb{R}\) is differentiable, has an \(L\)-Lipschitz continuous gradient (that is, \(f\) is \(L\)-smooth), and is \(\mu\)-strongly convex, with \(0<\mu<L\). We denote this class of functions by \(\mathcal{F}_{\mu,L}(\mathbb{R}^d)\). Equivalently, \(f\) satisfies

where \(\langle\cdot,\cdot\rangle\) is the Euclidean inner product and \(\lVert\cdot\rVert\) is its induced norm. For the quadratic bounds above, see [Bec17, Theorems 5.8(ii) and 5.24(ii)]. Strong convexity ensures that problem \(\eqref{eq:optimization-problem}\) has a unique minimizer \(x^\star\in\mathbb{R}^{d}\), characterized by \(\nabla f(x^\star)=0\) [Bec17, Theorems 3.63 and 5.25(a)]. Figure 1 illustrates these inequalities geometrically, showing how \(f\) is globally bounded above and below by quadratic functions tangent to it at any chosen point.

2. Fixed Lyapunov functions

Here we check a Lyapunov function chosen in advance, then use the tools developed here to search for one in Section 3.

2.1 The gradient method with inexact gradients

The standard gradient method provides a familiar starting point; here, we consider a slightly richer variant with inexact gradients. Given a constant step size \(\gamma>0\) and an initial point \(x^{0}\in\mathbb{R}^{d}\), consider

where \(\varepsilon\in[0,1)\) bounds the relative error in the gradient approximation. We seek convergence guarantees that hold for every sequence of approximations \(d^k\) satisfying this bound. When \(\varepsilon=0\), we recover the standard gradient method.

2.2 Lyapunov functions and contraction

To quantify convergence, we track a nonnegative scalar quantity such as squared distance to the solution \(\|x^{k}-x^{\star}\|^2\) or function-value suboptimality \(f(x^{k})-f(x^{\star})\). We call such a quantity a performance measure. Under the standing assumptions, it is natural to seek bounds of the form

and to determine which choices of \(\gamma\) ensure that the corresponding factors satisfy

In these examples, we use the performance measure itself as a Lyapunov function, establishing convergence by proving a contraction inequality for that same quantity. More generally, we consider nonnegative and nontrivial Lyapunov functions \(\mathcal{V}\) satisfying

and seek conditions ensuring \(0\leqslant\rho(\gamma,\varepsilon,\mu,L)<1\). Iterating the contraction inequality \eqref{eq:lyapunov-contraction} gives

and we say that \(\mathcal{V}\) converges linearly to zero with rate at most \(\rho(\gamma,\varepsilon,\mu,L)\). In particular, we seek the smallest factor for which \eqref{eq:lyapunov-contraction} holds.

We build on squared distance \(\|x^{k}-x^{\star}\|^2\) and function-value suboptimality \(f(x^{k})-f(x^{\star})\) in \(\eqref{eq:contraction-examples}\) by also using gradient information. Including the squared gradient norm \(\|\nabla f(x^{k})\|^2\) and the inner product \(\langle x^{k}-x^{\star},\nabla f(x^{k})\rangle\) allows us to account for the gradient’s magnitude and its alignment with \(x^{k}-x^{\star}\). We therefore consider the family

where \(Q=(Q_{i,j})_{1\leqslant i,j\leqslant2}\in\mathbb{R}^{2\times2}\) is symmetric and \(q\in\mathbb{R}\).

For fixed \(Q\) and \(q\), \(\mathcal{V}\) is linear in the squared norms, inner products, and function-value differences that become variables of the SDP derived in Section 2.3. More general mixed terms need not admit a linear representation in those variables. For that derivation, we group the two vectors \(x^{k}-x^{\star}\) and \(\nabla f(x^{k})\) and write

where \(\otimes\) is the Kronecker product and \(I_{d}\) is the identity matrix in \(\mathbb{R}^{d\times d}\).

Here, we fix \(Q\) and \(q\) so that \(\mathcal{V}\) is nonnegative and nontrivial. The two cases in \(\eqref{eq:contraction-examples}\) correspond to

2.3 Optimal contraction via semidefinite programming

To verify \(\eqref{eq:lyapunov-contraction}\) along every trajectory, it suffices to consider one step from an arbitrary point \(x^{k}\). The update rule and the coefficients \(Q\) and \(q\) are the same at every iteration, and any such step can be realized as the first step by setting \(x^{0}=x^{k}\). We therefore define the optimal one-step contraction factor by

We will show, step by step, how semidefinite programming gives an upper bound on this quantity and under which conditions the bound equals the optimal one-step contraction factor. We first fix a candidate \(\rho\) and derive a semidefinite feasibility test, then minimize \(\rho\) in this test. This feasibility formulation will also let us search for valid Lyapunov functions in Section 3, where their coefficients become additional variables.

Worst-case Lyapunov difference. For a fixed candidate \(\rho\), we maximize the difference \(\mathcal{V}(k+1)-\rho\mathcal{V}(k)\) over all admissible functions and updates:

The contraction inequality \(\eqref{eq:lyapunov-contraction}\) holds for every admissible function and update exactly when \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\leqslant0\). The optimal one-step contraction factor is therefore

We next derive a semidefinite test for the condition \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\leqslant0\). We will describe the possible function values and gradients using finitely many inequalities, collect vector inner products in a positive semidefinite matrix, and combine these inequalities to certify contraction.

Encoding inner products with a Gram matrix. Before we proceed, we introduce some convenient notation. Let \(\mathbb{S}^{n}\) and \(\mathbb{S}_{+}^{n}\) denote the sets of real-valued symmetric and positive semidefinite \(n\times n\) matrices, respectively. We define the Gramian function \(\mathtt{Gram}:(\mathbb{R}^{d})^{n}\to\mathbb{S}^{n}_{+}\) by

In particular, \((\mathtt{Gram}(\mathbf{z}))_{i,i}=\norm{z_i}^{2}\), so squared norms appear as diagonal entries. With this, we can write the Lyapunov function as

For \(A\in\mathbb{R}^{m\times n}\), \(M\in\mathbb{S}^{m}\), and \(\mathbf{z}\in(\mathbb{R}^{d})^{n}\), we will use the identities

Rewriting the objective. Our goal is to express the objective of \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\) in terms of the function-value differences \(f(x^{k})-f(x^{\star})\) and \(f(x^{k+1})-f(x^{\star})\), together with the Gram matrix \(\mathtt{Gram}(x^{k}-x^{\star},\nabla f(x^{k}),\nabla f(x^{k+1}),d^{k})\). For this purpose, we introduce the following matrices, whose columns correspond, from left to right, to the four arguments of this Gramian:

Applying the trace identity in \(\eqref{eq:gram-trace-identities}\) with these matrices gives the objective in the desired form

Optimizing over functions. Using this trace representation of the objective gives the following equivalent formulation of the problem defining \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\)

Here, \(f\) is an infinite-dimensional optimization variable, even though only its values and gradients at \(x^{k}\), \(x^{k+1}\), and \(x^{\star}\) appear in the objective and constraints. Introducing corresponding variables \(f^i\) and \(g^i\) therefore gives the equivalent formulation

At first glance, this may not seem like progress: we have replaced the infinite-dimensional variable by an existence constraint over the same function class. Have we actually made the problem any easier?

Determining whether the prescribed values and gradients can come from a single function in \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\) is an interpolation problem. For this class, interpolation conditions from the performance estimation literature let us replace the existence constraint by finitely many explicit inequalities [THG17]. The quadratic bounds \(\eqref{eq:quadratic-bounds}\), imposed only at the sampled points, are necessary but not sufficient for interpolation. Theorem 1 gives necessary and sufficient conditions for a function in \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\) to match all the prescribed values and gradients.

For an intuitive explanation of interpolation in the convex, possibly nonsmooth case, corresponding to \(\mu=0\) and \(L=+\infty\), see Adrien Taylor’s excellent post on Francis Bach’s blog [Tay20].

Applying Theorem 1 gives the following equivalent finite-dimensional formulation, where we have squared the relative-error bound and written each inequality with zero on the right

We have now eliminated the unknown function: only finitely many vector and scalar variables remain. For the three indices in \(\{k,k+1,\star\}\), there are six interpolation inequalities: one for each ordered pair of distinct indices. Swapping \(i\) and \(j\) generally gives a different inequality.

Rewriting the constraints. Substituting \(x^{k+1}=x^{k}-\gamma d^{k}\) and \(g^{\star}=0\) eliminates \(x^{k+1}\) and \(g^{\star}\) as optimization variables. The relative-error constraint becomes

For each interpolation inequality, we use \(A_{i,j}\in\mathbb{S}^{4}\) to represent the inner products and squared norms and \(a_{i,j}\in\mathbb{R}^{2}\) to select the function-value difference \(f^{j}-f^{i}\). The six inequalities can then be written as

for each pair of distinct indices \(i,j\in\{k,k+1,\star\}\).

Combining the rewritten objective and constraints gives the equivalent formulation

We have now obtained an equivalent quadratically constrained quadratic program (QCQP) for evaluating \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\), but it is generally nonconvex. The following result lets us replace the vectors by their Gram matrix, with a rank constraint accounting for the dimension \(d\).

Which matrices are Gram matrices? For \(G\in\mathbb{S}^{n}\), we have that

Primal SDP. Writing \(G=\mathtt{Gram}(x^{k}-x^{\star},g^{k},g^{k+1},d^{k})\) and \(F=\begin{bmatrix}f^{k}-f^{\star}\\f^{k+1}-f^{\star}\end{bmatrix}\in\mathbb{R}^{2}\), we obtain from \(\eqref{eq:gram-characterization}\) the equivalent formulation

The objective and scalar inequality constraints are linear in \(G\) and \(F\). Positive semidefiniteness is a convex constraint, so dropping the rank constraint gives the SDP

and hence an upper bound valid in every dimension \(d\)

For \(d\geqslant4\), the rank constraint is automatic, so

The SDP \(\mathcal{P}^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\) has one \(4\times4\) positive semidefinite matrix variable, two scalar variables, and seven scalar inequalities: six interpolation inequalities and one relative-error inequality. Its size is independent of \(d\).

Dual SDP. Next, we derive the Lagrange dual of \(\mathcal{P}^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\). We associate dual multipliers with the constraints as follows

for each pair of distinct indices \(i,j\in\{k,k+1,\star\}\). With the notation \(\lambda=(\lambda_{i,j})_{i,j\in\{k,k+1,\star\},\,i\ne j}\), the Lagrangian is

Taking the supremum of the Lagrangian over \(G\in\mathbb{S}_{+}^{4}\) and \(F\in\mathbb{R}^{2}\) gives the dual objective

Minimizing this dual objective over the nonnegative multipliers gives the dual SDP

Here, solving the dual SDP amounts to finding multipliers satisfying its constraints. The primal relaxation and weak duality (the primal SDP optimal value is no larger than the dual SDP optimal value) [BV04, Section 5.2.2] give

Thus, if we find such multipliers, the right-hand side is zero and \(\mathcal{P}_{d}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\leqslant0\).

Optional: Strong duality. For \(\varepsilon>0\), the primal and dual optimal values agree (known as strong duality),

This follows from the refined Slater condition [BV04, Section 5.2.3]. The verification is technical, so we include it below as optional reading for those interested in the details.

When \(\varepsilon=0\), we have \(d^{k}=g^{k}\), so the four-vector Gram matrix cannot be positive definite. We can use the smaller Gramian \(\mathtt{Gram}(x^{k}-x^{\star},g^{k},g^{k+1})\) to obtain a different primal–dual SDP pair for which strong duality holds. We do not study this reduced formulation further here.

Dual certificates. A feasible dual point gives nonnegative weights that combine the relative-error and interpolation inequalities into a proof of the desired Lyapunov inequality. To make this proof explicit, let \((\tau,\lambda)\) be a feasible point of \(\mathcal{D}^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\). Consider any step of the method \(\eqref{eq:relatively-inexact-gradient-method}\), and set

Since \(G\succcurlyeq0\), the constraints of \(\mathcal{D}^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\) give \(\mathcal{L}(G,F,\tau,\lambda)\leqslant0\). Thus, for this step,

The last inequality uses primal feasibility and nonnegative multipliers. Since \(x^k\) and \(d^k\) were arbitrary subject to the relative-error condition, the same feasible point \((\tau,\lambda)\) certifies \(\mathcal{V}(k+1)\leqslant\rho\mathcal{V}(k)\) at every iteration along every admissible trajectory, in every dimension. Conversely, for \(d\geqslant4\) and \(\varepsilon>0\), if \(\mathcal{V}(k+1)\leqslant\rho\mathcal{V}(k)\) holds for every admissible function and step, then \(\mathcal{D}^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L,\rho)\) has a feasible point \((\tau,\lambda)\). This follows from exactness of the primal formulation and strong duality.

Formal result. With \(Q\) and \(q\) fixed, minimizing \(\rho\) subject to the constraints of \(\mathcal{D}^{\mathrm{SDP}}\) is itself an SDP, since these constraints are affine in \(\rho,\tau,\lambda\). The derivation above gives Theorem 2.

2.4 Numerical results

The numerical results were obtained with MOSEK [Mos26].

We now compute \(\rho^{\mathrm{SDP}}(Q,q,\gamma,\varepsilon,\mu,L)\) in the case of the squared distance to the solution and the function-value suboptimality. Recall from \eqref{eq:lyapunov-examples} that

These two choices give

The resulting one-step guarantees are

Figure 2 compares these factors as the step size \(\gamma\) varies.

At the sampled parameter values, the numerical results suggest that function-value suboptimality gives a one-step contraction factor no larger than that for the squared distance to the solution. We therefore study the function-value case further by solving \(\mathcal{D}^{\mathrm{SDP}}(0,1,\gamma,\varepsilon,\mu,L,\rho)\) with the candidate factor suggested by fitting the numerical data in Figure 2:

Figure 3 shows the multipliers returned by the solver. These need not be unique; we look for simple patterns that can guide an explicit proof.

2.5 From numerical multipliers to an exact proof

Fixing all multipliers to the fitted expressions in Figure 3 preserves feasibility at the checked sample points. We now show directly that, for \(\varepsilon\in(0,1)\) and \(\gamma\in(0,\tfrac{2}{L(1+\varepsilon)})\), the factor

can be certified using the following multipliers, obtained by fitting the numerical data in Figure 3:

Indeed, following the dual-certificate argument in \(\eqref{eq:dual-certificate-bound}\), we combine the relative-error condition and the three interpolation inequalities below, using the multipliers indicated on the right:

The function-value terms in the weighted sum give the difference on the left below. After substituting \(x^{k+1}=x^k-\gamma d^k\) and \(\nabla f(x^\star)=0\), we move the remaining quadratic terms to the right and complete squares to obtain

which establishes \(f(x^{k+1})-f(x^{\star})\leqslant\rho_{\mathtt{func}}(f(x^k)-f(x^{\star}))\). The key point is how nonnegative weights combine known inequalities into a contraction bound, following the argument in \(\eqref{eq:dual-certificate-bound}\); the nonpositive sum of weighted squares in \eqref{eq:function-value-sos-bound} explicitly represents \(\mathcal{L}(G,F,\tau,\lambda)\).

The function-value contraction factor

is minimized when its two terms are equal, at

with minimum value

Tightness. For every \(d\in\mathbb{N}\) and \(\gamma\in(0,\tfrac{2}{L(1+\varepsilon)})\), there exist \(f\in\mathcal{F}_{\mu,L}(\mathbb{R}^d)\) and a step of \(\eqref{eq:relatively-inexact-gradient-method}\) with \(x^k\ne x^\star\) such that

Optional: Exact gradients. For \(\varepsilon=0\) and \(0<\gamma<\tfrac{2}{L}\), we recover the known sharp function-value bound for the gradient method [THG18, Theorem 3.3]:

This follows by letting \(\varepsilon\downarrow0\) in the contraction bound while keeping the exact-gradient step fixed.

Related literature. Gannot [Gan22, Remark 1.6] states \(f(x^{k+1})-f(x^\star)\leqslant\rho_{\mathtt{func}}(\gamma,\varepsilon,\mu,L)(f(x^k)-f(x^\star))\), drawing on a proof by Etienne de Klerk, François Glineur, and Adrien B. Taylor [dKGT20, Appendix B]. A proof for \(\gamma\in(0,\gamma_\star]\) is given by Uschmajew and Vandereycken [UV22, Theorem 4.2].

2.6 Optional: Reducing the certificate to two squared norms

3. Searching for Lyapunov functions

In Section 2, we compared the one-step contraction factors obtained from squared distance to the solution and function-value suboptimality. A contracting Lyapunov function that bounds a performance measure from above can certify a faster asymptotic rate than a direct contraction bound for that measure, or establish convergence even when the measure does not contract at every iteration. Here we search for such a Lyapunov function within a chosen family. We illustrate this search using stationary gradient methods with memory to solve \(\eqref{eq:optimization-problem}\), i.e., to minimize a function \(f\in\mathcal{F}_{\mu,L}(\mathbb{R}^d)\).

3.1 Stationary gradient methods with memory

Given constants \(\gamma>0\) and \(\alpha,\beta\in\mathbb{R}\), and initial points \(x^{0},x^{-1}\in\mathbb{R}^{d}\), consider

The method is called stationary because \(\gamma\), \(\alpha\), and \(\beta\) are constant across iterations. In the table below, we give some standard choices of \(\gamma\), \(\alpha\), and \(\beta\).

1. Gradient method: Optimal constant step size for worst-case distance contraction on \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\) [Pol63, Theorem 6; Nes18, Theorem 2.1.15].
2. Polyak's heavy-ball method: Optimal parameters for worst-case asymptotic convergence on quadratics in \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\) [Pol64, Theorem 9(3) and the note on p. 15].
3. Nesterov's method: Optimal dimension-independent decay exponent for convergence in distance on \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\), up to constant factors [Nes18, Theorem 2.2.1 and Equation (2.2.22); DT22, Corollary 3].
4. Triple momentum method [VSFL18, Definition 2]: Optimal dimension-independent decay exponent for convergence in distance on \(\mathcal{F}_{\mu,L}(\mathbb{R}^{d})\) [DT22, Corollary 3].

3.2 Choosing a family of Lyapunov functions

To establish a convergence rate for a nonnegative performance measure \(\mathcal{M}\), we look for a Lyapunov function \(\mathcal{V}\) and a contraction factor \(\rho\in[0,1)\) satisfying

for every \(k\geqslant0\), every \(f\in\mathcal{F}_{\mu,L}(\mathbb{R}^d)\), and every trajectory of \(\eqref{eq:gradient-method-with-memory}\).

Condition \(\eqref{eq:memory-lyapunov-contraction}\) requires the Lyapunov function to contract, but imposing it alone would allow the trivial choice \(\mathcal{V}\equiv0\), regardless of \(\mathcal{M}\). Condition \(\eqref{eq:memory-lyapunov-lower-bound}\) connects the Lyapunov function to the performance measure by requiring \(\mathcal{V}\) to bound \(\mathcal{M}\) from above. Iterating \(\eqref{eq:memory-lyapunov-contraction}\) and applying \(\eqref{eq:memory-lyapunov-lower-bound}\) then gives, for \(k\geqslant0\),

Thus, the performance measure \(\mathcal{M}\) converges linearly to zero with rate at most \(\rho\). We can therefore keep the performance measure \(\mathcal{M}\) fixed and search for a contracting upper bound \(\mathcal{V}\).

We restrict \(\mathcal{V}\) and \(\mathcal{M}\) to the following forms:

where \(Q,P\in\mathbb{S}^3\) and \(q,p\in\mathbb{R}\). We fix \(P,p\) to specify a nonnegative performance measure, i.e., \(\mathcal{M}(k)\geqslant0\), and search for \(Q,q\) such that \(\eqref{eq:memory-lyapunov-contraction}\) and \(\eqref{eq:memory-lyapunov-lower-bound}\) hold.

For example, this family includes

3.3 Optimal contraction via semidefinite programming

Adapting the argument of Section 2.3 separately to \(\eqref{eq:memory-lyapunov-contraction}\) and the inequality \(\mathcal{V}(k)\geqslant\mathcal{M}(k)\) in \(\eqref{eq:memory-lyapunov-lower-bound}\) gives Theorem 3. Nonnegativity of \(\mathcal{M}\) is ensured by our choice of \(P,p\). For completeness, we include a proof for the interested reader after the formal statement.

Searching for a Lyapunov function. We now let the parameters \(Q,q\) defining \(\mathcal{V}\) vary and minimize \(\rho\) subject to both tests in Theorem 3. The resulting problem has bilinear terms \(\rho Q\) and \(\rho q\), but becomes an SDP feasibility problem when \(\rho\) is fixed. We denote its optimal value by

Two observations allow us to use bisection with an SDP solver:

1. SDP feasibility. With \(\rho\) fixed, the constraints in \(\eqref{eq:memory-lyapunov-search}\) are affine in \(Q,q\) and the multipliers \(\lambda^{\mathrm{C1}},\lambda^{\mathrm{C2}}\), so we have an SDP feasibility problem.

2. Monotonicity. If \(\eqref{eq:memory-lyapunov-search}\) is feasible at \(\rho\), it is feasible at every \(\rho'\in(\rho,1)\).

   Show proofHide proof

   A feasible point gives \(\mathcal{V}(k)\geqslant\mathcal{M}(k)\geqslant0\). Hence

   \begin{equation*} \mathcal{V}(k+1)\leqslant\rho\mathcal{V}(k)\leqslant\rho'\mathcal{V}(k). \end{equation*}

   Since \(\eqref{eq:memory-lyapunov-search}\) does not depend on \(d\), we can take \(d=4\) and apply the necessity part of Theorem 3 to obtain C1 multipliers \(\lambda^{\mathrm{C1}}\) for \(\rho'\), with the same \(Q,q\). Together with the original C2 multipliers \(\lambda^{\mathrm{C2}}\), these give a feasible point of \(\eqref{eq:memory-lyapunov-search}\) at \(\rho'\).

We can therefore use bisection to search for the smallest contraction factor \(\rho\) certifiable by \(\eqref{eq:memory-lyapunov-search}\), which is also the smallest factor within the class of Lyapunov functions in \(\eqref{eq:memory-lyapunov-ansatz}\) when \(d\geqslant4\).

3.4 Numerical results

The numerical results were obtained with MOSEK [Mos26].

For \(\mathcal{M}(k)=\|x^k-x^\star\|^2\), we now use numerical bisection to estimate the optimal value of \(\eqref{eq:memory-lyapunov-search}\) for the four methods in Section 3.1, using the parameters given there. Figure 4 compares the results with and without the additional constraints \(Q\succcurlyeq0\) and \(q\geqslant0\).

For the parameter choices in Section 3.1, the optimal value \(\rho^\star\) of \(\eqref{eq:memory-lyapunov-search}\) depends only on the ratio \(L/\mu\), which explains the horizontal axis in Figure 4 and allows us to compute the curves with \(L=1\) and \(\mu\in(0,1)\).

Conclusions. Removing the additional constraints \(Q\succcurlyeq0\) and \(q\geqslant0\) can improve the certified contraction factors. For triple momentum at larger condition numbers, the search without these additional constraints still certifies contraction, while no factor below \(1\) is found when they are imposed. The restricted family includes familiar choices such as squared distance to the solution, function-value suboptimality, and squared gradient norm, illustrating the value of searching over a more flexible class of Lyapunov functions.

4. Going further

Both examples use semidefinite programming to combine valid inequalities into a Lyapunov certificate, either for a fixed Lyapunov function or while searching for one within a chosen family.

The theoretical framework in [UBTG25] extends these ideas beyond the contraction bounds considered here. It covers sublinear convergence guarantees and splitting methods, including proximal gradient and Douglas–Rachford.

AutoLyap, described in [UDGT+26], supports broader classes of optimization and inclusion problems, iteration-dependent methods and analyses, and richer families of Lyapunov functions. Users specify the problem class, method, and desired Lyapunov analysis through a high-level interface. The software constructs and numerically solves the corresponding SDPs. Its documentation provides worked examples.

Acknowledgements and use of generative AI. I thank Adrien Taylor and Adrien Banse for their helpful comments and suggestions. Large language models (LLMs) were used to help produce figures, correct typos, and improve wording.

@misc{upadhyaya2026lyapunov,
  author = {Upadhyaya, Manu},
  title = {Computer-aided {Lyapunov} analyses in optimization},
  year = {2026},
  howpublished = {Technical blog post},
  publisher = {Zenodo},
  doi = {10.5281/zenodo.23521024},
  url = {https://manuupadhyaya.github.io/blog/computer-aided-lyapunov-analyses-in-optimization/},
}

References

1. [Bec17]

   Amir Beck. First-Order Methods in Optimization. SIAM, 2017. doi:10.1137/1.9781611974997.

2. [BV04]

   Stephen Boyd and Lieven Vandenberghe. Convex Optimization. Cambridge University Press, 2004. URL.

3. [DT22]

   Yoel Drori and Adrien Taylor. On the oracle complexity of smooth strongly convex minimization. Journal of Complexity, 68:101590, 2022. doi:10.1016/j.jco.2021.101590.

4. [DT14]

   Yoel Drori and Marc Teboulle. Performance of first-order methods for smooth convex minimization: A novel approach. Mathematical Programming, 145(1–2):451–482, 2014. doi:10.1007/s10107-013-0653-0.

5. [Gan22]

   Oran Gannot. A frequency-domain analysis of inexact gradient methods. Mathematical Programming, 194:975–1016, 2022. doi:10.1007/s10107-021-01665-8.

6. [Gen98]

   James E. Gentle. Numerical Linear Algebra for Applications in Statistics. Springer, 1998. doi:10.1007/978-1-4612-0623-1.

7. [GMG+24]

   Baptiste Goujaud, Céline Moucer, François Glineur, Julien M. Hendrickx, Adrien B. Taylor, and Aymeric Dieuleveut. PEPit: Computer-assisted worst-case analyses of first-order optimization methods in Python. Mathematical Programming Computation, 16:337–367, 2024. doi:10.1007/s12532-024-00259-7.

8. [LRP16]

   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. doi:10.1137/15M1009597.

9. [MR97]

   Alexandre Megretski and Anders Rantzer. System analysis via integral quadratic constraints. IEEE Transactions on Automatic Control, 42(6):819–830, 1997. doi:10.1109/9.587335.

10. [MSE21]

    Simon Michalowsky, Carsten Scherer, and Christian Ebenbauer. Robust and structure exploiting optimisation algorithms: An integral quadratic constraint approach. International Journal of Control, 94(11):2956–2979, 2021. doi:10.1080/00207179.2020.1745286.

11. [Nes18]

    Yurii Nesterov. Lectures on Convex Optimization. 2nd ed., Springer, 2018. doi:10.1007/978-3-319-91578-4.

12. [Pol63]

    Boris T. Polyak. Gradient methods for the minimisation of functionals. USSR Computational Mathematics and Mathematical Physics, 3(4):864–878, 1963. doi:10.1016/0041-5553(63)90382-3.

13. [Pol64]

    Boris T. Polyak. Some methods of speeding up the convergence of iteration methods. USSR Computational Mathematics and Mathematical Physics, 4(5):1–17, 1964. doi:10.1016/0041-5553(64)90137-5.

14. [Tay20]

    Adrien Taylor. Computer-aided analyses in optimization. Machine Learning Research Blog, 2020. URL.

15. [TB19]

    Adrien Taylor and Francis Bach. Stochastic first-order methods: Non-asymptotic and computer-aided analyses via potential functions. Proceedings of the Thirty-Second Conference on Learning Theory, 99:2934–2992, 2019. URL.

16. [TVSL18]

    Adrien Taylor, Bryan Van Scoy, and Laurent Lessard. Lyapunov functions for first-order methods: Tight automated convergence guarantees. Proceedings of the 35th International Conference on Machine Learning, 80:4897–4906, 2018. URL.

17. [THG17]

    Adrien B. Taylor, Julien M. Hendrickx, and François Glineur. Smooth strongly convex interpolation and exact worst-case performance of first-order methods. Mathematical Programming, 161:307–345, 2017. doi:10.1007/s10107-016-1009-3.

18. [THG18]

    Adrien B. Taylor, Julien M. Hendrickx, and François Glineur. Exact worst-case convergence rates of the proximal gradient method for composite convex minimization. Journal of Optimization Theory and Applications, 178(2):455–476, 2018. doi:10.1007/s10957-018-1298-1. arXiv:1705.04398v4.

19. [TUG+26]

    Daniel Berg Thomsen, Manu Upadhyaya, Baptiste Goujaud, Aymeric Dieuleveut, and Adrien Taylor. Finding Simple Proofs for First-Order Optimization. 2026. arXiv:2607.08753v1.

20. [UBTG25]

    Manu Upadhyaya, Sebastian Banert, Adrien B. Taylor, and Pontus Giselsson. Automated tight Lyapunov analysis for first-order methods. Mathematical Programming, 209:133–170, 2025. doi:10.1007/s10107-024-02061-8.

21. [UDGT+26]

    Manu Upadhyaya, Shuvomoy Das Gupta, Adrien B. Taylor, Sebastian Banert, and Pontus Giselsson. The AutoLyap software suite for computer-assisted Lyapunov analyses of first-order methods. 2026. arXiv:2506.24076v2.

22. [UV22]

    André Uschmajew and Bart Vandereycken. A note on the optimal convergence rate of descent methods with fixed step sizes for smooth strongly convex functions. Journal of Optimization Theory and Applications, 194:364–373, 2022. doi:10.1007/s10957-022-02032-z.

23. [VSFL18]

    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, 2018. doi:10.1109/LCSYS.2017.2722406.

24. [Mos26]

    MOSEK ApS. MOSEK Fusion API for Python. Version 11.2. 2026. URL.

25. [dKGT20]

    Etienne de Klerk, François Glineur, and Adrien B. Taylor. Worst-case convergence analysis of inexact gradient and Newton methods through semidefinite programming performance estimation. SIAM Journal on Optimization, 30(3):2053–2082, 2020. doi:10.1137/19M1281368.