Skip to content

Commit

Permalink
Sampling
Browse files Browse the repository at this point in the history
  • Loading branch information
blegat committed Jul 11, 2024
1 parent 4af4b3e commit 6a3f122
Show file tree
Hide file tree
Showing 6 changed files with 150 additions and 7 deletions.
3 changes: 3 additions & 0 deletions docs/Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"
Expand All @@ -28,6 +30,7 @@ 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"
Expand Down
57 changes: 57 additions & 0 deletions docs/src/tutorials/Getting started/sampling.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
# # 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)

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

test_rand(scs, 2, MultivariateBases.Trigonometric)
test_rand(bmsos, 4, MultivariateBases.Trigonometric)
1 change: 1 addition & 0 deletions src/Bridges/Variable/Variable.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,5 +12,6 @@ include("psd2x2.jl")
include("scaled_diagonally_dominant.jl")
include("copositive_inner.jl")
include("kernel.jl")
include("lowrank.jl")

end
8 changes: 5 additions & 3 deletions src/Bridges/Variable/kernel.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
79 changes: 79 additions & 0 deletions src/Bridges/Variable/lowrank.jl
Original file line number Diff line number Diff line change
@@ -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
9 changes: 5 additions & 4 deletions src/variables.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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})
Expand Down

0 comments on commit 6a3f122

Please sign in to comment.