Overview

Algorithms are mathematical descriptions of computation. AlgorithmAnalysis.jl implements a domain-specific language (DSL) that enables users to represent algorithms symbolically with syntax that closely resembles their mathematical description, and then manipulate the algorithm both symbolically and numerically. We now provide an overview of this process; for more details please see the API.

Symbolics

Each mathematical object is a symbolic expression in some space. AlgorithmAnalysis defines several common spaces, such as R for the field of real numbers, Rⁿ for a real finite-dimensional vector space, and Prop for the set of propositions. Following standard mathematical notation, we can create variables in a space using the @alg macro as follows:

julia> @alg a ∈ R
julia> @alg u ∈ Rⁿ
julia> @alg p, q ∈ Prop

Common methods are implemented for each of these spaces. For instance,

julia> p ∧ q  # conjunction p and q, the symbol `∧` can be typed by `\wedge<tab>`p ∧ q
julia> a * u  # scale vector u by aa * u

To assign these expressions to other variables, we again use the @alg macro:

julia> @alg r = p ∧ q

We can also represent functions:

julia> @alg f ∈ functional(Rⁿ)

This constructs a functional Rⁿ → R, which can then be evaluated at vectors in Rⁿ to produce scalars in R:

julia> @alg a = f(u)

Depending on the type of expression, we can form other expressions using basic algebraic operations. With the real numbers R, for instance, we can construct other scalars using addition and multiplication:

julia> @alg a, b ∈ R
julia> 2a + b2 * a + b

Here, the literal constant 2 is promoted to a symbolic expression in R, which can be done explicitly with R(2), and assigned the constant value 2. The vector space Rⁿ has similar algebraic operations as well as an inner product and norm:

julia> @alg u, v ∈ Rⁿ
julia> u'(v)   # inner product of u and vu'(v)
julia> u^2     # squared norm of u‖u‖²

The inner product and norm produce scalars in R. The inner product uses the function evaluation notation u'(v) since the adjoint u' of a vector is a linear functional Rⁿ → R which can then be evaluated at vectors to produce scalars. The notation u^2 is not standard mathematical notation, but is a shorthand for ‖u‖².

AlgorithmAnalysis can represent constraints on expressions as propositions. For instance, the proposition that two expressions are equal is formed as:

julia> a == 0a = 0

We can also create inequality constraints on real scalars:

julia> @alg a ∈ R
julia> a ≥ 00 ≤ a

For symmetric matrices, we can specify that a matrix is semidefinite as:

julia> @alg A ∈ Sⁿ
julia> A ⪰ 0A ⪰ 0

Numerics

While symbolic manipulations enable us to specify algorithms, analyzing them often requires numeric computation. To mix symbolic and numeric computations, AlgorithmAnalysis uses scoped values to construct a local scope in which symbolic expressions are instantiated with numeric values. As a simple example, we can substitute a numeric value for a parameter:

@alg let
    a ∈ R
    with_parameters(Dict(a => 2)) do
        evaluate(a)  # 2
    end
    evaluate(a)  # a
end

This code constructs a scalar a, defines a local scope in which it has the value 2, and then evaluates it within this scope. Outside of this scope, however, a is still a symbolic variable.

We can also formulate symbolic optimization problems and then solve them numerically using JuMP. For instance, the following feasibility problem is (trivially) true:

@alg let
    x ∈ R
    feas = feasible((x ≥ 1) ∧ (x ≤ 2))
    with_numerics() do
        evaluate(feas)  # true
    end
end

This constructs a feasibility problem using feasible, defines a local scope in which symbolic expressions are evaluated using a JuMP model using with_numerics, and then solves the feasibility problem within this scope to obtain the boolean literal true. While these examples are quite simple, we can use the same ideas to solve more general linear programs:

@alg let
    x, y ∈ R
    c1 = 50x + 24y ≤ 2400
    c2 = 30x + 33y ≤ 2100
    c3 = x ≥ 45
    c4 = y ≥ 5
    cons = c1 ∧ c2 ∧ c3 ∧ c4
    obj = x + y - 50
    opt = maximize(obj, cons)
    with_numerics() do
        evaluate(opt) ≈ 1.25 && evaluate(x) ≈ 45.0 && evaluate(y) ≈ 6.25  # true
    end
end

Since with_numerics creates a local scope, expressions that are obtained numerically maintain their numeric values within this scope. The above code first solves the optimization problem to obtain the numeric values for x and y, and then compares these numeric values against the known solution to the problem. Beyond linear programs, we can also solve semidefinite programs in a similar way:

@alg let
    x ∈ R
    A = [2 x; x 2]
    opt = maximize(x, A ⪰ 0)
    with_numerics() do
        evaluate(opt) ≈ 2.0  # true
    end
end

Transformations

At this point, AlgorithmAnalysis can be viewed as a domain-specific language for optimization, similar to other DSLs such as JuMP and Convex. The main benefits come when combining symbolic and numeric computations together, which can create emergent behavior that can be useful in analyzing algorithms. To that end, consider the following optimization problem:

julia> @alg begin
           α, L ∈ R, x, xs ∈ Rⁿ, f ∈ differentiable_functional(Rⁿ)
       
           gs   = f'(xs)
           g    = f'(x)
           init = (x - xs)^2
           x⁺   = x - α * g
           f⁺   = f(x⁺)
           c1   = smooth_convex(f, L)
           c2   = gs^2 == zero(R)
           c3   = init ≤ one(R)
           con  = c1 ∧ c2 ∧ c3
           obj  = f⁺ - f(xs)
           opt  = maximize(obj, con)
       end

which prints as:

julia> optmaximize     f(x - α * ∇f(x)) - f(xs)
subject to   f ∈ SmoothConvex(L)
             ‖∇f(xs)‖² = 0
             ‖x - xs‖² ≤ 1

We can view all leaf nodes in this optimization problem using leaves:

julia> leaves(opt)Set{SymbolicUtils.BasicSymbolic} with 5 elements:
  f
  L
  x
  xs
  α

One of the leaf variables in this problem is the differentiable function f : Rⁿ → R, so the problem is intractable as written. However, we can symbolically simplify this optimization problem by replacing the $L$-smooth convex function $f$ with its interpolation conditions between all of the points at which it is evaluated:

julia> transformed_opt = smooth_convex_interpolation(opt)[ Info: Applying L-smooth convex interpolation to function f with evaluation points xs, x⁺, x
maximize     f(x⁺) - f(xs)
subject to   f(xs) + ∇f(xs)'(xs - xs) + 1 / 2 * L * ‖∇f(xs) - ∇f(xs)‖² ≤ f(xs)
             f(x⁺) + ∇f(x⁺)'(xs - (x - α * ∇f(x))) + 1 / 2 * L * ‖∇f(xs) - ∇f(x⁺)‖² ≤ f(xs)
             f(x) + ∇f(x)'(xs - x) + 1 / 2 * L * ‖∇f(xs) - ∇f(x)‖² ≤ f(xs)
             f(xs) + ∇f(xs)'((x - α * ∇f(x)) - xs) + 1 / 2 * L * ‖∇f(x⁺) - ∇f(xs)‖² ≤ f(x⁺)
             f(x⁺) + ∇f(x⁺)'((x - α * ∇f(x)) - (x - α * ∇f(x))) + 1 / 2 * L * ‖∇f(x⁺) - ∇f(x⁺)‖² ≤ f(x⁺)
             f(x) + ∇f(x)'((x - α * ∇f(x)) - x) + 1 / 2 * L * ‖∇f(x⁺) - ∇f(x)‖² ≤ f(x⁺)
             f(xs) + ∇f(xs)'(x - xs) + 1 / 2 * L * ‖∇f(x) - ∇f(xs)‖² ≤ f(x)
             f(x⁺) + ∇f(x⁺)'(x - (x - α * ∇f(x))) + 1 / 2 * L * ‖∇f(x) - ∇f(x⁺)‖² ≤ f(x)
             f(x) + ∇f(x)'(x - x) + 1 / 2 * L * ‖∇f(x) - ∇f(x)‖² ≤ f(x)
             ‖∇f(xs)‖² = 0
             ‖x - xs‖² ≤ 1

This transformation has completely removed the function f from the problem. While some of the expressions appear to still depend on f, these have all been flattened to leaf expressions:

julia> leaves(transformed_opt)Set{SymbolicUtils.BasicSymbolic} with 10 elements:
  ∇f(x)
  f(x⁺)
  L
  f(xs)
  ∇f(x⁺)
  x
  f(x)
  xs
  ∇f(xs)
  α

For instance, the symbol ∇f(x) represents a leaf node in Rⁿ, and not the evaluation of the gradient of f at x. The optimization problem is still intractable though, as it involves vectors in Rⁿ whose dimension is arbitrarily large and therefore cannot be instantiated with numeric values. However, these vectors only appear in the optimization problem as inner products. A necessary and sufficient condition on these scalar values to be the inner product of such vectors is that the Gram matrix of all combinations of their inner products is positive semidefinite. We can therefore replace these vectors with their (flattened) inner products subject to the constraint that the Gram matrix is positive semidefinite. This is done automatically using:

julia> final_opt = gram_transformation(transformed_opt)[ Info: Applying Gram transformation to vector space Rⁿ with vectors ∇f(x⁺), xs, ∇f(x), x, ∇f(xs)
maximize     f(x⁺) - f(xs)
subject to   f(xs) + ⟨xs,∇f(xs)⟩ - ⟨xs,∇f(xs)⟩ + 1 / 2 * L * ((‖∇f(xs)‖² - ‖∇f(xs)‖²) - (‖∇f(xs)‖² - ‖∇f(xs)‖²)) ≤ f(xs)
             f(x⁺) + ⟨xs,∇f(x⁺)⟩ - (⟨x,∇f(x⁺)⟩ - α * ⟨∇f(x),∇f(x⁺)⟩) + 1 / 2 * L * ((‖∇f(xs)‖² - ⟨∇f(xs),∇f(x⁺)⟩) - (⟨∇f(xs),∇f(x⁺)⟩ - ‖∇f(x⁺)‖²)) ≤ f(xs)
             f(x) + ⟨xs,∇f(x)⟩ - ⟨x,∇f(x)⟩ + 1 / 2 * L * ((‖∇f(xs)‖² - ⟨∇f(x),∇f(xs)⟩) - (⟨∇f(x),∇f(xs)⟩ - ‖∇f(x)‖²)) ≤ f(xs)
             f(xs) + (⟨x,∇f(xs)⟩ - α * ⟨∇f(x),∇f(xs)⟩) - ⟨xs,∇f(xs)⟩ + 1 / 2 * L * ((‖∇f(x⁺)‖² - ⟨∇f(xs),∇f(x⁺)⟩) - (⟨∇f(xs),∇f(x⁺)⟩ - ‖∇f(xs)‖²)) ≤ f(x⁺)
             f(x⁺) + (⟨x,∇f(x⁺)⟩ - α * ⟨∇f(x),∇f(x⁺)⟩) - (⟨x,∇f(x⁺)⟩ - α * ⟨∇f(x),∇f(x⁺)⟩) + 1 / 2 * L * ((‖∇f(x⁺)‖² - ‖∇f(x⁺)‖²) - (‖∇f(x⁺)‖² - ‖∇f(x⁺)‖²)) ≤ f(x⁺)
             f(x) + (⟨x,∇f(x)⟩ - α * ‖∇f(x)‖²) - ⟨x,∇f(x)⟩ + 1 / 2 * L * ((‖∇f(x⁺)‖² - ⟨∇f(x),∇f(x⁺)⟩) - (⟨∇f(x),∇f(x⁺)⟩ - ‖∇f(x)‖²)) ≤ f(x⁺)
             f(xs) + ⟨x,∇f(xs)⟩ - ⟨xs,∇f(xs)⟩ + 1 / 2 * L * ((‖∇f(x)‖² - ⟨∇f(x),∇f(xs)⟩) - (⟨∇f(x),∇f(xs)⟩ - ‖∇f(xs)‖²)) ≤ f(x)
             f(x⁺) + ⟨x,∇f(x⁺)⟩ - (⟨x,∇f(x⁺)⟩ - α * ⟨∇f(x),∇f(x⁺)⟩) + 1 / 2 * L * ((‖∇f(x)‖² - ⟨∇f(x),∇f(x⁺)⟩) - (⟨∇f(x),∇f(x⁺)⟩ - ‖∇f(x⁺)‖²)) ≤ f(x)
             f(x) + ⟨x,∇f(x)⟩ - ⟨x,∇f(x)⟩ + 1 / 2 * L * ((‖∇f(x)‖² - ‖∇f(x)‖²) - (‖∇f(x)‖² - ‖∇f(x)‖²)) ≤ f(x)
             ‖∇f(xs)‖² = 0
             (‖x‖² - ⟨x,xs⟩) - (⟨x,xs⟩ - ‖xs‖²) ≤ 1
             [‖∇f(x⁺)‖² ⟨xs,∇f(x⁺)⟩ ⟨∇f(x),∇f(x⁺)⟩ ⟨x,∇f(x⁺)⟩ ⟨∇f(xs),∇f(x⁺)⟩; ⟨xs,∇f(x⁺)⟩ ‖xs‖² ⟨xs,∇f(x)⟩ ⟨x,xs⟩ ⟨xs,∇f(xs)⟩; ⟨∇f(x),∇f(x⁺)⟩ ⟨xs,∇f(x)⟩ ‖∇f(x)‖² ⟨x,∇f(x)⟩ ⟨∇f(x),∇f(xs)⟩; ⟨x,∇f(x⁺)⟩ ⟨x,xs⟩ ⟨x,∇f(x)⟩ ‖x‖² ⟨x,∇f(xs)⟩; ⟨∇f(xs),∇f(x⁺)⟩ ⟨xs,∇f(xs)⟩ ⟨∇f(x),∇f(xs)⟩ ⟨x,∇f(xs)⟩ ‖∇f(xs)‖²] ⪰ 0

The first 9 constraints are from the interpolation conditions, the next two are from the original problem, and the last constraint is that the Gram matrix of inner products is positive semidefinite. Again, this appears to still depend on vectors, but notice that all the inner products are now leaf expressions:

julia> leaves(final_opt)Set{SymbolicUtils.BasicSymbolic} with 20 elements:
  ⟨∇f(x),∇f(xs)⟩
  ⟨∇f(x),∇f(x⁺)⟩
  ⟨xs,∇f(x⁺)⟩
  ⟨x,∇f(xs)⟩
  ⟨x,xs⟩
  ‖∇f(x⁺)‖²
  α
  f(x⁺)
  ‖x‖²
  ‖xs‖²
  f(x)
  ‖∇f(x)‖²
  ⟨xs,∇f(x)⟩
  ⟨xs,∇f(xs)⟩
  ⟨x,∇f(x)⟩
  ‖∇f(xs)‖²
  ⟨∇f(xs),∇f(x⁺)⟩
  ⟨x,∇f(x⁺)⟩
  L
  f(xs)

All of these quantities are symbolic variables in R, and not function evaluations. To help distinguish between inner products of vectors and their flattened scalar variables, the inner product between u and v is printed as u'(v) while the corresponding flattened scalar leaf prints as ⟨u,v⟩ when the vectors are different and ‖u‖² when they are the same. The optimization problem has now been transformed to a semidefinite program (once the parameters α and L are fixed) and can therefore be solved numerically:

julia> with_numerics(parameters = Dict(α => 0.075, L => 10.0)) do
           evaluate(final_opt)
       end2.000040600426323

While we manually selected which transformations to apply and in what order (smooth convex interpolation followed by the Gram transformation), this process can be automated. These transformations and more are implemented in simplify, which automatically reconstructs the final semidefinite program from the original optimization problem as:

julia> final_opt = simplify(opt);[ Info: Applying L-smooth convex interpolation to function f with evaluation points xs, x⁺, x
[ Info: Applying Gram transformation to vector space Rⁿ with vectors ∇f(x⁺), xs, ∇f(x), x, ∇f(xs)

Also, numerically evaluating a symbolic expression does not modify the expression (which is immutable). Therefore, an expression can be evaluated again with different parameters:

julia> with_numerics(parameters = Dict(α => 0.075, L => 10.0)) do
           evaluate(final_opt)
       end2.000040600426323
julia> with_numerics(parameters = Dict(α => 1/3, L => 3.0)) do
           evaluate(final_opt)
       end0.5000111673778475

For more details on this approach to algorithm analysis, see the performance estimation section of the manual.

Iterative algorithms as dynamical systems

We can interpret iterative algorithms as dynamical systems, which can then be analyzed using tools from control theory. Dynamical systems have a state that evolves over time. At each point in time, the state completely describes the behavior of the system. While time may either be continuous or discrete, algorithms are naturally represented as discrete-time systems in which the discrete time index represents the iteration of the algorithm. To describe the algorithm, we must therefore describe how the state changes between each iteration. When the update has the same form at each iteration, the system is time invariant and we only need to specify the update at some (generic) iteration.

We represent a state transition as:

julia> @alg x₁, x₂ ∈ R
julia> @alg t = x₁ → x₂

This defines a state transition t (which is a proposition) that specifies that x₁ is a component of the algorithm state whose value at the next iteration is x₂. An algorithm may have more than one component in the state, and any node on the left side of a transition is considered part of the state. For instance, a state with multiple components is

julia> @alg y₁, y₂ ∈ Rⁿ
julia> @alg t = (x₁ → x₂) ∧ (y₁ → y₂)

This specifies that the state consists of both x₁ and y₁ whose values at the next iteration are x₂ and y₂, respectively.

As a more interesting example, consider the following algorithm, which represents gradient descent applied to a differentiable function:

julia> @alg begin
           α, μ, L ∈ R
           x, xs ∈ Rⁿ
           f ∈ differentiable_functional(Rⁿ)
           gs = f'(xs)
           g  = f'(x)
           x₊ = x - α * g
           t1 = x → x₊
           t2 = xs → xs
           t3 = f → f
           c1 = sector_bounded(f, μ, L)
           c2 = gs^2 == zero(R)
           con = t1 ∧ t2 ∧ t3 ∧ c1 ∧ c2
           perf = (x - xs)^2
       end

The main algorithmic update is the first transition, t1 = x → x₊ where x represents the iterate of the algorithm at some generic point in time and x₊ = x - α * f'(x) its value at the next iteration. The other two transitions are for the stationary point xs and the objective function f, which are specified to be time invariant by the transitions t2 = xs → xs and t3 = f → f.

Given these symbolic expressions, we can symbolically construct the problem of searching for a certificate that guarantees convergence of the algorithm with a particular rate:

julia> @alg ρ ∈ R
julia> prob = certify(ρ, perf, con)Lyapunov-based stability certification
  Rate:        ρ
  Performance: ‖x - xs‖²
  Constraints: x → (x - α * ∇f(x))
               xs → xs
               f → f
               f ∈ SectorBounded(μ, L)
               ‖∇f(xs)‖² = 0

This high-level symbolic object represents searching for a Lyapunov-based stability certificate for the given rate, performance measure, and constraint. Before formulating the search for a Lyapunov function that certifies stability, however, we must first apply several transformations to make the problem tractable. As before, the problem contains the function f as a variable, which cannot be used as a numeric decision variable. Instead, we replace the function with the sector bound:

julia> transformed_prob1 = sector_bounded_interpolation(prob)[ Info: Applying [μ,L] sector bounded interpolation to function f with evaluation points x, xs
Lyapunov-based stability certification
  Rate:        ρ
  Performance: ‖x - xs‖²
  Constraints: x → (x - α * ∇f(x))
               xs → xs
               ‖∇f(xs)‖² = 0
               (∇f(x) - μ * x)'(∇f(x) - L * x) ≤ 0
               (∇f(xs) - μ * xs)'(∇f(xs) - L * xs) ≤ 0

The function f has now been completely removed from the problem:

julia> leaves(transformed_prob1)Set{SymbolicUtils.BasicSymbolic} with 7 elements:
  ∇f(x)
  ρ
  L
  x
  xs
  ∇f(xs)
  α

Some of these leaves are vectors in Rⁿ, which also cannot be used as numeric decision variables. Replacing these vectors with the constraint that their Gram matrix is positive semidefinite yields the transformed problem:

julia> transformed_prob2 = gram_transformation(transformed_prob1)[ Info: Applying Gram transformation to vector space Rⁿ with vectors xs, ∇f(x), x, ∇f(xs)
Lyapunov-based stability certification
  Rate:        ρ
  Performance: (‖x‖² - ⟨x,xs⟩) - (⟨x,xs⟩ - ‖xs‖²)
  Constraints: ‖∇f(xs)‖² = 0
               (‖∇f(x)‖² - L * ⟨x,∇f(x)⟩) - (μ * ⟨x,∇f(x)⟩ - μ * L * ‖x‖²) ≤ 0
               (‖∇f(xs)‖² - L * ⟨xs,∇f(xs)⟩) - (μ * ⟨xs,∇f(xs)⟩ - μ * L * ‖xs‖²) ≤ 0
               ‖xs‖² → ‖xs‖²
               ⟨x,xs⟩ → (⟨x,xs⟩ - α * ⟨xs,∇f(x)⟩)
               ‖x‖² → ((‖x‖² - α * ⟨x,∇f(x)⟩) - (α * ⟨x,∇f(x)⟩ - α * α * ‖∇f(x)‖²))
               [‖xs‖² ⟨xs,∇f(x)⟩ ⟨x,xs⟩ ⟨xs,∇f(xs)⟩; ⟨xs,∇f(x)⟩ ‖∇f(x)‖² ⟨x,∇f(x)⟩ ⟨∇f(x),∇f(xs)⟩; ⟨x,xs⟩ ⟨x,∇f(x)⟩ ‖x‖² ⟨x,∇f(xs)⟩; ⟨xs,∇f(xs)⟩ ⟨∇f(x),∇f(xs)⟩ ⟨x,∇f(xs)⟩ ‖∇f(xs)‖²] ⪰ 0

All inner products have now been flattened to scalar leaf variables:

julia> leaves(transformed_prob2)Set{SymbolicUtils.BasicSymbolic} with 14 elements:
  ⟨∇f(x),∇f(xs)⟩
  ρ
  μ
  ⟨x,xs⟩
  ⟨x,∇f(xs)⟩
  α
  ‖x‖²
  ‖xs‖²
  ‖∇f(x)‖²
  ⟨xs,∇f(x)⟩
  ⟨x,∇f(x)⟩
  ⟨xs,∇f(xs)⟩
  ‖∇f(xs)‖²
  L

Also, recall that the state transitions in the original problem involved vectors like x and xs. Before flattening inner products, these transitions were first propagated to all applicable inner products, which resulted in the following three transitions:

‖xs‖² → ‖xs‖²
⟨x,xs⟩ → (⟨x,xs⟩ - α * ⟨xs,∇f(x)⟩)
‖x‖² → ((‖x‖² - α * ⟨x,∇f(x)⟩) - (α * ⟨x,∇f(x)⟩ - α * α * ‖∇f(x)‖²))

These transitions were able to be inferred from those of x and xs. Other inner products, such as ⟨x,∇f(xs)⟩, however, do not appear as states with a transition since at least one component of the inner product is not a state and therefore has no available transition (in this case, the gradient ∇f(xs)).

Now that all expressions in the problem are scalars in R, we can formulate the search for a Lyapunov function that certifies convergence with the given rate. Doing so yields the final semidefinite program:

julia> transformed_prob3 = with_parameters(Dict(ρ => 0.81, α => 0.1, μ => 1.0, L => 10.0)) do
           lyapunov_transformation(transformed_prob2)
       end[ Info: Applying Lyapunov transformation to performance measure performance with rate ρ
feasible
subject to   0 ≤ λ₂
             0 ≤ λ₃
             [λ₄ λ₅ λ₆ λ₇; λ₅ λ₈ λ₉ λ₁₀; λ₆ λ₉ λ₁₁ λ₁₂; λ₇ λ₁₀ λ₁₂ λ₁₃] ⪰ 0
             0 ≤ λ₁₅
             0 ≤ λ₁₆
             [λ₁₇ λ₁₈ λ₁₉ λ₂₀; λ₁₈ λ₂₁ λ₂₂ λ₂₃; λ₁₉ λ₂₂ λ₂₄ λ₂₅; λ₂₀ λ₂₃ λ₂₅ λ₂₆] ⪰ 0
             1 + (-1) * θ₃ + (-1) * μ * L * λ₂ + 1 * λ₁₁ = 0
             1 * λ₁₀ + 1 * λ₁₀ = 0
             1 + (-1) * θ₁ + (-1) * μ * L * λ₃ + 1 * λ₄ = 0
             1 * L * λ₂ + 1 * μ * λ₂ + 1 * λ₉ + 1 * λ₉ = 0
             1 * L * λ₃ + 1 * μ * λ₃ + 1 * λ₇ + 1 * λ₇ = 0
             (-1) * λ₂ + 1 * λ₈ = 0
             1 * λ₅ + 1 * λ₅ = 0
             1 * λ₁ + (-1) * λ₃ + 1 * λ₁₃ = 0
             -1 + -1 + (-1) * θ₂ + 1 * λ₆ + 1 * λ₆ = 0
             1 * λ₁₂ + 1 * λ₁₂ = 0
             1 * 0 * λ₂ + 1 * 0 * λ₃ = 0
             1 * θ₃ + (-1) * θ₃ * ρ + (-1) * μ * L * λ₁₅ + 1 * λ₂₄ = 0
             1 * λ₂₃ + 1 * λ₂₃ = 0
             1 * θ₁ + (-1) * θ₁ * ρ + (-1) * μ * L * λ₁₆ + 1 * λ₁₇ = 0
             (-1) * α * θ₃ + (-1) * α * θ₃ + 1 * L * λ₁₅ + 1 * μ * λ₁₅ + 1 * λ₂₂ + 1 * λ₂₂ = 0
             1 * L * λ₁₆ + 1 * μ * λ₁₆ + 1 * λ₂₀ + 1 * λ₂₀ = 0
             1 * α * α * θ₃ + (-1) * λ₁₅ + 1 * λ₂₁ = 0
             (-1) * α * θ₂ + 1 * λ₁₈ + 1 * λ₁₈ = 0
             1 * λ₁₄ + (-1) * λ₁₆ + 1 * λ₂₆ = 0
             1 * θ₂ + (-1) * θ₂ * ρ + 1 * λ₁₉ + 1 * λ₁₉ = 0
             1 * λ₂₅ + 1 * λ₂₅ = 0
             (-1) * 0 * ρ + 1 * 0 * λ₁₅ + 1 * 0 * λ₁₆ = 0

While this problem appears quite complex, it is a semidefinite program that is systematically constructed from the transformed problem.

Now that the problem of searching for a Lyapunov certificate of convergence has been formulated as a semidefinite program, we can evaluate it numerically as before:

julia> with_numerics(parameters = Dict(ρ => 0.81, α => 0.1, μ => 1.0, L => 10.0)) do
           evaluate(transformed_prob3)
       endtrue

For this problem, the symbolic simplification process required applying sector_bounded_interpolation followed by gram_transformation and then lyapunov_transformation. Again, this process is automated with simplify, which systematically constructs the final semidefinite program directly from the original problem:

julia> transformed_prob = with_parameters(Dict(ρ => 0.81, α => 0.1, μ => 1.0, L => 10.0)) do
           simplify(prob)
       end;[ Info: Applying [μ,L] sector bounded interpolation to function f with evaluation points x, xs
[ Info: Applying Gram transformation to vector space Rⁿ with vectors xs, ∇f(x), x, ∇f(xs)
[ Info: Applying Lyapunov transformation to performance measure performance with rate ρ

While certify constructed the search for a Lyapunov certificate for a particular rate, we often want to compute the fastest (smallest) rate that can be certified. This can be constructed as:

julia> opt = rate(perf, con)Lyapunov-based stability certification
  Rate:        ρₒₚₜ
  Performance: ‖x - xs‖²
  Constraints: x → (x - α * ∇f(x))
               xs → xs
               f → f
               f ∈ SectorBounded(μ, L)
               ‖∇f(xs)‖² = 0

Evaluating this node performs a bisection search to find the smallest rate for which certify is true:

julia> with_parameters(Dict(α => 0.1, μ => 1.0, L => 10.0)) do
           evaluate(simplify(opt))
       end[ Info: Applying [μ,L] sector bounded interpolation to function f with evaluation points x, xs
[ Info: Applying Gram transformation to vector space Rⁿ with vectors xs, ∇f(x), x, ∇f(xs)
[ Info: Applying Lyapunov transformation to performance measure performance with rate nothing
[ Info: Bisecting over ρ in [0.0, 1.0] with tolerance 1.0e-8
0.8100000023841858

For more details on this approach to algorithm analysis, see the Lyapunov analysis section of the manual.