diff --git a/docs/Project.toml b/docs/Project.toml index 9a1959d0c..bd2c769ff 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,9 +1,12 @@ [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" 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" @@ -12,6 +15,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 +32,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] diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl new file mode 100644 index 000000000..47edac90c --- /dev/null +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -0,0 +1,141 @@ +# # 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 +import DSDP +dsdp = DSDP.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(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 +model = test(scs, false); +model = test(sdplr, false); +test(hypatia, false); +#test(bmsos, true) +model = test(dsdp, false); + +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) + #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 + 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 = 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 + +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/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..1cecb628a 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,8 @@ 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) + if mono in set.basis + t = set.basis[mono] first[t] = k if is_diag MOI.Utilities.operate_output_index!( @@ -139,6 +140,8 @@ function MOI.Bridges.Constraint.bridge_constraint( inv(T(2)) * scalars[t], ) end + else + @warn("$mono not in basis") end end end @@ -167,14 +170,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..3fbb097b9 100644 --- a/src/Bridges/Variable/Variable.jl +++ b/src/Bridges/Variable/Variable.jl @@ -12,5 +12,15 @@ include("psd2x2.jl") include("scaled_diagonally_dominant.jl") include("copositive_inner.jl") include("kernel.jl") +include("lowrank.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/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..9a53a5bd1 --- /dev/null +++ b/src/Bridges/Variable/lowrank.jl @@ -0,0 +1,93 @@ +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 LowRankBridge{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 + +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 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})