Notebook

Univariate Solver¶

Contributed by: Benoît Legat

When using an SDP solver, the Sum-of-Squares constraint is bridged into a semidefinite constraint. This reformulation is done only if the solver does not support the Sum-of-Squares constraint. We show in this tutorial how to define a solver that supports such constraint. The same model with be then solved with or without reformulation depending on the solver.

We consider a solver that would support finding the SOS decomposition of nonnegative univariate polynomials.

In [1]:

module MyUnivariateSolver

import LinearAlgebra
import MathOptInterface as MOI
import MultivariatePolynomials as MP
import MultivariateBases as MB
import SumOfSquares as SOS

function decompose(p::MP.AbstractPolynomial, tol=1e-6)
    vars = MP.effective_variables(p)
    if length(vars) != 1
        error("$p is not univariate")
    end
    x = first(vars)
    lead = MP.leading_coefficient(p)
    if !isone(lead)
        p = p / lead
    end
    deg = MP.maxdegree(p)
    if isodd(deg)
        return
    end
    d = div(deg, 2)
    companion = zeros(2d, 2d)
    for i in 0:(2d-1)
        if i > 0
            companion[i + 1, i] = 1.0
        end
        companion[i + 1, end] = -MP.coefficient(p, x^i)
    end
    F = LinearAlgebra.schur(complex(companion))
    q = one(p)
    i = 1
    while i <= length(F.values)
        root = F.values[i]
        q *= (x - root)
        if !isapprox(real(root), real(F.values[i+1]), rtol=tol, atol=tol)
            return # Cannot happen for complex conjugate root so it means that we have a root which does not have an even multiplicity This means that the polynomial is not nonnegative
        end
        i += 2
    end
    q1 = MP.map_coefficients(real, q)
    q2 = MP.map_coefficients(imag, q)
    return SOS.SOSDecomposition([MB.algebra_element(q1), MB.algebra_element(q2)])
end

mutable struct Optimizer <: MOI.AbstractOptimizer
    p::Union{Nothing,MP.AbstractPolynomial}
    decomposition::Union{Nothing,SOS.SOSDecomposition}
    tol::Float64
    function Optimizer()
        return new(nothing, nothing, 1e-6)
    end
end

MOI.is_empty(optimizer::Optimizer) = optimizer.p === nothing
function MOI.empty!(optimizer::Optimizer)
    optimizer.p = nothing
    return
end

function MOI.supports_constraint(::Optimizer, ::Type{<:MOI.VectorAffineFunction}, ::Type{<:SOS.SOSPolynomialSet{SOS.FullSpace}})
    return true
end
function MOI.add_constraint(optimizer::Optimizer, func::MOI.VectorAffineFunction, set::SOS.SOSPolynomialSet{SOS.FullSpace})
    if optimizer.p !== nothing
        error("Only one constraint is supported")
    end
    if !isempty(func.terms)
        error("Only supports constant polynomials")
    end
    optimizer.p = MP.polynomial(MB.algebra_element(func.constants, set.basis))
    return MOI.ConstraintIndex{typeof(func),typeof(set)}(1) # There will be only ever one constraint so the index does not matter.
end

MOI.supports_incremental_interface(::Optimizer) = true
function MOI.copy_to(optimizer::Optimizer, model::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(optimizer, model)
end
function MOI.optimize!(optimizer::Optimizer)
    optimizer.decomposition = decompose(optimizer.p, optimizer.tol)
end

function MOI.get(optimizer::Optimizer, ::MOI.TerminationStatus)
    if optimizer.decomposition === nothing
        return MOI.INFEASIBLE
    else
        return MOI.OPTIMAL
    end
end

function MOI.get(optimizer::Optimizer, ::MOI.PrimalStatus)
    if optimizer.decomposition === nothing
        return MOI.NO_SOLUTION
    else
        return MOI.FEASIBLE_POINT
    end
end

function MOI.get(optimizer::Optimizer, ::SOS.SOSDecompositionAttribute, ::MOI.ConstraintIndex)
    return optimizer.decomposition
end

end

Out[1]:

Main.var"##17010".MyUnivariateSolver

We define the same function both for our new solver and SDP solvers.

In [2]:

using SumOfSquares
function decompose(p, solver)
    model = Model(solver)
    con = @constraint(model, p in SOSCone())
    optimize!(model)
    @assert primal_status(model) == MOI.FEASIBLE_POINT
    return sos_decomposition(con)
end

Out[2]:

decompose (generic function with 1 method)

We consider the following univariate polynomial from Blekherman2012; Example 3.35.

Using our solver we find the following decomposition in two squares.

In [3]:

using DynamicPolynomials
@polyvar x
p = x^4 + 4x^3 + 6x^2 + 4x + 5
dec = decompose(p, MyUnivariateSolver.Optimizer)

Out[3]:

(-0.9999999999999989·1 + 2.0000000000000018·x + 1.0·x^2)^2 + (-2.000000000000003·1 - 2.0·x)^2

We can verify as follows that it gives a correct decomposition:

In [4]:

polynomial(dec)

Out[4]:

$5.000000000000011 + 4.000000000000013x + 6.000000000000009x^{2} + 4.0000000000000036x^{3} + x^{4}$

We can also use a semidefinite solver:

In [5]:

import CSDP
dec = decompose(p, CSDP.Optimizer)

Iter: 17 Ap: 5.20e-01 Pobj:  0.0000000e+00 Ad: 6.52e-01 Dobj:  2.9486799e-05 
Success: SDP solved
Primal objective value: 0.0000000e+00 
Dual objective value: 0.0000000e+00 
Relative primal infeasibility: 9.76e-09 
Relative dual infeasibility: 1.89e-47 
Real Relative Gap: 0.00e+00 
XZ Relative Gap: 3.58e-44 
DIMACS error measures: 1.18e-08 0.00e+00 1.89e-47 0.00e+00 0.00e+00 3.58e-44
CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00 
Iter:  1 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.6555556e+01 
Iter:  2 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  4.6844813e+00

Out[5]:

(-1.4271293634654956·1 - 2.4208468386942927·x - 0.608679682768493·x^2)^2 + (1.7193818199786142·1 - 0.8390365176838647·x - 0.6942919122813989·x^2)^2 + (0.08383279228081496·1 - 0.14597477689220884·x + 0.3840153438672542·x^2)^2

The decomposition is different, it is the sum of 3 squares. However, it is also valid:

In [6]:

polynomial(dec)

Out[6]:

$5.000000000000001 + 3.9999999999999956x + 5.999999999999998x^{2} + 4.0x^{3} + 1.0000000000000004x^{4}$

This notebook was generated using Literate.jl.