Add WeightedSOSCone (#351)

* Add WeightedSOSCone

* Fix

* Fix format

* Fix

* fix

* Add docstrings in docs

* Fix

* Fix

* rename

* Add tests

* Fix format

* Update ci.yml

* include tests

* MOI master

Author: blegat

.github/workflows/ci.yml

@@ -32,6 +32,14 @@ jobs:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
       - uses: julia-actions/cache@v1
+      - name: MOI
+        shell: julia --project=@. {0}
+        run: |
+          using Pkg
+          Pkg.add([
+              PackageSpec(name="MultivariateBases", rev="master"),
+              PackageSpec(name="MathOptInterface", rev="master"),
+          ])
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
         with:

docs/src/reference/standard_form.md

@@ -43,6 +43,12 @@ SumOfSquares.PolyJuMP.SAGE.AGEBridge
 
 ## MOI Sets
 
+Sum-of-Squares cones:
+```@docs
+SumOfSquares.SOSPolynomialSet
+SumOfSquares.WeightedSOSCone
+```
+
 Special cases of positive semidefinite cones:
 ```@docs
 SumOfSquares.EmptyCone

docs/src/references.bib

@@ -61,6 +61,19 @@ @Article{Hesse1973
   publisher = {Institute for Operations Research and the Management Sciences (INFORMS)},
 }
 
+@article{Kapelevich2023,
+  title = {Sum of Squares Generalizations for Conic Sets},
+  author = {Kapelevich, Lea and Coey, Chris and Vielma, Juan Pablo},
+  year = {2023},
+  month = may,
+  journal = {Mathematical Programming},
+  volume = {199},
+  number = {1},
+  pages = {1417--1429},
+  issn = {1436-4646},
+  doi = {10.1007/s10107-022-01831-6},
+}
+
 @Book{Lasserre2009,
   author    = {Lasserre, Jean Bernard},
   publisher = {Imperial College Press},
@@ -98,6 +111,18 @@ @Article{Murty1987
   publisher = {Springer-Verlag Berlin, Heidelberg}
 }
 
+@Article{Papp2017,
+  title = {Sum-of-Squares Optimization without Semidefinite Programming},
+  author = {Papp, D. and Y{\i}ld{\i}z, S.},
+  year = {2017},
+  month = dec,
+  journal = {ArXiv e-prints},
+  eprint = {1712.01792},
+  primaryclass = {math.OC},
+  adsnote = {Provided by the SAO/NASA Astrophysics Data System},
+  archiveprefix = {arxiv},
+}
+
 @Article{Parrilo2003,
   title = {Semidefinite Programming Relaxations for Semialgebraic Problems},
   author = {Parrilo, Pablo A},

src/Bridges/Constraint/utilities.jl

@@ -14,7 +14,7 @@ function union_set_types(MCT)
     }
 end
 function constrained_variable_types(MCT)
-    return [
+    return Tuple{Type}[
         (typeof(SOS.matrix_cone(MCT, 0)),),
         (typeof(SOS.matrix_cone(MCT, 1)),),
         (typeof(SOS.matrix_cone(MCT, 2)),),

src/Bridges/Variable/Variable.jl

@@ -1,11 +1,13 @@
 module Variable
 
+import MutableArithmetics as MA
 import MathOptInterface as MOI
-
+import MultivariatePolynomials as MP
 import SumOfSquares as SOS
 
 include("psd2x2.jl")
 include("scaled_diagonally_dominant.jl")
 include("copositive_inner.jl")
+include("kernel.jl")
 
 end

src/Bridges/Variable/kernel.jl

@@ -0,0 +1,125 @@
+struct KernelBridge{T,M} <: MOI.Bridges.Variable.AbstractBridge
+    affine::Vector{MOI.ScalarAffineFunction{T}}
+    variables::Vector{Vector{MOI.VariableIndex}}
+    constraints::Vector{MOI.ConstraintIndex{MOI.VectorOfVariables}}
+    set::SOS.WeightedSOSCone{M}
+end
+
+function MOI.Bridges.Variable.bridge_constrained_variable(
+    ::Type{KernelBridge{T}},
+    model::MOI.ModelLike,
+    set::SOS.WeightedSOSCone{M},
+) where {T,M}
+    variables = Vector{MOI.VariableIndex}[]
+    constraints = MOI.ConstraintIndex{MOI.VectorOfVariables}[]
+    acc = MA.Zero()
+    for (gram_basis, weight) in zip(set.gram_bases, set.weights)
+        gram, vars, con = SOS.add_gram_matrix(model, M, gram_basis, T)
+        push!(variables, vars)
+        push!(constraints, con)
+        acc = MA.add_mul!!(acc, weight, gram)
+    end
+    affine = MP.coefficients(acc, set.basis)
+    return KernelBridge{T,M}(affine, variables, constraints, set)
+end
+
+function MOI.Bridges.Variable.supports_constrained_variable(
+    ::Type{<:KernelBridge},
+    ::Type{<:SOS.WeightedSOSCone},
+)
+    return true
+end
+
+function MOI.Bridges.added_constrained_variable_types(
+    ::Type{KernelBridge{T,M}},
+) where {T,M}
+    return SOS.Bridges.Constraint.constrained_variable_types(M)
+end
+
+function MOI.Bridges.added_constraint_types(::Type{<:KernelBridge})
+    return Tuple{Type,Type}[]
+end
+
+# Attributes, Bridge acting as a model
+function MOI.get(bridge::KernelBridge, ::MOI.NumberOfVariables)
+    return sum(length, bridge.variables)
+end
+
+function MOI.get(bridge::KernelBridge, ::MOI.ListOfVariableIndices)
+    return reduce(vcat, bridge.variables)
+end
+
+function MOI.get(
+    bridge::KernelBridge,
+    ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S},
+) where {S<:MOI.AbstractVectorSet}
+    return count(bridge.constraints) do ci
+        return ci isa MOI.ConstraintIndex{MOI.VectorOfVariables,S}
+    end
+end
+
+function MOI.get(
+    bridge::KernelBridge,
+    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S},
+) where {S}
+    return [
+        ci for ci in bridge.constraints if
+        ci isa MOI.ConstraintIndex{MOI.VectorOfVariables,S}
+    ]
+end
+
+# Indices
+function MOI.delete(model::MOI.ModelLike, bridge::KernelBridge)
+    for vars in bridge.variables
+        MOI.delete(model, vars)
+    end
+    return
+end
+
+# Attributes, Bridge acting as a constraint
+
+function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::KernelBridge)
+    return bridge.set
+end
+
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintPrimal,
+    bridge::KernelBridge,
+)
+    return [
+        MOI.get(
+            model,
+            MOI.VariablePrimal(attr.result_index),
+            bridge,
+            MOI.Bridges.IndexInVector(i),
+        ) for i in eachindex(bridge.affine)
+    ]
+end
+
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.VariablePrimal,
+    bridge::KernelBridge,
+    i::MOI.Bridges.IndexInVector,
+)
+    return MOI.Utilities.eval_variable(bridge.affine[i.value]) do
+        return vi -> MOI.get(model, MOI.VariablePrimal(attr.result_index), vi)
+    end
+end
+
+function MOI.Bridges.bridged_function(
+    bridge::KernelBridge,
+    i::MOI.Bridges.IndexInVector,
+)
+    return bridge.affine[i.value]
+end
+
+function MOI.Bridges.Variable.unbridged_map(
+    bridge::KernelBridge{T},
+    coefs::Vector{MOI.VariableIndex},
+) where {T}
+    F = MOI.ScalarAffineFunction{T}
+    map = Pair{MOI.VariableIndex,F}[]
+    return nothing
+end

src/sets.jl

@@ -108,6 +108,66 @@ function Base.copy(
     return set
 end
 
+"""
+    struct WeightedSOSCone{
+        M,
+        B<:AbstractPolynomialBasis,
+        G<:AbstractPolynomialBasis,
+        W<:MP.AbstractPolynomialLike,
+    } <: MOI.AbstractVectorSet
+        basis::B
+        gram_bases::Vector{G}
+        weights::Vector{W}
+    end
+
+The weighted sum-of-squares cone is the set of vectors of coefficients `a` in `basis`
+that are the sum of `weights[i]` multiplied by a gram matrix with basis `gram_bases[i]`.
+The matrix cone type `M` is used to decide in which cone the gram matrix is constrained,
+e.g., `MOI.PositiveSemidefiniteConeTriangle` or
+[`SumOfSquares.ScaledDiagonallyDominantConeTriangle`](@ref).
+
+See [Papp2017; Section 1.1](@cite) and [Kapelevich2023; Section 1](@cite).
+"""
+struct WeightedSOSCone{
+    M,
+    B<:AbstractPolynomialBasis,
+    G<:AbstractPolynomialBasis,
+    W<:MP.AbstractPolynomialLike,
+} <: MOI.AbstractVectorSet
+    basis::B
+    gram_bases::Vector{G}
+    weights::Vector{W}
+end
+function WeightedSOSCone{M}(
+    basis::AbstractPolynomialBasis,
+    gram_bases::Vector{G},
+    weights::Vector{W},
+) where {M,G<:AbstractPolynomialBasis,W<:MP.AbstractPolynomialLike}
+    return WeightedSOSCone{M,typeof(basis),G,W}(basis, gram_bases, weights)
+end
+MOI.dimension(set::WeightedSOSCone) = length(set.basis)
+Base.copy(set::WeightedSOSCone) = set
+function Base.:(==)(a::WeightedSOSCone, b::WeightedSOSCone)
+    return a.basis == b.basis &&
+           a.gram_bases == b.gram_bases &&
+           a.weights == b.weights
+end
+
+"""
+    struct SOSPolynomialSet{
+        DT<:AbstractSemialgebraicSet,
+        MT<:MP.AbstractMonomial,
+        MVT<:AbstractVector{MT},
+        CT<:Certificate.AbstractCertificate,
+    } <: MOI.AbstractVectorSet
+        domain::DT
+        monomials::MVT
+        certificate::CT
+    end
+
+The sum-of-squares cone is the set of vectors of coefficients `a` for monomials `monomials`
+for which the polynomial is a sum-of-squares `certificate` over `domain`.
+"""
 struct SOSPolynomialSet{
     DT<:AbstractSemialgebraicSet,
     MT<:MP.AbstractMonomial,

test/Bridges/Bridges.jl

@@ -0,0 +1,12 @@
+using Test
+
+@testset "$(dir)" for dir in ["Variable"]
+    @testset "$(file)" for file in readdir(joinpath(@__DIR__, dir))
+        if !endswith(file, ".jl")
+            continue
+        elseif dir == "." && file == "Bridges.jl"
+            continue
+        end
+        include(joinpath(@__DIR__, dir, file))
+    end
+end

test/Bridges/Variable/kernel.jl

@@ -0,0 +1,46 @@
+module TestVariableKernel
+
+using Test
+using DynamicPolynomials
+using SumOfSquares
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function test_runtests()
+    @polyvar x y
+    MOI.Bridges.runtests(
+        SumOfSquares.Bridges.Variable.KernelBridge,
+        model -> begin
+            MOI.add_constrained_variables(
+                model,
+                SumOfSquares.WeightedSOSCone{
+                    MOI.PositiveSemidefiniteConeTriangle,
+                }(
+                    MonomialBasis([x^4, x^3 * y, x^2 * y^2, y^4]),
+                    [MonomialBasis([x^2, y^2, x * y])],
+                    [1.0 * x^0 * y^0],
+                ),
+            )
+        end,
+        model -> begin
+            Q, _ = MOI.add_constrained_variables(
+                model,
+                MOI.PositiveSemidefiniteConeTriangle(3),
+            )
+        end,
+    )
+    return
+end
+
+end  # module
+
+TestVariableKernel.runtests()

test/runtests.jl

@@ -32,6 +32,7 @@ include("gram_matrix.jl")
 
 include("variable.jl")
 include("constraint.jl")
+include("Bridges/Bridges.jl")
 
 include("Mock/mock_tests.jl")