From 6a3f1223364c45c4ee0ec42007664c12b0f8a347 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 11 Jul 2024 23:26:43 +0200 Subject: [PATCH 1/6] Sampling --- docs/Project.toml | 3 + .../src/tutorials/Getting started/sampling.jl | 57 +++++++++++++ src/Bridges/Variable/Variable.jl | 1 + src/Bridges/Variable/kernel.jl | 8 +- src/Bridges/Variable/lowrank.jl | 79 +++++++++++++++++++ src/variables.jl | 9 ++- 6 files changed, 150 insertions(+), 7 deletions(-) create mode 100644 docs/src/tutorials/Getting started/sampling.jl create mode 100644 src/Bridges/Variable/lowrank.jl diff --git a/docs/Project.toml b/docs/Project.toml index 9a1959d0c..562a86e28 100644 --- a/docs/Project.toml +++ b/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,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" diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl new file mode 100644 index 000000000..f634f75f1 --- /dev/null +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -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) diff --git a/src/Bridges/Variable/Variable.jl b/src/Bridges/Variable/Variable.jl index 4a154fd2d..36efa95bc 100644 --- a/src/Bridges/Variable/Variable.jl +++ b/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 diff --git a/src/Bridges/Variable/kernel.jl b/src/Bridges/Variable/kernel.jl index cf438f7a4..1b53f1c7a 100644 --- a/src/Bridges/Variable/kernel.jl +++ b/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( diff --git a/src/Bridges/Variable/lowrank.jl b/src/Bridges/Variable/lowrank.jl new file mode 100644 index 000000000..9708abbfe --- /dev/null +++ b/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 diff --git a/src/variables.jl b/src/variables.jl index e3b0eeaf7..c5b3d1f97 100644 --- a/src/variables.jl +++ b/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}) From faf17503ca271dee6638f049e09b99ab60a5918f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 16 Jul 2024 10:04:28 +0200 Subject: [PATCH 2/6] Add chsh example --- .../Noncommutative and Hermitian/chsh.jl | 116 ++++++++++++ src/Bridges/Bridges.jl | 5 + src/Bridges/Constraint/Constraint.jl | 9 + src/Bridges/Constraint/image.jl | 24 ++- src/Bridges/Variable/Variable.jl | 9 + src/Bridges/Variable/lagrange.jl | 178 ++++++++++++++++++ 6 files changed, 328 insertions(+), 13 deletions(-) create mode 100644 docs/src/tutorials/Noncommutative and Hermitian/chsh.jl create mode 100644 src/Bridges/Variable/lagrange.jl diff --git a/docs/src/tutorials/Noncommutative and Hermitian/chsh.jl b/docs/src/tutorials/Noncommutative and Hermitian/chsh.jl new file mode 100644 index 000000000..f691465eb --- /dev/null +++ b/docs/src/tutorials/Noncommutative and Hermitian/chsh.jl @@ -0,0 +1,116 @@ +import QuantumStuff: trace_monoid, Monoids +using StarAlgebras + +M, A, C = trace_monoid(2, 2, A=:A, C=:C) + +RM = let M = M, A = A, C = C, level = 4 + A_l, sizesA = Monoids.wlmetric_ball(A, radius=level) + C_l, sizesC = Monoids.wlmetric_ball(C, radius=level) + + # starAlg(M, 1, half = unique!([a*c for a in A_l for c in C_l])) + + @time words, sizes = Monoids.wlmetric_ball( + unique!([a * c for a in A_l for c in C_l]); + radius=2, + ) + @info "Sizes of generated balls:" (A, C, combined) = (sizesA, sizesC, sizes) + + b = @time StarAlgebras.FixedBasis(words, StarAlgebras.DiracMStructure(*), (UInt32(sizes[1]), UInt32(sizes[1]))) + StarAlgebra(M, b) +end + +A = RM.(A) +C = RM.(C) +chsh = A[1] * C[1] + A[1] * C[2] + A[2] * C[1] - A[2] * C[2] + +import StarAlgebras as SA +struct Full{B} <: SA.ImplicitBasis{B,B} end +Base.getindex(::Full{B}, b::B) where {B} = b +import MultivariateBases as MB +MB.implicit_basis(::SA.FixedBasis{B}) where {B} = Full{B}() +MB.algebra(b::Full{B}) where {B} = SA.StarAlgebra(M, b) +SA.mstructure(::Full) = SA.DiracMStructure(*) + +b = basis(chsh) +import StarAlgebras as SA +f = SA.AlgebraElement( + SA.SparseCoefficients( + [b[k] for (k, _) in SA.nonzero_pairs(coeffs(chsh))], + [v for (_, v) in SA.nonzero_pairs(coeffs(chsh))], + ), + SA.StarAlgebra( + parent(chsh).object, + Full{eltype(b)}() + ), +) +n = size(b.table, 1) +gram_basis = @time StarAlgebras.FixedBasis(b.elts[1:n], StarAlgebras.DiracMStructure(*)); +one(f) +SA.coeffs(f, b) +using SumOfSquares +function SumOfSquares._term_element(a, mono::Monoids.MonoidElement) + SA.AlgebraElement( + SA.SparseCoefficients((mono,), (a,)), + MB.algebra(Full{typeof(mono)}()), + ) +end + +cone = SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( + b, + [gram_basis], + [one(f)], +) +import SCS +scs = optimizer_with_attributes( + SCS.Optimizer, + "eps_abs" => 1e-9, + "eps_rel" => 1e-9, + "max_iters" => 1000, +) + +import Dualization +#model = Model(Dualization.dual_optimizer(scs)) +model = Model(scs) +SumOfSquares.Bridges.add_all_bridges(backend(model).optimizer, Float64) +MOI.Bridges.remove_bridge(backend(model).optimizer, SumOfSquares.Bridges.Constraint.ImageBridge{Float64}) +@variable(model, λ) +@objective(model, Min, λ) +@constraint(model, SA.coeffs(λ * one(f) - f, b) in cone); +optimize!(model) +solution_summary(model) +objective_value(model) ≈ 2√2 +function _add!(f, psd, model, F, S) + append!(psd, [ + f(MOI.get(model, MOI.ConstraintSet(), ci)) + for ci in MOI.get(model, MOI.ListOfConstraintIndices{F,S}()) + ], + ) +end +function summary(model) + free = MOI.get(model, MOI.NumberOfVariables()) + psd = Int[] + F = MOI.VectorOfVariables + S = MOI.PositiveSemidefiniteConeTriangle + _add!(MOI.side_dimension, psd, model, F, S) + S = SCS.ScaledPSDCone + _add!(MOI.side_dimension, psd, model, F, S) + free -= sum(psd, init = 0) do d + div(d * (d + 1), 2) + end + F = MOI.VectorAffineFunction{Float64} + S = MOI.PositiveSemidefiniteConeTriangle + _add!(MOI.side_dimension, psd, model, F, S) + S = SCS.ScaledPSDCone + _add!(MOI.side_dimension, psd, model, F, S) + eq = Int[] + F = MOI.VectorAffineFunction{Float64} + S = MOI.Zeros + _add!(MOI.dimension, eq, model, F, S) + F = MOI.ScalarAffineFunction{Float64} + S = MOI.EqualTo{Float64} + _add!(MOI.dimension, eq, model, F, S) + return free, psd, sum(eq, init = 0) +end +summary(backend(model).optimizer.model) +summary(backend(model).optimizer.model.optimizer.dual_problem.dual_model.model) +print_active_bridges(model) diff --git a/src/Bridges/Bridges.jl b/src/Bridges/Bridges.jl index c03efeb62..ae6f820f8 100644 --- a/src/Bridges/Bridges.jl +++ b/src/Bridges/Bridges.jl @@ -6,6 +6,11 @@ import SumOfSquares as SOS include("Variable/Variable.jl") include("Constraint/Constraint.jl") +function add_all_bridges(model, ::Type{T}) where {T} + Variable.add_all_bridges(model, T) + Constraint.add_all_bridges(model, T) +end + function MOI.get( model::MOI.ModelLike, attr::Union{ diff --git a/src/Bridges/Constraint/Constraint.jl b/src/Bridges/Constraint/Constraint.jl index 3ea7a7542..27459f459 100644 --- a/src/Bridges/Constraint/Constraint.jl +++ b/src/Bridges/Constraint/Constraint.jl @@ -26,4 +26,13 @@ include("image.jl") include("sos_polynomial.jl") include("sos_polynomial_in_semialgebraic_set.jl") +function add_all_bridges(model, ::Type{T}) where {T} + MOI.Bridges.add_bridge(model, EmptyBridge{T}) + MOI.Bridges.add_bridge(model, PositiveSemidefinite2x2Bridge{T}) + MOI.Bridges.add_bridge(model, DiagonallyDominantBridge{T}) + MOI.Bridges.add_bridge(model, ImageBridge{T}) + MOI.Bridges.add_bridge(model, SOSPolynomialBridge{T}) + MOI.Bridges.add_bridge(model, SOSPolynomialInSemialgebraicSetBridge{T}) +end + end diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index b7281ff86..954dd4f6b 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -81,17 +81,18 @@ function MOI.Bridges.Constraint.bridge_constraint( @assert MOI.output_dimension(g) == length(set.basis) scalars = MOI.Utilities.scalarize(g) k = 0 - found = Dict{eltype(set.basis.monomials),Int}() + found = Dict{eltype(set.basis),Int}() first = Union{Nothing,Int}[nothing for _ in eachindex(scalars)] variables = MOI.VariableIndex[] constraints = MOI.ConstraintIndex{F}[] for (gram_basis, weight) in zip(set.gram_bases, set.weights) + # TODO don't ignore weight cone = SOS.matrix_cone(M, length(gram_basis)) f = MOI.Utilities.zero_with_output_dimension(F, MOI.dimension(cone)) - for j in eachindex(gram_basis.monomials) + for j in eachindex(gram_basis) for i in 1:j k += 1 - mono = gram_basis.monomials[i] * gram_basis.monomials[j] + mono = SA.star(gram_basis[i]) * gram_basis[j] is_diag = i == j if haskey(found, mono) var = MOI.add_variable(model) @@ -119,8 +120,9 @@ function MOI.Bridges.Constraint.bridge_constraint( MOI.Utilities.operate_output_index!(-, T, k, f, var) else found[mono] = k - t = MB.monomial_index(set.basis, mono) - if !isnothing(t) + @show mono, set.basis[mono] + if mono in set.basis + t = set.basis[mono] first[t] = k if is_diag MOI.Utilities.operate_output_index!( @@ -139,6 +141,8 @@ function MOI.Bridges.Constraint.bridge_constraint( inv(T(2)) * scalars[t], ) end + else + @warn("$mono not in basis") end end end @@ -167,14 +171,8 @@ end function MOI.supports_constraint( ::Type{ImageBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, - ::Type{ - <:SOS.WeightedSOSCone{ - M, - <:MB.SubBasis{MB.Monomial}, - <:MB.SubBasis{MB.Monomial}, - }, - }, -) where {T,M} + ::Type{<:SOS.WeightedSOSCone}, +) where {T} return true end diff --git a/src/Bridges/Variable/Variable.jl b/src/Bridges/Variable/Variable.jl index 4a154fd2d..7893c3dae 100644 --- a/src/Bridges/Variable/Variable.jl +++ b/src/Bridges/Variable/Variable.jl @@ -13,4 +13,13 @@ include("scaled_diagonally_dominant.jl") include("copositive_inner.jl") include("kernel.jl") +function add_all_bridges(model, ::Type{T}) where {T} + MOI.Bridges.add_bridge(model, PositiveSemidefinite2x2Bridge{T}) + MOI.Bridges.add_bridge(model, ScaledDiagonallyDominantBridge{T}) + MOI.Bridges.add_bridge(model, CopositiveInnerBridge{T}) + MOI.Bridges.add_bridge(model, KernelBridge{T}) + MOI.Bridges.add_bridge(model, LowRankBridge{T}) + return +end + end diff --git a/src/Bridges/Variable/lagrange.jl b/src/Bridges/Variable/lagrange.jl new file mode 100644 index 000000000..ae8de253a --- /dev/null +++ b/src/Bridges/Variable/lagrange.jl @@ -0,0 +1,178 @@ +struct LagrangeBridge{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,M}}, + model::MOI.ModelLike, + set::SOS.WeightedSOSCone{M}, +) where {T,M} + variables = Vector{MOI.VariableIndex}[] + constraints = MOI.ConstraintIndex{MOI.VectorOfVariables}[] + acc = zero( + MOI.ScalarAffineFunction{T}, + MB.algebra(MB.implicit_basis(set.basis)), + ) + 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) + MA.operate!(SA.UnsafeAddMul(*), acc, gram, weight) + end + MA.operate!(SA.canonical, SA.coeffs(acc)) + return KernelBridge{T,M}( + SA.coeffs(acc, set.basis), + 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 + +function MOI.Bridges.Variable.concrete_bridge_type( + ::Type{<:KernelBridge{T}}, + ::Type{<:SOS.WeightedSOSCone{M}}, +) where {T,M} + return KernelBridge{T,M} +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_variables(bridge.affine[i.value]) do vi + return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi) + end +end + +function MOI.get( + model::MOI.ModelLike, + attr::SOS.GramMatrixAttribute, + bridge::KernelBridge{T,M}, +) where {T,M} + SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) + return SOS.build_gram_matrix( + convert( + Vector{T}, + MOI.get( + model, + MOI.VariablePrimal(attr.result_index), + bridge.variables[attr.multiplier_index], + ), + ), + bridge.set.gram_bases[attr.multiplier_index], + M, + T, + ) +end + +function MOI.get( + model::MOI.ModelLike, + attr::SOS.MomentMatrixAttribute, + bridge::KernelBridge{T,M}, +) where {T,M} + SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) + return SOS.build_moment_matrix( + convert( + Vector{T}, + MOI.get( + model, + MOI.ConstraintDual(attr.result_index), + bridge.constraints[attr.multiplier_index], + ), + ), + bridge.set.gram_bases[attr.multiplier_index], + ) +end + +function MOI.Bridges.bridged_function( + bridge::KernelBridge, + i::MOI.Bridges.IndexInVector, +) + return bridge.affine[i.value] +end + +function MOI.Bridges.Variable.unbridged_map( + ::KernelBridge{T}, + ::Vector{MOI.VariableIndex}, +) where {T} + return nothing +end From ad05db5dd122f5bb173b24dc4b3f103a4cbf2195 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 16 Jul 2024 10:10:21 +0200 Subject: [PATCH 3/6] Add random --- .../src/tutorials/Getting started/sampling.jl | 25 +++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl index f634f75f1..51a773087 100644 --- a/docs/src/tutorials/Getting started/sampling.jl +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -44,6 +44,24 @@ 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) @@ -53,5 +71,8 @@ function test_rand(solver, d, B) return solve_time(model) end -test_rand(scs, 2, MultivariateBases.Trigonometric) -test_rand(bmsos, 4, MultivariateBases.Trigonometric) +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) From 3a79a65260858cf1d9c0d7be4ea63e3d8908e9f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 16 Jul 2024 10:10:30 +0200 Subject: [PATCH 4/6] Add TrigPolys --- docs/Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/Project.toml b/docs/Project.toml index 562a86e28..34fae70a2 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -35,6 +35,7 @@ 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] From 2281da7d096784454e48dd570610db4d98eb5f28 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 19 Aug 2024 18:00:31 +0200 Subject: [PATCH 5/6] up --- docs/Project.toml | 2 + .../src/tutorials/Getting started/sampling.jl | 91 ++++++++++++++++--- src/Bridges/Constraint/image.jl | 1 - src/Bridges/Variable/lagrange.jl | 38 ++++---- src/Bridges/Variable/lowrank.jl | 16 +++- 5 files changed, 113 insertions(+), 35 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 34fae70a2..bd2c769ff 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,6 +5,8 @@ CSDP = "0a46da34-8e4b-519e-b418-48813639ff34" Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b" +DSDP = "2714ae6b-e930-5b4e-9c21-d0bacf577842" +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl index 51a773087..47edac90c 100644 --- a/docs/src/tutorials/Getting started/sampling.jl +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -23,6 +23,8 @@ scs = SCS.Optimizer sdplr = optimizer_with_attributes(SDPLR.Optimizer, "maxrank" => (m, n) -> 4) hypatia = Hypatia.Optimizer bmsos = BMSOS.Optimizer +import DSDP +dsdp = DSDP.Optimizer function test(solver, feas::Bool) model = Model(solver) set_silent(model) @@ -34,15 +36,22 @@ function test(solver, feas::Bool) 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 + @test primal_status(model) == MOI.FEASIBLE_POINT + @show objective_value(model) + @show value(γ) + if !feas + if !isapprox(value(γ), -6, rtol=1e-3) + @warn("$(value(γ)) != -6") + end + #@test value(γ) ≈ -6 rtol=1e-4 end + return model end -test(scs) -test(sdplr) -test(hypatia) -test(bmsos, true) +model = test(scs, false); +model = test(sdplr, false); +test(hypatia, false); +#test(bmsos, true) +model = test(dsdp, false); import Random import TrigPolys @@ -50,7 +59,9 @@ import TrigPolys 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 + #N = 10000000 + N = 1000000 + p - minimum(TrigPolys.evaluate(TrigPolys.pad_to(p, N))) + n * tol a = zeros(2n + 1) a[1] = p.a0 a[2:2:2n] = p.ac @@ -65,14 +76,66 @@ 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))) + p = if B == MB.Trigonometric + random_positive_poly(d) + else + MB.algebra_element(rand(2d+1), MB.SubBasis{B}(monomials(x, 0:2d))) + end @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) +using DataFrames +df = DataFrame(solver=String[], degree=Int[], time=Float64[]) + +function timing(solver, d, dual::Bool = false) + name = MOI.get(MOI.instantiate(solver), MOI.SolverName()) + if dual + solver = dual_optimizer(solver) + end + time = test_rand(solver, d, MB.Trigonometric) + push!(df, (name, d, time)) +end + +using Dualization +using MosekTools +mosek = Mosek.Optimizer + +d = 400 +timing(bmsos, d) +timing(hypatia, d) +timing(dsdp, d) +timing(mosek, d) +#timing(sdplr, d) +#timing(dual_optimizer(scs), d) +#timing(scs, d, true) +#timing(mosek, d, true) + +using Printf + +function table(degs, solvers) + print("| |") + for solver in solvers + print(" $solver |") + end + println() + for _ in 0:length(solvers) + print("|-----") + end + println("|") + for deg in degs + print("| ", deg, " |") + for solver in solvers + times = subset(df, :solver => c -> c .== Ref(solver), :degree => d -> d .== deg).time + if isempty(times) + print(" |") + else + @printf(" %.3e |", minimum(times)) + end + end + println() + end +end + +table([100, 200, 300, 400, 500, 800], ["BMSOS", "Hypatia", "DSDP", "Mosek"]) diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index 954dd4f6b..1cecb628a 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -120,7 +120,6 @@ function MOI.Bridges.Constraint.bridge_constraint( MOI.Utilities.operate_output_index!(-, T, k, f, var) else found[mono] = k - @show mono, set.basis[mono] if mono in set.basis t = set.basis[mono] first[t] = k diff --git a/src/Bridges/Variable/lagrange.jl b/src/Bridges/Variable/lagrange.jl index ae8de253a..9ffe70068 100644 --- a/src/Bridges/Variable/lagrange.jl +++ b/src/Bridges/Variable/lagrange.jl @@ -6,7 +6,7 @@ struct LagrangeBridge{T,M} <: MOI.Bridges.Variable.AbstractBridge end function MOI.Bridges.Variable.bridge_constrained_variable( - ::Type{KernelBridge{T,M}}, + ::Type{LagrangeBridge{T,M}}, model::MOI.ModelLike, set::SOS.WeightedSOSCone{M}, ) where {T,M} @@ -23,7 +23,7 @@ function MOI.Bridges.Variable.bridge_constrained_variable( MA.operate!(SA.UnsafeAddMul(*), acc, gram, weight) end MA.operate!(SA.canonical, SA.coeffs(acc)) - return KernelBridge{T,M}( + return LagrangeBridge{T,M}( SA.coeffs(acc, set.basis), variables, constraints, @@ -32,40 +32,40 @@ function MOI.Bridges.Variable.bridge_constrained_variable( end function MOI.Bridges.Variable.supports_constrained_variable( - ::Type{<:KernelBridge}, + ::Type{<:LagrangeBridge}, ::Type{<:SOS.WeightedSOSCone}, ) return true end function MOI.Bridges.added_constrained_variable_types( - ::Type{KernelBridge{T,M}}, + ::Type{LagrangeBridge{T,M}}, ) where {T,M} return SOS.Bridges.Constraint.constrained_variable_types(M) end -function MOI.Bridges.added_constraint_types(::Type{<:KernelBridge}) +function MOI.Bridges.added_constraint_types(::Type{<:LagrangeBridge}) return Tuple{Type,Type}[] end function MOI.Bridges.Variable.concrete_bridge_type( - ::Type{<:KernelBridge{T}}, + ::Type{<:LagrangeBridge{T}}, ::Type{<:SOS.WeightedSOSCone{M}}, ) where {T,M} - return KernelBridge{T,M} + return LagrangeBridge{T,M} end # Attributes, Bridge acting as a model -function MOI.get(bridge::KernelBridge, ::MOI.NumberOfVariables) +function MOI.get(bridge::LagrangeBridge, ::MOI.NumberOfVariables) return sum(length, bridge.variables) end -function MOI.get(bridge::KernelBridge, ::MOI.ListOfVariableIndices) +function MOI.get(bridge::LagrangeBridge, ::MOI.ListOfVariableIndices) return reduce(vcat, bridge.variables) end function MOI.get( - bridge::KernelBridge, + bridge::LagrangeBridge, ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S}, ) where {S<:MOI.AbstractVectorSet} return count(bridge.constraints) do ci @@ -74,7 +74,7 @@ function MOI.get( end function MOI.get( - bridge::KernelBridge, + bridge::LagrangeBridge, ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}, ) where {S} return [ @@ -84,7 +84,7 @@ function MOI.get( end # Indices -function MOI.delete(model::MOI.ModelLike, bridge::KernelBridge) +function MOI.delete(model::MOI.ModelLike, bridge::LagrangeBridge) for vars in bridge.variables MOI.delete(model, vars) end @@ -93,14 +93,14 @@ end # Attributes, Bridge acting as a constraint -function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::KernelBridge) +function MOI.get(::MOI.ModelLike, ::MOI.ConstraintSet, bridge::LagrangeBridge) return bridge.set end function MOI.get( model::MOI.ModelLike, attr::MOI.ConstraintPrimal, - bridge::KernelBridge, + bridge::LagrangeBridge, ) return [ MOI.get( @@ -115,7 +115,7 @@ end function MOI.get( model::MOI.ModelLike, attr::MOI.VariablePrimal, - bridge::KernelBridge, + bridge::LagrangeBridge, i::MOI.Bridges.IndexInVector, ) return MOI.Utilities.eval_variables(bridge.affine[i.value]) do vi @@ -126,7 +126,7 @@ end function MOI.get( model::MOI.ModelLike, attr::SOS.GramMatrixAttribute, - bridge::KernelBridge{T,M}, + bridge::LagrangeBridge{T,M}, ) where {T,M} SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) return SOS.build_gram_matrix( @@ -147,7 +147,7 @@ end function MOI.get( model::MOI.ModelLike, attr::SOS.MomentMatrixAttribute, - bridge::KernelBridge{T,M}, + bridge::LagrangeBridge{T,M}, ) where {T,M} SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) return SOS.build_moment_matrix( @@ -164,14 +164,14 @@ function MOI.get( end function MOI.Bridges.bridged_function( - bridge::KernelBridge, + bridge::LagrangeBridge, i::MOI.Bridges.IndexInVector, ) return bridge.affine[i.value] end function MOI.Bridges.Variable.unbridged_map( - ::KernelBridge{T}, + ::LagrangeBridge{T}, ::Vector{MOI.VariableIndex}, ) where {T} return nothing diff --git a/src/Bridges/Variable/lowrank.jl b/src/Bridges/Variable/lowrank.jl index 9708abbfe..9a53a5bd1 100644 --- a/src/Bridges/Variable/lowrank.jl +++ b/src/Bridges/Variable/lowrank.jl @@ -33,7 +33,7 @@ function MOI.Bridges.Variable.bridge_constrained_variable( ), ) end - return KernelBridge{T,M}( + return LowRankBridge{T,M}( [ MOI.ScalarAffineFunction( [MOI.ScalarAffineTerm(one(T), variables[i][j]) for i in eachindex(set.gram_bases)], @@ -77,3 +77,17 @@ function MOI.Bridges.Variable.concrete_bridge_type( ) where {T,M} return LowRankBridge{T,M} end + +function MOI.Bridges.bridged_function( + bridge::LowRankBridge, + i::MOI.Bridges.IndexInVector, +) + return bridge.affine[i.value] +end + +function MOI.Bridges.Variable.unbridged_map( + ::LowRankBridge, + ::Vector{MOI.VariableIndex}, +) + return nothing +end From 9f22461c868c586a462c8ac837d097bf3c323d10 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 19 Aug 2024 18:01:19 +0200 Subject: [PATCH 6/6] Delete lagrange --- src/Bridges/Variable/lagrange.jl | 178 ------------------------------- 1 file changed, 178 deletions(-) delete mode 100644 src/Bridges/Variable/lagrange.jl diff --git a/src/Bridges/Variable/lagrange.jl b/src/Bridges/Variable/lagrange.jl deleted file mode 100644 index 9ffe70068..000000000 --- a/src/Bridges/Variable/lagrange.jl +++ /dev/null @@ -1,178 +0,0 @@ -struct LagrangeBridge{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{LagrangeBridge{T,M}}, - model::MOI.ModelLike, - set::SOS.WeightedSOSCone{M}, -) where {T,M} - variables = Vector{MOI.VariableIndex}[] - constraints = MOI.ConstraintIndex{MOI.VectorOfVariables}[] - acc = zero( - MOI.ScalarAffineFunction{T}, - MB.algebra(MB.implicit_basis(set.basis)), - ) - 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) - MA.operate!(SA.UnsafeAddMul(*), acc, gram, weight) - end - MA.operate!(SA.canonical, SA.coeffs(acc)) - return LagrangeBridge{T,M}( - SA.coeffs(acc, set.basis), - variables, - constraints, - set, - ) -end - -function MOI.Bridges.Variable.supports_constrained_variable( - ::Type{<:LagrangeBridge}, - ::Type{<:SOS.WeightedSOSCone}, -) - return true -end - -function MOI.Bridges.added_constrained_variable_types( - ::Type{LagrangeBridge{T,M}}, -) where {T,M} - return SOS.Bridges.Constraint.constrained_variable_types(M) -end - -function MOI.Bridges.added_constraint_types(::Type{<:LagrangeBridge}) - return Tuple{Type,Type}[] -end - -function MOI.Bridges.Variable.concrete_bridge_type( - ::Type{<:LagrangeBridge{T}}, - ::Type{<:SOS.WeightedSOSCone{M}}, -) where {T,M} - return LagrangeBridge{T,M} -end - -# Attributes, Bridge acting as a model -function MOI.get(bridge::LagrangeBridge, ::MOI.NumberOfVariables) - return sum(length, bridge.variables) -end - -function MOI.get(bridge::LagrangeBridge, ::MOI.ListOfVariableIndices) - return reduce(vcat, bridge.variables) -end - -function MOI.get( - bridge::LagrangeBridge, - ::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::LagrangeBridge, - ::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::LagrangeBridge) - 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::LagrangeBridge) - return bridge.set -end - -function MOI.get( - model::MOI.ModelLike, - attr::MOI.ConstraintPrimal, - bridge::LagrangeBridge, -) - 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::LagrangeBridge, - i::MOI.Bridges.IndexInVector, -) - return MOI.Utilities.eval_variables(bridge.affine[i.value]) do vi - return MOI.get(model, MOI.VariablePrimal(attr.result_index), vi) - end -end - -function MOI.get( - model::MOI.ModelLike, - attr::SOS.GramMatrixAttribute, - bridge::LagrangeBridge{T,M}, -) where {T,M} - SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) - return SOS.build_gram_matrix( - convert( - Vector{T}, - MOI.get( - model, - MOI.VariablePrimal(attr.result_index), - bridge.variables[attr.multiplier_index], - ), - ), - bridge.set.gram_bases[attr.multiplier_index], - M, - T, - ) -end - -function MOI.get( - model::MOI.ModelLike, - attr::SOS.MomentMatrixAttribute, - bridge::LagrangeBridge{T,M}, -) where {T,M} - SOS.check_multiplier_index_bounds(attr, eachindex(bridge.constraints)) - return SOS.build_moment_matrix( - convert( - Vector{T}, - MOI.get( - model, - MOI.ConstraintDual(attr.result_index), - bridge.constraints[attr.multiplier_index], - ), - ), - bridge.set.gram_bases[attr.multiplier_index], - ) -end - -function MOI.Bridges.bridged_function( - bridge::LagrangeBridge, - i::MOI.Bridges.IndexInVector, -) - return bridge.affine[i.value] -end - -function MOI.Bridges.Variable.unbridged_map( - ::LagrangeBridge{T}, - ::Vector{MOI.VariableIndex}, -) where {T} - return nothing -end