Add tutorial on writing SOS solver

Author: blegat

docs/Project.toml

@@ -10,6 +10,7 @@ GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"

docs/make.jl

@@ -13,6 +13,7 @@ const _TUTORIAL_SUBDIR = [
     "Noncommutative and Hermitian",
     "Sparsity",
     "Symmetry",
+    "Extension",
 ]
 
 function literate_directory(dir)

docs/src/constraints.md

@@ -324,6 +324,7 @@ SumOfSquares.LagrangianMultipliers
 lagrangian_multipliers
 SOSDecomposition
 SOSDecompositionWithDomain
+SOSDecompositionAttribute
 sos_decomposition
 ```
 

docs/src/tutorials/Extension/tmp.jl

@@ -0,0 +1,61 @@
+function decompose(p::MP.AbstractPolynomial, tol)
+    vars = MP.effective_variables(p)
+    if length(vars) != 1
+        error("$p is not univariate")
+    end
+    x = first(vars)
+    lead = MP.leadingcoefficient(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
+    #display(companion)
+    F = LinearAlgebra.schur(complex(companion))
+    #display(F.Z * F.T * F.Z')
+    σ = vcat(1:2:(2d - 1), 2:2:2d)
+    S = LinearAlgebra.ordschur(F, isodd.(1:2d))
+    #U = S.Z[σ, σ]
+    #S.Z * S.T * S.Z'
+    #display(σ)
+    #display(S.T)
+    #S.T[σ, σ]
+    vt = S.Z[1, 1:d]
+    I = 1:d
+    J = d .+ I
+    U11 = S.Z[I, I]
+    U22 = S.Z[J, J]
+    U12 = S.Z[I, J]
+    display(S.Z)
+    return U12
+    U21 = S.Z[J, I]
+    vt = U22[:, 1]
+    U = [U11 U12; U21 U22]
+    display(S.T)
+    display(U * S.T * U')
+    display(S.Z * S.T * S.Z')
+    display(U * U')
+    @show lead
+    q = lead * U12' \ vt
+    q = vec(vt' * inv(U12))
+    @show q
+    q1 = lead * x^d + sum([-real(q[i + 1]) * x^i for i in 0:(d-1)])
+    q2 = sum([-imag(q[i + 1]) * x^i for i in 0:(d-1)])
+    q_1 = Polynomial([-real.(q); 1.0])
+    q_2 = Polynomial([-imag.(q); 0.0])
+    @show q1^2 + q2^2
+    @show q_1^2 + q_2^2
+    @show Polynomials.roots(q_1 + im * q_2)
+    @show Polynomials.roots(q_1 - im * q_2)
+    return q1, q2
+end

docs/src/tutorials/Extension/univariate_solver.jl

@@ -0,0 +1,164 @@
+# # Univariate Solver
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Extension/univariate_solver.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Extension/univariate_solver.ipynb)
+# **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.
+
+using Test #src
+
+module MyUnivariateSolver
+
+import LinearAlgebra
+import MathOptInterface
+const MOI = MathOptInterface
+import MultivariatePolynomials
+const MP = MultivariatePolynomials
+import SumOfSquares
+const SOS = SumOfSquares
+
+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.leadingcoefficient(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)
+            # 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
+            return
+        end
+        i += 2
+    end
+    q1 = MP.mapcoefficientsnz(real, q)
+    q2 = MP.mapcoefficientsnz(imag, q)
+    return SOS.SOSDecomposition([q1, 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(func.constants, set.monomials)
+    # There will be only ever one constraint so the index does not matter.
+    return MOI.ConstraintIndex{typeof(func),typeof(set)}(1)
+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
+
+# We define the same function both for our new solver and SDP solvers.
+
+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
+
+# We consider the following univariate polynomial from
+# Example 3.35 of [BPT12].
+#
+# [BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
+# *Semidefinite Optimization and Convex Algebraic Geometry*.
+# Society for Industrial and Applied Mathematics, **2012**.
+
+# Using our solver we find the following decomposition in two squares.
+
+using DynamicPolynomials
+@polyvar x
+p = x^4 + 4x^3 + 6x^2 + 4x + 5
+dec = decompose(p, MyUnivariateSolver.Optimizer)
+
+# We can verify as follows that it gives a correct decomposition:
+
+@test polynomial(dec) ≈ p #src
+polynomial(dec)
+
+# We can also use a semidefinite solver:
+
+import CSDP
+dec = decompose(p, CSDP.Optimizer)
+@test length(dec.ps) == 3 #src
+@test polynomial(dec) ≈ p #src
+
+# The decomposition is different, it is the sum of 3 squares.
+# However, it is also valid:
+
+polynomial(dec)

src/SumOfSquares.jl

@@ -17,6 +17,7 @@ Reexport.@reexport using SemialgebraicSets
 Reexport.@reexport using MultivariateMoments
 
 include("gram_matrix.jl")
+include("sosdec.jl")
 
 using PolyJuMP
 abstract type SOSLikeCone <: PolyJuMP.PolynomialSet end
@@ -53,7 +54,6 @@ include("Bridges/Bridges.jl")
 
 Reexport.@reexport using JuMP
 
-include("sosdec.jl")
 include("utilities.jl")
 include("variable.jl")
 include("constraint.jl")

src/attributes.jl

@@ -19,6 +19,45 @@ struct GramMatrixAttribute <: MOI.AbstractConstraintAttribute
 end
 GramMatrixAttribute() = GramMatrixAttribute(1)
 
+"""
+    struct SOSDecompositionAttribute
+        ranktol::Real
+        dec::MultivariateMoments.LowRankChol
+    end
+
+A constraint attribute for the [`SOSDecomposition`](@ref) of a constraint.
+By default, it is computed using
+`SOSDecomposition(gram, ranktol, dec)` where `gram` is the value of the
+[`GramMatrixAttribute`](@ref).
+"""
+struct SOSDecompositionAttribute <: MOI.AbstractConstraintAttribute
+    ranktol::Real
+    dec::MultivariateMoments.LowRankChol
+    result_index::Int
+end
+function SOSDecompositionAttribute(ranktol::Real, dec::MultivariateMoments.LowRankChol)
+    return SOSDecompositionAttribute(ranktol, dec, 1)
+end
+
+function MOI.get_fallback(
+    model::MOI.ModelLike,
+    attr::SOSDecompositionAttribute,
+    ci::MOI.ConstraintIndex,
+)
+    gram = MOI.get(model, attr.result_index, ci)
+    return SOSDecomposition(gram, attr.ranktol, attr.dec)
+end
+# TODO bridges should redirect to `MOI.get_fallback` as well so that
+# we can just use `Union{MOI.ConstraintIndex,MOI.Bridges.AbstractBridge}` above:
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::SOSDecompositionAttribute,
+    bridge::Union{Bridges.Constraint.SOSPolynomialBridge,Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge},
+)
+    gram = MOI.get(model, GramMatrixAttribute(attr.result_index), bridge)
+    return SOSDecomposition(gram, attr.ranktol, attr.dec)
+end
+
 """
     MomentMatrixAttribute(N)
     MomentMatrixAttribute()
@@ -49,6 +88,7 @@ LagrangianMultipliers() = LagrangianMultipliers(1)
 # optimize, even of `CertificateBasis` which is set befor optimize.
 function MOI.is_set_by_optimize(::Union{CertificateBasis,
                                         GramMatrixAttribute,
+                                        SOSDecompositionAttribute,
                                         MomentMatrixAttribute,
                                         LagrangianMultipliers})
     return true
@@ -64,7 +104,7 @@ function MOI.Bridges.Constraint.invariant_under_function_conversion(::Union{
 end
 
 # This is type piracy but we tolerate it.
-const ObjectWithoutIndex = AbstractGramMatrix{<:MOI.Utilities.ObjectWithoutIndex}
+const ObjectWithoutIndex = Union{AbstractGramMatrix{<:MOI.Utilities.ObjectWithoutIndex},SOSDecomposition{<:MOI.Utilities.ObjectWithoutIndex}}
 const ObjectOrTupleWithoutIndex = Union{ObjectWithoutIndex, Tuple{Vararg{ObjectWithoutIndex}}}
 const ObjectOrTupleOrArrayWithoutIndex = Union{ObjectOrTupleWithoutIndex, AbstractArray{<:ObjectOrTupleWithoutIndex}}
 MOI.Utilities.map_indices(::Function, x::ObjectOrTupleOrArrayWithoutIndex) = x

src/constraint.jl

@@ -226,6 +226,27 @@ function gram_matrix(cref::JuMP.ConstraintRef)
     return MOI.get(cref.model, GramMatrixAttribute(), cref)
 end
 
+"""
+    sos_decomposition(cref::JuMP.ConstraintRef)
+
+Return the [`SOSDecompositionAttribute`](@ref) of `cref`.
+"""
+function sos_decomposition(cref::JuMP.ConstraintRef, ranktol::Real=0.0,
+                           dec::MultivariateMoments.LowRankChol=SVDChol())
+    return MOI.get(cref.model, SOSDecompositionAttribute(ranktol, dec), cref)
+end
+
+"""
+    sos_decomposition(cref::JuMP.ConstraintRef, K<:AbstractBasicSemialgebraicSet)
+
+Return representation in the quadraic module associated with K.
+"""
+function sos_decomposition(cref::JuMP.ConstraintRef, K::AbstractBasicSemialgebraicSet, args...)
+    lm = SOSDecomposition.(lagrangian_multipliers(cref), args...)
+    gm = sos_decomposition(cref, args...)
+    return SOSDecompositionWithDomain(gm, lm, K)
+end
+
 """
     moment_matrix(cref::JuMP.ConstraintRef)
 

src/sosdec.jl

@@ -93,15 +93,6 @@ function MP.polynomial(decomp::SOSDecomposition, T::Type)
     return MP.polynomial(MP.polynomial(decomp), T)
 end
 
-"""
-    function sos_decomposition(cref::JuMP.ConstraintRef, ranktol::Float64, dec::MultivariateMoments.lowrankchol)
-
-Return representation as a sum of squares.
-"""
-function sos_decomposition(cref::JuMP.ConstraintRef, args...)
-    return SOSDecomposition(gram_matrix(cref), args...)
-end
-
 """
     struct SOSDecompositionWithDomain{T, PT, S}
 
@@ -145,14 +136,3 @@ end
 function Base.isapprox(p::SOSDecompositionWithDomain, q::SOSDecompositionWithDomain; kwargs...)
     return isapprox(p.sos, q.sos) && all(isapprox.(p.sosj, q.sosj)) && p.domain == q.domain
 end
-
-"""
-    sos_decomposition(cref::JuMP.ConstraintRef, K<:AbstractBasicSemialgebraicSet)
-
-Return representation in the quadraic module associated with K.
-"""
-function sos_decomposition(cref::JuMP.ConstraintRef, K::AbstractBasicSemialgebraicSet, args...)
-    lm = SOSDecomposition.(lagrangian_multipliers(cref), args...)
-    gm = sos_decomposition(cref, args...)
-    return SOSDecompositionWithDomain(gm, lm, K)
-end