Sampling

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]

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"])

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)

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{

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

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
 

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