Add SAGE Optimizer

README.md

@@ -63,7 +63,7 @@ This allows mixing SAGE and SOS constraints in the same model.
 
 ### Polynomial optimization
 
-PolyJuMP also allows solving polynomial optimization problems using the `QCQP` and `KKT` solvers.
+PolyJuMP also allows solving polynomial optimization problems using the `QCQP`, `KKT` and `SAGE` solvers.
 Polynomial optimization problems do not involve any symbolic variables from DynamicPolynomials or TypedPolynomials,
 instead all variables are JuMP decision variables.
 
@@ -89,6 +89,53 @@ model = Model(optimizer_with_attributes(
 ))
 ```
 
+The `SAGE` solver computes a bound on the optimal objective value using the
+SAGE relaxation of the problem [CP16, MCW21]: it certifies the nonnegativity
+of the Lagrangian with the SAGE cone, using one SAGE multiplier for each
+inequality constraint and one free polynomial multiplier for each equality
+constraint.
+It is parametrized by an inner solver for the resulting relative entropy
+program. For instance, to compute a lower bound on the minimum of the Motzkin
+polynomial with `ECOS.Optimizer` as inner solver, use:
+```julia
+using JuMP, PolyJuMP, ECOS
+model = Model(() -> PolyJuMP.SAGE.Optimizer(ECOS.Optimizer))
+@variable(model, x)
+@variable(model, y)
+@objective(model, Min, x^4 * y^2 + x^2 * y^4 + 1 - 3 * x^2 * y^2)
+optimize!(model)
+objective_bound(model) # ≈ 0
+```
+The bound is returned as `MOI.ObjectiveBound`. In addition, candidate
+solutions are recovered from the dual of the SAGE constraint, which is a
+vector of pseudo-moments, following [MCW21, Section 4.2] (see also its
+reference implementation `poly_solrec` in
+[sageopt](https://github.com/rileyjmurray/sageopt)); `result_count(model)`
+gives the number of candidates found, sorted by feasibility and objective
+value. In the example above, the four minimizers `(±1, ±1)` are recovered:
+```julia
+value(x; result = 1), value(y; result = 1) # ≈ (1, 1)
+```
+The maximum degree of the multiplier of a constraint is chosen with the
+`PolyJuMP.MultiplierMaxdegree` constraint attribute:
+```julia
+@constraint(model, con, x^2 >= 1)
+MOI.set(model, PolyJuMP.MultiplierMaxdegree(), con, 2)
+```
+The `SumOfSquares.Optimizer` of [SumOfSquares.jl](https://github.com/jump-dev/SumOfSquares.jl)
+is the analogous solver certifying the nonnegativity of the Lagrangian with
+the SOS cone instead; increasing the `PolyJuMP.MultiplierMaxdegree` attributes
+then gives the higher levels of the Lasserre hierarchy.
+
+[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.
+*Signomials and polynomial optimization via relative entropy and partial
+dualization.* Mathematical Programming Computation 13 (2021): 257-295.
+[arXiv:1907.00814](https://arxiv.org/abs/1907.00814)
+
 ## Documentation
 
 Documentation for `PolyJuMP.jl` is included in the

src/KKT/KKT.jl

@@ -13,7 +13,7 @@ Base.@kwdef mutable struct Options{T}
     feasibility_tolerance::T = Base.rtoldefault(T)
 end
 
-mutable struct Optimizer{T} <: MOI.AbstractOptimizer
+mutable struct Optimizer{T} <: PolyJuMP.AbstractPolynomialOptimizer{T}
     model::PolyJuMP.Model{T}
     options::Options{T}
     # Result
@@ -40,34 +40,6 @@ Optimizer() = Optimizer{Float64}()
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "PolyJuMP.KKT"
 
-MOI.is_empty(model::Optimizer) = MOI.is_empty(model.model)
-
-function MOI.empty!(model::Optimizer)
-    MOI.empty!(model.model)
-    invalidate_solutions!(model)
-    return
-end
-
-function MOI.supports(model::Optimizer, attr::MOI.AbstractModelAttribute)
-    return MOI.supports(model.model, attr)
-end
-
-function MOI.set(model::Optimizer, attr::MOI.AbstractModelAttribute, value)
-    MOI.set(model.model, attr, value)
-    invalidate_solutions!(model)
-    return
-end
-
-function MOI.get(
-    model::Optimizer,
-    attr::Union{
-        MOI.AbstractModelAttribute,
-        MOI.Bridges.ListOfNonstandardBridges,
-    },
-)
-    return MOI.get(model.model, attr)
-end
-
 function MOI.supports(::Optimizer{T}, attr::MOI.RawOptimizerAttribute) where {T}
     return hasfield(Options{T}, Symbol(attr.name))
 end
@@ -87,44 +59,14 @@ function MOI.get(model::Optimizer, attr::MOI.RawOptimizerAttribute)
     return getfield(model.options, Symbol(attr.name))
 end
 
-function invalidate_solutions!(model::Optimizer)
+function PolyJuMP._invalidate!(model::Optimizer)
     empty!(model.solutions)
     model.solve_time = NaN
     model.termination_status = MOI.OPTIMIZE_NOT_CALLED
     model.raw_status = ""
     return
 end
 
-MOI.is_valid(model::Optimizer, i::MOI.Index) = MOI.is_valid(model.model, i)
-
-function MOI.add_variable(model::Optimizer)
-    invalidate_solutions!(model)
-    return MOI.add_variable(model.model)
-end
-
-function MOI.supports_constraint(
-    model::Optimizer,
-    ::Type{F},
-    ::Type{S},
-) where {F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
-    return MOI.supports_constraint(model.model, F, S)
-end
-
-function MOI.add_constraint(
-    model::Optimizer,
-    func::MOI.AbstractFunction,
-    set::MOI.AbstractSet,
-)
-    ci = MOI.add_constraint(model.model, func, set)
-    invalidate_solutions!(model)
-    return ci
-end
-
-MOI.supports_incremental_interface(::Optimizer) = true
-function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike)
-    return MOI.Utilities.default_copy_to(dest, src)
-end
-
 function _add_to_system(system, lagrangian, ::SS.FullSpace, ::Bool)
     return lagrangian
 end
@@ -176,7 +118,7 @@ function _square(x::Vector{T}, n) where {T}
     return T[(i + n in eachindex(x)) ? x[i] : x[i]^2 for i in eachindex(x)]
 end
 
-function _optimize!(model::Optimizer{T}) where {T}
+function PolyJuMP._optimize!(model::Optimizer{T}) where {T}
     if isnothing(model.options.solver)
         system = SS.AlgebraicSet{T,PolyJuMP.PolyType{T}}()
     else
@@ -242,12 +184,6 @@ function _optimize!(model::Optimizer{T}) where {T}
     return
 end
 
-function MOI.optimize!(model::Optimizer)
-    return model.solve_time = @elapsed _optimize!(model)
-end
-
-MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time
-
 function MOI.get(model::Optimizer, ::MOI.RawStatusString)
     if model.termination_status === MOI.OPTIMIZE_NOT_CALLED
         return "`optimize!` has not yet been called"

src/PolyJuMP.jl

@@ -33,6 +33,7 @@ include("data.jl")
 include("default.jl")
 
 include("model.jl")
+include("optimizer.jl")
 include("KKT/KKT.jl")
 include("QCQP/QCQP.jl")
 include("SAGE/SAGE.jl")

src/SAGE/SAGE.jl

@@ -182,4 +182,6 @@ function PolyJuMP.bridges(
     )]
 end
 
+include("optimizer.jl")
+
 end

src/SAGE/bridges/age.jl

@@ -103,6 +103,32 @@ function MOI.Bridges.added_constraint_types(
     return [(F, MOI.EqualTo{T}), (G, MOI.RelativeEntropyCone)]
 end
 
+# The coefficients `c` of the AGE constraint only appear in the relative
+# entropy constraint `(c_k + ∑ν, c_{-k}, ν) ∈ RelativeEntropyCone` where `k`
+# is the index of the distinguished monomial. The dual is the adjoint of this
+# map applied to the dual `(u, v, w)` of the relative entropy constraint,
+# that is, `u` for the entry `k` and `v` for the other entries. This is the
+# dual AGE cone characterization given by the conic duality of the relative
+# entropy formulation [MCW21, (2)].
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::AGEBridge,
+)
+    dual = MOI.get(model, attr, bridge.relative_entropy_constraint)
+    m = div(length(dual) + 1, 2)
+    v = Vector{eltype(dual)}(undef, m)
+    v[bridge.k] = dual[1]
+    j = 1
+    for i in 1:m
+        if i != bridge.k
+            j += 1
+            v[i] = dual[j]
+        end
+    end
+    return v
+end
+
 function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{<:AGEBridge{T}},
     H::Type{<:MOI.AbstractVectorFunction},

src/SAGE/bridges/sage.jl

@@ -69,6 +69,21 @@ function MOI.Bridges.Constraint.concrete_bridge_type(
     return SAGEBridge{T,F,G}
 end
 
+# The signomial SAGE constraint `func ∈ SAGE` is reformulated into the
+# equality constraints `∑_k ν[k, i] - func_i = 0` in which `func` appears
+# with coefficient `-1` so the dual is `-μ` where `μ` is the dual of these
+# equality constraints; this is the adjoint of the reformulation map.
+# Since the SAGE cone is the sum of the AGE cones, its dual is the
+# intersection of the duals of the AGE cones and indeed, at the optimum, the
+# dual of each constraint `ν[k, :] ∈ AGE` also equals `-μ`.
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::SAGEBridge,
+)
+    return [-MOI.get(model, attr, ci) for ci in bridge.equality_constraints]
+end
+
 function MOI.get(
     model::MOI.ModelLike,
     attr::DecompositionAttribute,

src/SAGE/bridges/signomial.jl

@@ -11,30 +11,46 @@ https://arxiv.org/abs/1810.01614
 Mathematical Programming Computation 13 (2021): 257-295.
 https://arxiv.org/pdf/1907.00814.pdf
 """
-struct SignomialsBridge{T,S,P,F} <: MOI.Bridges.Constraint.AbstractBridge
+struct SignomialsBridge{T,S,P,F,G} <: MOI.Bridges.Constraint.AbstractBridge
+    # Indices of the rows `i` of `set.α` that have an odd entry
+    odd::Vector{Int}
+    # For each odd row `i`, the constraint `vi - g[i] ≤ 0`
+    lower::Vector{MOI.ConstraintIndex{G,MOI.LessThan{T}}}
+    # For each odd row `i`, the constraint `vi + g[i] ≤ 0`
+    upper::Vector{MOI.ConstraintIndex{G,MOI.LessThan{T}}}
     constraint::MOI.ConstraintIndex{F,S}
 end
 
 function MOI.Bridges.Constraint.bridge_constraint(
-    ::Type{SignomialsBridge{T,S,P,F}},
+    ::Type{SignomialsBridge{T,S,P,F,G}},
     model,
     func::F,
     set,
-) where {T,S,P,F}
+) where {T,S,P,F,G}
     g = MOI.Utilities.scalarize(func)
+    odd = Int[]
+    lower = MOI.ConstraintIndex{G,MOI.LessThan{T}}[]
+    upper = MOI.ConstraintIndex{G,MOI.LessThan{T}}[]
     for i in eachindex(g)
         if any(isodd, set.α[i, :])
             vi = MOI.add_variable(model)
+            push!(odd, i)
             # vi ≤ -|g[i]|
-            MOI.Utilities.normalize_and_add_constraint(
-                model,
-                one(T) * vi - g[i],
-                MOI.LessThan(zero(T)),
+            push!(
+                lower,
+                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)),
+            push!(
+                upper,
+                MOI.Utilities.normalize_and_add_constraint(
+                    model,
+                    one(T) * vi + g[i],
+                    MOI.LessThan(zero(T)),
+                ),
             )
             g[i] = vi
         end
@@ -44,7 +60,7 @@ function MOI.Bridges.Constraint.bridge_constraint(
         MOI.Utilities.vectorize(g),
         Cone(Signomials(set.cone.monomial), set.α),
     )
-    return SignomialsBridge{T,S,P,F}(constraint)
+    return SignomialsBridge{T,S,P,F,G}(odd, lower, upper, constraint)
 end
 
 function MOI.supports_constraint(
@@ -62,17 +78,23 @@ function MOI.Bridges.added_constrained_variable_types(
 end
 
 function MOI.Bridges.added_constraint_types(
-    ::Type{<:SignomialsBridge{T,S,P,F}},
-) where {T,S,P,F}
-    return [(F, S)]
+    ::Type{<:SignomialsBridge{T,S,P,F,G}},
+) where {T,S,P,F,G}
+    return [(F, S), (G, MOI.LessThan{T})]
 end
 
 function MOI.Bridges.Constraint.concrete_bridge_type(
     ::Type{<:SignomialsBridge{T}},
     F::Type{<:MOI.AbstractVectorFunction},
     P::Type{Cone{Polynomials{M}}},
 ) where {T,M}
-    return SignomialsBridge{T,Cone{Signomials{M}},P,F}
+    G = MOI.Utilities.promote_operation(
+        -,
+        T,
+        MOI.ScalarAffineFunction{T},
+        MOI.Utilities.scalar_type(F),
+    )
+    return SignomialsBridge{T,Cone{Signomials{M}},P,F,G}
 end
 
 function MOI.get(
@@ -82,3 +104,26 @@ function MOI.get(
 )
     return MOI.get(model, attr, bridge.constraint)
 end
+
+# The dual of the polynomial SAGE constraint is the adjoint of the linear map
+# used in the reformulation, applied to the duals of the constraints created
+# by the bridge. For a row `i` with even exponents, `g[i]` only appears in row
+# `i` of the signomial constraint so the dual is the corresponding entry `w[i]`
+# of its dual `w`. For an odd row, `g[i]` appears with coefficient `-1` in
+# `lower[i]` and `+1` in `upper[i]` so the dual is the difference of the duals
+# of these two constraints. The result `v` satisfies `|v[i]| ≤ w[i]` for odd
+# rows, which matches the characterization of the dual of the polynomial SAGE
+# cone in terms of the dual of the signomial SAGE cone of [MCW20].
+function MOI.get(
+    model::MOI.ModelLike,
+    attr::MOI.ConstraintDual,
+    bridge::SignomialsBridge,
+)
+    v = MOI.get(model, attr, bridge.constraint)
+    for (j, i) in enumerate(bridge.odd)
+        v[i] =
+            MOI.get(model, attr, bridge.upper[j]) -
+            MOI.get(model, attr, bridge.lower[j])
+    end
+    return v
+end