diff --git a/src/PolyJuMP.jl b/src/PolyJuMP.jl index 8b5452e..8b017bf 100644 --- a/src/PolyJuMP.jl +++ b/src/PolyJuMP.jl @@ -34,5 +34,6 @@ include("default.jl") include("model.jl") include("KKT/KKT.jl") include("QCQP/QCQP.jl") +include("RelativeEntropy/RelativeEntropy.jl") end # module diff --git a/src/RelativeEntropy/RelativeEntropy.jl b/src/RelativeEntropy/RelativeEntropy.jl new file mode 100644 index 0000000..b0fef8c --- /dev/null +++ b/src/RelativeEntropy/RelativeEntropy.jl @@ -0,0 +1,214 @@ +module RelativeEntropy + +import MutableArithmetics as MA +import MultivariatePolynomials as MP +import MathOptInterface as MOI +import JuMP +import PolyJuMP + +abstract type AbstractAGECone <: MOI.AbstractVectorSet end + +""" + struct SignomialSAGECone <: MOI.AbstractVectorSet + α::Matrix{Int} + end + +**S**ums of **A**M/**G**M **E**xponential for signomials. +""" +struct SignomialSAGECone <: AbstractAGECone + α::Matrix{Int} +end + +""" + struct PolynomialSAGECone <: MOI.AbstractVectorSet + α::Matrix{Int} + end + +**S**ums of **A**M/**G**M **E**xponential for polynomials. +""" +struct PolynomialSAGECone <: AbstractAGECone + α::Matrix{Int} +end + +struct SignomialAGECone <: AbstractAGECone + α::Matrix{Int} + k::Int +end + +struct PolynomialAGECone <: AbstractAGECone + α::Matrix{Int} + k::Int +end + +MOI.dimension(set::AbstractAGECone) = size(set.α, 1) +Base.copy(set::AbstractAGECone) = set + +function _exponents_matrix(monos) + α = Matrix{Int}(undef, length(monos), MP.nvariables(monos)) + for (i, mono) in enumerate(monos) + exp = MP.exponents(mono) + for j in eachindex(exp) + α[i, j] = exp[j] + end + end + return α +end + +struct SignomialSAGESet <: PolyJuMP.PolynomialSet end +function JuMP.reshape_set(::SignomialSAGECone, ::PolyJuMP.PolynomialShape) + return SignomialSAGESet() +end +function JuMP.moi_set(::SignomialSAGESet, monos) + return SignomialSAGECone(_exponents_matrix(monos)) +end + +struct PolynomialSAGESet <: PolyJuMP.PolynomialSet end +function JuMP.reshape_set(::PolynomialSAGECone, ::PolyJuMP.PolynomialShape) + return PolynomialSAGESet() +end +function JuMP.moi_set(::PolynomialSAGESet, monos) + return PolynomialSAGECone(_exponents_matrix(monos)) +end + +struct SignomialAGESet{MT<:MP.AbstractMonomial} <: PolyJuMP.PolynomialSet + monomial::MT +end +function JuMP.reshape_set( + set::SignomialAGECone, + shape::PolyJuMP.PolynomialShape, +) + return SignomialAGESet(shape.monomials[set.k]) +end +function JuMP.moi_set(set::SignomialAGESet, monos) + k = findfirst(isequal(set.monomial), monos) + return SignomialAGECone(_exponents_matrix(monos), k) +end + +struct PolynomialAGESet{MT<:MP.AbstractMonomial} <: PolyJuMP.PolynomialSet + monomial::MT +end +function JuMP.reshape_set( + set::PolynomialAGECone, + shape::PolyJuMP.PolynomialShape, +) + return PolynomialAGESet(shape.monomials[set.k]) +end +function JuMP.moi_set(set::PolynomialAGESet, monos) + k = findfirst(isequal(set.monomial), monos) + return PolynomialAGECone(_exponents_matrix(monos), k) +end + +function setdefaults!(data::PolyJuMP.Data) + PolyJuMP.setdefault!(data, PolyJuMP.NonNegPoly, PolynomialSAGESet) + return +end + +function JuMP.build_constraint( + _error::Function, + p, + set::Union{ + SignomialSAGESet, + PolynomialSAGESet, + SignomialAGESet, + PolynomialAGESet, + }; + kws..., +) + coefs = PolyJuMP.non_constant_coefficients(p) + monos = MP.monomials(p) + cone = JuMP.moi_set(set, monos) + shape = PolyJuMP.PolynomialShape(monos) + return PolyJuMP.bridgeable( + JuMP.VectorConstraint(coefs, cone, shape), + JuMP.moi_function_type(typeof(coefs)), + typeof(cone), + ) +end + +const APL{T} = MP.AbstractPolynomialLike{T} + +""" + struct Decomposition{T, PT} + +Represents a SAGE decomposition. +""" +struct Decomposition{T,P<:APL{T}} <: APL{T} + polynomials::Vector{P} + function Decomposition(ps::Vector{P}) where {T,P<:APL{T}} + return new{T,P}(ps) + end +end + +function Base.show(io::IO, p::Decomposition) + for (i, q) in enumerate(p.polynomials) + print(io, "(") + print(io, q) + print(io, ")") + if i != length(p.polynomials) + print(io, " + ") + end + end +end + +function MP.polynomial(d::Decomposition) + return sum(d.polynomials) +end + +""" + struct DecompositionAttribute{T} <: MOI.AbstractConstraintAttribute + tol::T + end + +A constraint attribute for the [`Decomposition`](@ref) of a constraint. +""" +struct DecompositionAttribute{T} <: MOI.AbstractConstraintAttribute + tol::T + result_index::Int +end +function DecompositionAttribute(tol::Real) + return DecompositionAttribute(tol, 1) +end + +function decomposition( + con_ref::JuMP.ConstraintRef; + tol::Real, + result_index::Int = 1, +) + monos = con_ref.shape.monomials + attr = DecompositionAttribute(tol, result_index) + return Decomposition([ + MP.polynomial(a, monos) for + a in MOI.get(JuMP.owner_model(con_ref), attr, con_ref) + ]) +end + +MOI.is_set_by_optimize(::DecompositionAttribute) = true + +include("bridges/sage.jl") + +function PolyJuMP.bridges( + F::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialSAGECone}, +) + return [(SAGEBridge, PolyJuMP._coef_type(F))] +end + +include("bridges/age.jl") + +function PolyJuMP.bridges( + F::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialAGECone}, +) + return [(AGEBridge, PolyJuMP._coef_type(F))] +end + +include("bridges/signomial.jl") + +function PolyJuMP.bridges( + F::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:Union{PolynomialSAGECone,PolynomialAGECone}}, +) + return [(SignomialBridge, PolyJuMP._coef_type(F))] +end + +end diff --git a/src/RelativeEntropy/bridges/age.jl b/src/RelativeEntropy/bridges/age.jl new file mode 100644 index 0000000..b8da176 --- /dev/null +++ b/src/RelativeEntropy/bridges/age.jl @@ -0,0 +1,116 @@ +""" + AGEBridge{T,F,G,H} <: MOI.Bridges.Constraint.AbstractBridge + +The nonnegativity of `x ≥ 0` in +``` +∑ ci x^αi ≥ -c0 x^α0 +``` +can be reformulated as +``` +∑ ci exp(αi'y) ≥ -β exp(α0'y) +``` +In any case, it is shown to be equivalent to +``` +∃ ν ≥ 0 s.t. D(ν, exp(1)*c) ≤ β, ∑ νi αi = α0 ∑ νi [CP16, (3)] +``` +where `N(ν, λ) = ∑ νj log(νj/λj)` is the relative entropy function. +The constant `exp(1)` can also be moved out of `D` into +``` +∃ ν ≥ 0 s.t. D(ν, c) - ∑ νi ≤ β, ∑ νi αi = α0 ∑ νi [MCW21, (2)] +``` +The relative entropy cone in MOI is `(u, v, w)` such that `D(w, v) ≥ u`. +Therefore, we can either encode `(β, exp(1)*c, ν)` [CP16, (3)] or +`(β + ∑ νi, c, ν)` [MCW21, (2)]. +In this bridge, we use the second formulation. + +!!! note + A direct application of the Arithmetic-Geometric mean inequality gives + ``` + ∃ ν ≥ 0 s.t. D(ν, exp(1)*c) ≤ -log(-β), ∑ νi αi = α0, ∑ νi = 1 [CP16, (4)] + ``` + which is not jointly convex in (ν, c, β). + The key to get the convex formulation [CP16, (3)] or [MCW21, (2)] instead is to + use the convex conjugacy between the exponential and the negative entropy functions [CP16, (iv)]. + +[CP16] Chandrasekaran, Venkat, and Parikshit Shah. +"Relative entropy relaxations for signomial optimization." +SIAM Journal on Optimization 26.2 (2016): 1147-1173. +[MCW21] Murray, Riley, Venkat Chandrasekaran, and Adam Wierman. +"Signomial and polynomial optimization via relative entropy and partial dualization." +Mathematical Programming Computation 13 (2021): 257-295. +https://arxiv.org/pdf/1907.00814.pdf +""" +struct AGEBridge{T,F,G,H} <: MOI.Bridges.Constraint.AbstractBridge + k::Int + equality_constraints::Vector{MOI.ConstraintIndex{F,MOI.EqualTo{T}}} + relative_entropy_constraint::MOI.ConstraintIndex{G,MOI.RelativeEntropyCone} +end + +function MOI.Bridges.Constraint.bridge_constraint( + ::Type{AGEBridge{T,F,G,H}}, + model, + func::H, + set::SignomialAGECone, +) where {T,F,G,H} + m = size(set.α, 1) + ν = MOI.add_variables(model, m - 1) + sumν = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(one(T), ν), zero(T)) + ceq = map(1:size(set.α, 2)) do var + f = zero(typeof(sumν)) + j = 0 + for i in 1:m + if i == set.k + MA.sub_mul!!(f, convert(T, set.α[i, var]), sumν) + else + j += 1 + MA.add_mul!!(f, convert(T, set.α[i, var]), ν[j]) + end + end + return MOI.add_constraint(model, f, MOI.EqualTo(zero(T))) + end + scalars = MOI.Utilities.eachscalar(func) + f = MOI.Utilities.operate( + vcat, + T, + scalars[set.k] + sumν, + scalars[setdiff(1:m, set.k)], + MOI.VectorOfVariables(ν), + ) + relative_entropy_constraint = + MOI.add_constraint(model, f, MOI.RelativeEntropyCone(2m - 1)) + return AGEBridge{T,F,G,H}(set.k, ceq, relative_entropy_constraint) +end + +function MOI.supports_constraint( + ::Type{<:AGEBridge{T}}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialAGECone}, +) where {T} + return true +end + +function MOI.Bridges.added_constrained_variable_types(::Type{<:AGEBridge}) + return Tuple{Type}[] +end + +function MOI.Bridges.added_constraint_types( + ::Type{<:AGEBridge{T,F,G}}, +) where {T,F,G} + return [(F, MOI.EqualTo{T}), (G, MOI.RelativeEntropyCone)] +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{<:AGEBridge{T}}, + H::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialAGECone}, +) where {T} + S = MOI.Utilities.scalar_type(H) + F = MOI.Utilities.promote_operation( + +, + T, + S, + MOI.ScalarAffineFunction{Float64}, + ) + G = MOI.Utilities.promote_operation(vcat, T, S, F) + return AGEBridge{T,F,G,H} +end diff --git a/src/RelativeEntropy/bridges/sage.jl b/src/RelativeEntropy/bridges/sage.jl new file mode 100644 index 0000000..6a89996 --- /dev/null +++ b/src/RelativeEntropy/bridges/sage.jl @@ -0,0 +1,90 @@ +struct SAGEBridge{T,F,G} <: MOI.Bridges.Constraint.AbstractBridge + ν::Matrix{MOI.VariableIndex} + age_constraints::Vector{ + MOI.ConstraintIndex{MOI.VectorOfVariables,SignomialAGECone}, + } + equality_constraints::Vector{MOI.ConstraintIndex{F,MOI.EqualTo{T}}} +end + +function MOI.Bridges.Constraint.bridge_constraint( + ::Type{SAGEBridge{T,F,G}}, + model, + func::G, + set::SignomialSAGECone, +) where {T,F,G} + m = size(set.α, 1) + ν = Matrix{MOI.VariableIndex}(undef, m, m) + A = SignomialAGECone + age_constraints = + Vector{MOI.ConstraintIndex{MOI.VectorOfVariables,A}}(undef, m) + for k in 1:m + ν[k, :], age_constraints[k] = + MOI.add_constrained_variables(model, A(set.α, k)) + end + scalars = MOI.Utilities.eachscalar(func) + equality_constraints = map(1:m) do i + f = MOI.ScalarAffineFunction( + MOI.ScalarAffineTerm.(one(T), ν[:, i]), + zero(T), + ) + return MOI.Utilities.normalize_and_add_constraint( + model, + MA.sub!!(f, scalars[i]), + MOI.EqualTo(zero(T)), + ) + end + return SAGEBridge{T,F,G}(ν, age_constraints, equality_constraints) +end + +function MOI.supports_constraint( + ::Type{<:SAGEBridge{T}}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialSAGECone}, +) where {T} + return true +end + +function MOI.Bridges.added_constrained_variable_types(::Type{<:SAGEBridge}) + return Tuple{Type}[(SignomialAGECone,)] +end + +function MOI.Bridges.added_constraint_types( + ::Type{<:SAGEBridge{T,F}}, +) where {T,F} + return Tuple{Type,Type}[(F, MOI.EqualTo{T})] +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{<:SAGEBridge{T}}, + G::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:SignomialSAGECone}, +) where {T} + S = MOI.Utilities.scalar_type(G) + F = MOI.Utilities.promote_operation( + -, + T, + MOI.ScalarAffineFunction{Float64}, + S, + ) + return SAGEBridge{T,F,G} +end + +function MOI.get( + model::MOI.ModelLike, + attr::DecompositionAttribute, + bridge::SAGEBridge, +) + return filter!( + !isempty, + [ + filter( + x -> !isapprox(x, zero(x), atol = attr.tol), + MOI.get( + model, + MOI.VariablePrimal(attr.result_index), + bridge.ν[k, :], + ), + ) for k in axes(bridge.ν, 1) + ], + ) +end diff --git a/src/RelativeEntropy/bridges/signomial.jl b/src/RelativeEntropy/bridges/signomial.jl new file mode 100644 index 0000000..e7bb45a --- /dev/null +++ b/src/RelativeEntropy/bridges/signomial.jl @@ -0,0 +1,85 @@ +""" + SignomialBridge{T,S,P,F} <: MOI.Bridges.Constraint.AbstractBridge + +We use the Signomial Representative `SR` equation of [MCW21]. + +[MCW20] Riley Murray, Venkat Chandrasekaran, Adam Wierman +"Newton Polytopes and Relative Entropy Optimization" +https://arxiv.org/abs/1810.01614 +[MCW21] Murray, Riley, Venkat Chandrasekaran, and Adam Wierman. +"Signomial and polynomial optimization via relative entropy and partial dualization." +Mathematical Programming Computation 13 (2021): 257-295. +https://arxiv.org/pdf/1907.00814.pdf +""" +struct SignomialBridge{T,S,P,F} <: MOI.Bridges.Constraint.AbstractBridge + constraint::MOI.ConstraintIndex{F,S} +end + +_signomial(set::PolynomialSAGECone) = SignomialSAGECone(set.α) +_signomial(::Type{PolynomialSAGECone}) = SignomialSAGECone +_signomial(set::PolynomialAGECone) = SignomialAGECone(set.α, set.k) +_signomial(::Type{PolynomialAGECone}) = SignomialAGECone + +function MOI.Bridges.Constraint.bridge_constraint( + ::Type{SignomialBridge{T,S,P,F}}, + model, + func::F, + set, +) where {T,S,P,F} + g = MOI.Utilities.scalarize(func) + for i in eachindex(g) + if any(isodd, set.α[i, :]) + vi = MOI.add_variable(model) + # vi ≤ -|g[i]| + MOI.Utilities.normalize_and_add_constraint( + model, + one(T) * vi - g[i], + MOI.LessThan(zero(T)), + ) + MOI.Utilities.normalize_and_add_constraint( + model, + one(T) * vi + g[i], + MOI.LessThan(zero(T)), + ) + g[i] = vi + end + end + constraint = + MOI.add_constraint(model, MOI.Utilities.vectorize(g), _signomial(set)) + return SignomialBridge{T,S,P,F}(constraint) +end + +function MOI.supports_constraint( + ::Type{<:SignomialBridge{T}}, + ::Type{<:MOI.AbstractVectorFunction}, + ::Type{<:Union{PolynomialSAGECone,PolynomialAGECone}}, +) where {T} + return true +end + +function MOI.Bridges.added_constrained_variable_types(::Type{<:SignomialBridge}) + return Tuple{Type}[(MOI.Reals,)] +end + +function MOI.Bridges.added_constraint_types( + ::Type{<:SignomialBridge{T,S,P,F}}, +) where {T,S,P,F} + return [(F, S)] +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{<:SignomialBridge{T}}, + F::Type{<:MOI.AbstractVectorFunction}, + P::Type{<:Union{PolynomialSAGECone,PolynomialAGECone}}, +) where {T} + S = _signomial(P) + return SignomialBridge{T,S,P,F} +end + +function MOI.get( + model::MOI.ModelLike, + attr::DecompositionAttribute, + bridge::SignomialBridge, +) + return MOI.get(model, attr, bridge.constraint) +end diff --git a/test/Project.toml b/test/Project.toml index 2b094bf..125b4cf 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -1,5 +1,6 @@ [deps] DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +ECOS = "e2685f51-7e38-5353-a97d-a921fd2c8199" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/test/relative_entropy.jl b/test/relative_entropy.jl new file mode 100644 index 0000000..c32a946 --- /dev/null +++ b/test/relative_entropy.jl @@ -0,0 +1,106 @@ +module TestRelativeEntropy + +using Test + +import MultivariatePolynomials as MP +using SemialgebraicSets + +import DynamicPolynomials +import ECOS + +using JuMP +using PolyJuMP + +function _test_motzkin(x, y, T, solver, set, feasible, square, neg) + model = Model(solver) + a = square ? x^2 : x + b = square ? y^2 : y + PolyJuMP.setpolymodule!(model, PolyJuMP.RelativeEntropy) + if neg + motzkin = -a^2 * b - a * b^2 + one(T) - 3a * b + else + motzkin = a^2 * b + a * b^2 + one(T) - 3a * b + end + con_ref = @constraint(model, motzkin in set) + optimize!(model) + if feasible + @test termination_status(model) == MOI.OPTIMAL + @test primal_status(model) == MOI.FEASIBLE_POINT + if set isa Union{ + PolyJuMP.RelativeEntropy.SignomialSAGESet, + PolyJuMP.RelativeEntropy.PolynomialSAGESet, + } + d = PolyJuMP.RelativeEntropy.decomposition(con_ref; tol = 1e-6) + p = MP.polynomial(d) + if set isa PolyJuMP.RelativeEntropy.SignomialSAGESet + @test p ≈ motzkin atol = 1e-6 + else + for m in MP.monomials(p - motzkin) + @test MP.coefficient(p, m) ≈ MP.coefficient(motzkin, m) atol = + 1e-6 + end + end + end + else + @test termination_status(model) == MOI.INFEASIBLE + end +end + +function test_motzkin(x, y, T, solver) + set = PolyJuMP.RelativeEntropy.SignomialAGESet(x^2 * y^2) + _test_motzkin(x, y, T, solver, set, true, true, false) + set = PolyJuMP.RelativeEntropy.SignomialAGESet(x * y) + _test_motzkin(x, y, T, solver, set, true, false, false) + set = PolyJuMP.RelativeEntropy.SignomialAGESet(x^4 * y^2) + _test_motzkin(x, y, T, solver, set, false, true, false) + set = PolyJuMP.RelativeEntropy.SignomialAGESet(x^2 * y) + _test_motzkin(x, y, T, solver, set, false, false, false) + set = PolyJuMP.RelativeEntropy.SignomialSAGESet() + _test_motzkin(x, y, T, solver, set, true, true, false) + _test_motzkin(x, y, T, solver, set, true, false, false) + _test_motzkin(x, y, T, solver, set, false, true, true) + _test_motzkin(x, y, T, solver, set, false, false, true) + set = PolyJuMP.RelativeEntropy.PolynomialAGESet(x^2 * y^2) + _test_motzkin(x, y, T, solver, set, true, true, false) + set = PolyJuMP.RelativeEntropy.PolynomialAGESet(x * y) + _test_motzkin(x, y, T, solver, set, false, false, false) + set = PolyJuMP.RelativeEntropy.PolynomialAGESet(x^4 * y^2) + _test_motzkin(x, y, T, solver, set, false, true, false) + set = PolyJuMP.RelativeEntropy.PolynomialAGESet(x^2 * y) + _test_motzkin(x, y, T, solver, set, false, false, false) + set = PolyJuMP.RelativeEntropy.PolynomialSAGESet() + _test_motzkin(x, y, T, solver, set, true, true, false) + set = PolyJuMP.RelativeEntropy.PolynomialSAGESet() + _test_motzkin(x, y, T, solver, set, false, false, false) + return +end + +import ECOS +const SOLVERS = + [optimizer_with_attributes(ECOS.Optimizer, MOI.Silent() => true)] + +function runtests(x, y, T) + for name in names(@__MODULE__; all = true) + if startswith("$name", "test_") + @testset "$(name) $solver)" for solver in SOLVERS + getfield(@__MODULE__, name)(x, y, T, solver) + end + end + end +end + +end # module + +using Test + +import DynamicPolynomials +@testset "DynamicPolynomials" begin + DynamicPolynomials.@polyvar(x, y) + TestRelativeEntropy.runtests(x, y, Float64) +end + +import TypedPolynomials +@testset "DynamicPolynomials" begin + TypedPolynomials.@polyvar(x, y) + TestRelativeEntropy.runtests(x, y, Float64) +end diff --git a/test/runtests.jl b/test/runtests.jl index 5083959..77fe131 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,3 +6,4 @@ include("functions.jl") include("Mock/mock_tests.jl") include("kkt.jl") +include("relative_entropy.jl")