Sampling

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+BMSOS = "288039e9-afd1-4944-b7ae-3ac2e056f6e9"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CSDP = "0a46da34-8e4b-519e-b418-48813639ff34"
 Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e"
@@ -12,6 +13,7 @@ Dualization = "191a621a-6537-11e9-281d-650236a99e60"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
+Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2"
 ImplicitPlots = "55ecb840-b828-11e9-1645-43f4a9f9ace7"
 Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
@@ -28,10 +30,12 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13"
+SDPLR = "56161740-ea4e-4253-9d15-43c62ff94d95"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
 SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"
 SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2"
+TrigPolys = "bbdedc48-cb31-4a37-9fe3-b015aecc8dd3"
 TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
 
 [compat]

docs/src/tutorials/Getting started/sampling.jl

@@ -0,0 +1,78 @@
+# # Sampling basis
+
+#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/sampling.ipynb)
+#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/sampling.ipynb)
+# **Contributed by**: Benoît Legat
+
+using Test #src
+using DynamicPolynomials
+using SumOfSquares
+import MultivariateBases as MB
+import SDPLR
+import Hypatia
+import SCS
+import BMSOS
+
+# In this tutorial, we show how to use a different polynomial basis
+# for enforcing the equality between the polynomial and its Sum-of-Squares decomposition.
+
+@polyvar x
+p = x^4 - 4x^3 - 2x^2 + 12x + 3
+
+scs = SCS.Optimizer
+sdplr = optimizer_with_attributes(SDPLR.Optimizer, "maxrank" => (m, n) -> 4)
+hypatia = Hypatia.Optimizer
+bmsos = BMSOS.Optimizer
+function test(solver, feas::Bool)
+    model = Model(solver)
+    set_silent(model)
+    if feas
+        γ = -6
+    else
+        @variable(model, γ)
+        @objective(model, Max, γ)
+    end
+    @constraint(model, p - γ in SOSCone(), zero_basis = BoxSampling([-1], [1]))
+    optimize!(model)
+    @test primal_status == MOI.FEASIBLE_POINT
+    if !feasible
+        @test value(γ) ≈ -6 rtol=1e-4
+    end
+end
+test(scs)
+test(sdplr)
+test(hypatia)
+test(bmsos, true)
+
+import Random
+import TrigPolys
+# See https://codeocean.com/capsule/8311748/tree/v1
+function random_positive_poly(n; tol=1e-5)
+    Random.seed!(0)
+    p = TrigPolys.random_trig_poly(n)
+    p - minimum(TrigPolys.evaluate(TrigPolys.pad_to(p, 10000000))) + n * tol
+    a = zeros(2n + 1)
+    a[1] = p.a0
+    a[2:2:2n] = p.ac
+    a[3:2:(2n+1)] = p.as
+    return MB.algebra_element(
+        a,
+        MB.SubBasis{MB.Trigonometric}(monomials(x, 0:2n)),
+    )
+end
+random_positive_poly(20)
+
+function test_rand(solver, d, B)
+    model = Model(solver)
+    set_silent(model)
+    p = MB.algebra_element(rand(2d+1), MB.SubBasis{B}(monomials(x, 0:2d)))
+    @constraint(model, p in SOSCone(), zero_basis = BoxSampling([-1], [1]))
+    optimize!(model)
+    return solve_time(model)
+end
+
+d = 10
+test_rand(scs, d, MultivariateBases.Trigonometric)
+test_rand(hypatia, d, MultivariateBases.Trigonometric)
+test_rand(sdplr, d, MultivariateBases.Trigonometric)
+test_rand(bmsos, d, MultivariateBases.Trigonometric)

src/Bridges/Variable/Variable.jl

@@ -12,5 +12,6 @@ include("psd2x2.jl")
 include("scaled_diagonally_dominant.jl")
 include("copositive_inner.jl")
 include("kernel.jl")
+include("lowrank.jl")
 
 end

src/Bridges/Variable/kernel.jl

@@ -33,9 +33,11 @@ end
 
 function MOI.Bridges.Variable.supports_constrained_variable(
     ::Type{<:KernelBridge},
-    ::Type{<:SOS.WeightedSOSCone},
-)
-    return true
+    ::Type{<:SOS.WeightedSOSCone{M,B}},
+) where {M,B}
+    # Could be made to work but doesn't work yet so it's best to use
+    # `LowRankBridge` which can then be bridged to classical PSD by MOI's bridge
+    return !(B <: MB.LagrangeBasis)
 end
 
 function MOI.Bridges.added_constrained_variable_types(

src/Bridges/Variable/lowrank.jl

@@ -0,0 +1,79 @@
+struct LowRankBridge{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
+
+import LinearAlgebra
+
+function MOI.Bridges.Variable.bridge_constrained_variable(
+    ::Type{LowRankBridge{T,M}},
+    model::MOI.ModelLike,
+    set::SOS.WeightedSOSCone{M},
+) where {T,M}
+    variables = Vector{Vector{MOI.VariableIndex}}(undef, length(set.gram_bases))
+    constraints = Vector{MOI.ConstraintIndex{MOI.VectorOfVariables}}(undef, length(set.gram_bases))
+    for i in eachindex(set.gram_bases)
+        U = MB.transformation_to(set.gram_bases[i], set.basis)
+        weights = SA.coeffs(set.weights[i], set.basis)
+        variables[i], constraints[i] = MOI.add_constrained_variables(
+            model,
+            MOI.SetWithDotProducts(
+                SOS.matrix_cone(M, length(set.gram_bases[i])),
+                [
+                    MOI.TriangleVectorization(
+                        MOI.LowRankMatrix(
+                            [weights[j]],
+                            reshape(U[j, :], size(U, 2), 1),
+                        )
+                    )
+                    for j in eachindex(set.basis)
+                ],
+            ),
+        )
+    end
+    return KernelBridge{T,M}(
+        [
+            MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(one(T), variables[i][j]) for i in eachindex(set.gram_bases)],
+                zero(T),
+            )
+            for j in eachindex(set.basis)
+        ],
+        variables,
+        constraints,
+        set,
+    )
+end
+
+function MOI.Bridges.Variable.supports_constrained_variable(
+    ::Type{<:LowRankBridge},
+    ::Type{<:SOS.WeightedSOSCone{M,B}},
+) where {M,B}
+    # Could be made to work for non-LagrangeBasis but it's not low rank in
+    # we we'll need a high bridge cost for them so that it's only UnsafeAddMul
+    # if the other bridges are removed
+    return B <: MB.LagrangeBasis
+end
+
+function MOI.Bridges.added_constrained_variable_types(
+    ::Type{LowRankBridge{T,M}},
+) where {T,M}
+    return Tuple{Type}[
+        (MOI.SetWithDotProducts{S[1],MOI.TriangleVectorization{MOI.LowRankMatrix{T}}},)
+        for S in SOS.Bridges.Constraint.constrained_variable_types(M)
+        if S[1] == MOI.PositiveSemidefiniteConeTriangle # FIXME hack
+    ]
+end
+
+function MOI.Bridges.added_constraint_types(::Type{<:LowRankBridge})
+    return Tuple{Type,Type}[]
+end
+
+function MOI.Bridges.Variable.concrete_bridge_type(
+    ::Type{<:LowRankBridge{T}},
+    ::Type{<:SOS.WeightedSOSCone{M}},
+) where {T,M}
+    return LowRankBridge{T,M}
+end

src/variables.jl

@@ -9,10 +9,11 @@ function _bridge_coefficient_type(::Type{SOSPolynomialSet{D,B,C}}) where {D,B,C}
 end
 
 function PolyJuMP.bridges(S::Type{<:WeightedSOSCone})
-    return Tuple{Type,Type}[(
-        Bridges.Variable.KernelBridge,
-        _bridge_coefficient_type(S),
-    )]
+    T = _bridge_coefficient_type(S)
+    return Tuple{Type,Type}[
+        (Bridges.Variable.KernelBridge, T),
+        (Bridges.Variable.LowRankBridge, T),
+    ]
 end
 
 function PolyJuMP.bridges(::Type{<:PositiveSemidefinite2x2ConeTriangle})