Using SymbolicAnalysis.jl for convexity certificates

In this tutorial, we will show how to use automatic convexity certification of the optimization problem using SymbolicAnalysis.jl.

This works with the structural_analysis keyword argument to OptimizationProblem. This tells the package to try to trace through the objective and constraints with symbolic variables (for more details on this look at the Symbolics documentation). This relies on the Disciplined Programming approach hence neccessitates the use of "atoms" from the SymbolicAnalysis.jl package.

We'll use a simple example to illustrate the convexity structure certification process.

using SymbolicAnalysis, Zygote, LinearAlgebra, Optimization

function f(x, p = nothing)
    return exp(x[1]) + x[1]^2
end

optf = OptimizationFunction(f, Optimization.AutoForwardDiff())
prob = OptimizationProblem(optf, [0.4], structural_analysis = true)

sol = solve(prob, Optimization.LBFGS(), maxiters = 1000)

retcode: Success
u: 1-element Vector{Float64}:
 -0.3517365907643902

The result can be accessed as the analysis_results field of the solution.

sol.cache.analysis_results.objective

SymbolicAnalysis.AnalysisResult(SymbolicAnalysis.Convex, SymbolicAnalysis.Positive, nothing)

Relatedly you can enable structural analysis in Riemannian optimization problems (supported only on the SPD manifold).

We'll look at the Riemannian center of mass of SPD matrices which is known to be a Geodesically Convex problem on the SPD manifold.

using Optimization, OptimizationManopt, Symbolics, Manifolds, Random, LinearAlgebra,
      SymbolicAnalysis

M = SymmetricPositiveDefinite(5)
m = 100
σ = 0.005
q = Matrix{Float64}(LinearAlgebra.I(5)) .+ 2.0

data2 = [exp(M, q, σ * rand(M; vector_at = q)) for i in 1:m];

f(x, p = nothing) = sum(SymbolicAnalysis.distance(M, data2[i], x)^2 for i in 1:5)
optf = OptimizationFunction(f, Optimization.AutoZygote())
prob = OptimizationProblem(optf, data2[1]; manifold = M, structural_analysis = true)

opt = OptimizationManopt.GradientDescentOptimizer()
sol = solve(prob, opt, maxiters = 100)

retcode: Failure
u: 5×5 Matrix{Float64}:
 2.99572  1.99688  1.99689  1.99583  1.99731
 1.99688  2.99767  1.99756  1.9971   1.99827
 1.99689  1.99756  2.99724  1.99676  1.99835
 1.99583  1.9971   1.99676  2.99603  1.99761
 1.99731  1.99827  1.99835  1.99761  2.9986