From 57da88c261b3a821e7bff1f7c0e4fa58927f064e Mon Sep 17 00:00:00 2001 From: Giovanni Karra Date: Tue, 20 Oct 2026 13:38:34 +0200 Subject: [PATCH] Add polynomial preprocessing bridges and QCQP.optimizer Add `ToReformulatedPolynomialBridge` (constraint and objective), which substitutes high-degree subexpressions of a `MOI.ScalarNonlinearFunction` with bounded auxiliary variables before conversion to a polynomial, and `QCQP.optimizer(model; preprocessing)` stacking these bridges on top of `QCQP.Optimizer`. `ToPolynomialBridge` now only bridges `MOI.ScalarNonlinearFunction` and non-integer exponents now throw instead of being rounded. --- src/Bridges/Bridges.jl | 4 + src/Bridges/Constraint/Constraint.jl | 1 + .../Constraint/to_reformulated_polynomial.jl | 79 +++++ src/Bridges/Objective/Objective.jl | 4 +- .../Objective/to_reformulated_polynomial.jl | 58 ++++ src/Bridges/Utilities.jl | 287 ++++++++++++++++++ src/QCQP/MOI_wrapper.jl | 16 + src/nl_to_polynomial.jl | 5 +- test/qcqp.jl | 33 ++ 9 files changed, 483 insertions(+), 4 deletions(-) create mode 100644 src/Bridges/Constraint/to_reformulated_polynomial.jl create mode 100644 src/Bridges/Objective/to_reformulated_polynomial.jl create mode 100644 src/Bridges/Utilities.jl diff --git a/src/Bridges/Bridges.jl b/src/Bridges/Bridges.jl index 656bfda..d9e3ac5 100644 --- a/src/Bridges/Bridges.jl +++ b/src/Bridges/Bridges.jl @@ -1,5 +1,9 @@ module Bridges +import MathOptInterface as MOI +using PolyJuMP + +include("Utilities.jl") # Constraint bridges include("Constraint/Constraint.jl") # Objective bridges diff --git a/src/Bridges/Constraint/Constraint.jl b/src/Bridges/Constraint/Constraint.jl index 87a0707..7889cb7 100644 --- a/src/Bridges/Constraint/Constraint.jl +++ b/src/Bridges/Constraint/Constraint.jl @@ -13,5 +13,6 @@ include("zero_polynomial.jl") include("zero_polynomial_in_algebraic_set.jl") include("plus_minus.jl") include("to_polynomial.jl") +include("to_reformulated_polynomial.jl") end # module diff --git a/src/Bridges/Constraint/to_reformulated_polynomial.jl b/src/Bridges/Constraint/to_reformulated_polynomial.jl new file mode 100644 index 0000000..eddd7d9 --- /dev/null +++ b/src/Bridges/Constraint/to_reformulated_polynomial.jl @@ -0,0 +1,79 @@ +""" + ToReformulatedPolynomialBridge{T,S} <: Bridges.Constraint.AbstractBridge + +`ToReformulatedPolynomialBridge` implements a sparsity-preserving +reformulation of `f(x) in S` into `g(x, y) in S` where `f` is a +`MOI.ScalarNonlinearFunction`: products and powers of sums are replaced by +auxiliary variables `y`, so that converting `g` into a polynomial does not +expand them. + +Since source and target are both `MOI.ScalarNonlinearFunction`, this bridge +doesn't have any effect when added to a `MOI.Bridges.LazyBridgeOptimizer`; use it in a +`MOI.Bridges.Constraint.SingleBridgeOptimizer` on top of an optimizer converting +`MOI.ScalarNonlinearFunction` into polynomials, e.g., [`PolyJuMP.QCQP.Optimizer`](@ref). + +## Source node + +`ToReformulatedPolynomialBridge` supports: + + * `MOI.ScalarNonlinearFunction` in `S` + +## Target nodes + +`ToReformulatedPolynomialBridge` creates: + + * `MOI.ScalarNonlinearFunction` in `S` + * `MOI.ScalarNonlinearFunction` in `MOI.EqualTo{T}`, defining the + auxiliary variables + * `MOI.VariableIndex` in `MOI.Interval{T}`, the auxiliary variables +""" +struct ToReformulatedPolynomialBridge{T,S} <: + MOI.Bridges.Constraint.AbstractFunctionConversionBridge{MOI.ScalarNonlinearFunction,S} + constraint::MOI.ConstraintIndex{MOI.ScalarNonlinearFunction,S} +end + +function MOI.Bridges.Constraint.bridge_constraint( + ::Type{ToReformulatedPolynomialBridge{T,S}}, + model::MOI.ModelLike, + func::MOI.ScalarNonlinearFunction, + set::S, +) where {T,S} + reformulated = PolyJuMP.Bridges._to_reformulated_polynomial!(model, T, func) + constraint = MOI.add_constraint(model, reformulated, set) + return ToReformulatedPolynomialBridge{T,S}(constraint) +end + +function MOI.Bridges.Objective.supports_objective_function( + ::Type{ToReformulatedPolynomialBridge{T}}, + ::Type{F}, +) where {T,F<:MOI.ScalarNonlinearFunction} + return MOI.Utilities.is_coefficient_type(F, T) +end + + +function MOI.supports_constraint( + ::Type{ToReformulatedPolynomialBridge{T}}, + ::Type{MOI.ScalarNonlinearFunction}, + ::Type{<:MOI.AbstractScalarSet}, +) where {T} + return true +end + +function MOI.Bridges.Constraint.concrete_bridge_type( + ::Type{<:ToReformulatedPolynomialBridge{T}}, + ::Type{MOI.ScalarNonlinearFunction}, + ::Type{S}, +) where {T,S<:MOI.AbstractScalarSet} + return ToReformulatedPolynomialBridge{T,S} +end + + +function MOI.get( + model::MOI.ModelLike, + attr::MOI.Bridges.ObjectiveFunctionValue, + ::ToReformulatedPolynomialBridge{T}, +) where {T} + F = FuncType{T} + attr_f = MOI.Bridges.ObjectiveFunctionValue{F}(attr.result_index) + return MOI.get(model, attr_f) +end diff --git a/src/Bridges/Objective/Objective.jl b/src/Bridges/Objective/Objective.jl index 34d2f44..2d0d25c 100644 --- a/src/Bridges/Objective/Objective.jl +++ b/src/Bridges/Objective/Objective.jl @@ -1,10 +1,10 @@ module Objective -import MathOptInterface -const MOI = MathOptInterface +import MathOptInterface as MOI using PolyJuMP include("to_polynomial.jl") +include("to_reformulated_polynomial.jl") end # module diff --git a/src/Bridges/Objective/to_reformulated_polynomial.jl b/src/Bridges/Objective/to_reformulated_polynomial.jl new file mode 100644 index 0000000..83a8c69 --- /dev/null +++ b/src/Bridges/Objective/to_reformulated_polynomial.jl @@ -0,0 +1,58 @@ +""" + ToReformulatedPolynomialBridge{T} <: Bridges.Objective.AbstractBridge + +`ToReformulatedPolynomialBridge` implements a sparsity-preserving +reformulation of the objective `f(x)` into `g(x, y)` where `f` is a +`MOI.ScalarNonlinearFunction`: products and powers of sums are replaced by +auxiliary variables `y`, so that converting `g` into a polynomial does not +expand them. + +Since source and target are both `MOI.ScalarNonlinearFunction`, this bridge +doesn't have any effect when added to a `MOI.Bridges.LazyBridgeOptimizer`; use it in a +`MOI.Bridges.Objective.SingleBridgeOptimizer` on top of an optimizer converting +`MOI.ScalarNonlinearFunction` into polynomials, e.g., [`PolyJuMP.QCQP.Optimizer`](@ref). + +## Source node + +`ToReformulatedPolynomialBridge` supports: + + * `MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}` + +## Target nodes + +`ToReformulatedPolynomialBridge` creates: + + * `MOI.ObjectiveFunction{MOI.ScalarNonlinearFunction}` + * `MOI.ScalarNonlinearFunction` in `MOI.EqualTo{T}`, defining the auxiliary + variables + * `MOI.VariableIndex` in `MOI.Interval{T}`, the auxiliary variables +""" +struct ToReformulatedPolynomialBridge{T} <: MOI.Bridges.Objective.AbstractBridge end + +function MOI.Bridges.Objective.bridge_objective( + ::Type{ToReformulatedPolynomialBridge{T}}, + model::MOI.ModelLike, + func::MOI.ScalarNonlinearFunction, +) where {T} + reformulated = PolyJuMP.Bridges._to_reformulated_polynomial!(model, T, func) + F = MOI.ScalarNonlinearFunction + MOI.set(model, MOI.ObjectiveFunction{F}(), reformulated) + return ToReformulatedPolynomialBridge{T}() +end + +function MOI.Bridges.Objective.supports_objective_function( + ::Type{ToReformulatedPolynomialBridge{T}}, + ::Type{F}, +) where {T,F<:MOI.ScalarNonlinearFunction} + return MOI.Utilities.is_coefficient_type(F, T) +end + +function MOI.get( + model::MOI.ModelLike, + attr::MOI.Bridges.ObjectiveFunctionValue, + ::ToReformulatedPolynomialBridge{T}, +) where {T} + F = FuncType{T} + attr_f = MOI.Bridges.ObjectiveFunctionValue{F}(attr.result_index) + return MOI.get(model, attr_f) +end diff --git a/src/Bridges/Utilities.jl b/src/Bridges/Utilities.jl new file mode 100644 index 0000000..977e1da --- /dev/null +++ b/src/Bridges/Utilities.jl @@ -0,0 +1,287 @@ +import MultivariatePolynomials as MP +import DynamicPolynomials +using IntervalArithmetic +const IA = IntervalArithmetic + +const VariableOrder = + DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder} +const MonomialOrder = MP.Graded{MP.LexOrder} +const VarType = DynamicPolynomials.Variable{VariableOrder,MonomialOrder} +const PolyType{T} = DynamicPolynomials.Polynomial{VariableOrder,MonomialOrder,T} +const FuncType{T} = PolyJuMP.ScalarPolynomialFunction{T,PolyType{T}} + +# Auxiliary variables created by the sparsity-preserving reformulation, keyed +# by `_expr_key` of the expression they replace. +const ExprMap = Dict{Any,MOI.VariableIndex} + +_expr_key(x::MOI.VariableIndex) = x +_expr_key(x::Number) = x +_expr_key(x) = _expr_key(convert(MOI.ScalarNonlinearFunction, x)) +function _expr_key(f::MOI.ScalarNonlinearFunction) + args = Any[_expr_key(arg) for arg in f.args] + if f.head in (:+, :*) + sort!(args; by = string) + end + return (f.head, args...) +end + + +""" + _to_reformulated_polynomial!(model, T, expr) + +Returns `expr` where products and powers of sums whose expansion would be +large are replaced by auxiliary variables : + +* in a product of degree larger than 2 with several factors containing a + sum, these factors are substituted from the highest degree to the lowest + until only one remains, since expanding the product of one sum with + monomials does not increase the density much; +* the base of a power is substituted unless the power has degree less than + 3, or degree 3 with at most 2 variables. + +The auxiliary variables are added to `model` with bounds computed by +interval arithmetic and their defining equality constraints. +""" +function _to_reformulated_polynomial!( + model::MOI.ModelLike, + ::Type{T}, + expr, +) where {T} + d = ExprMap() + return _to_reformulated_polynomial!(model, T, d, expr) +end + +function _to_reformulated_polynomial!( + model::MOI.ModelLike, + ::Type{T}, + d::ExprMap, + expr, +) where {T} + return _to_reformulated_polynomial!(model, T, d, convert(MOI.ScalarNonlinearFunction, expr)) +end + +function _to_reformulated_polynomial!( + ::MOI.ModelLike, + ::Type{T}, + ::ExprMap, + expr::Union{MOI.VariableIndex,Number}, +) where {T} + return expr +end + +function _to_reformulated_polynomial!( + model::MOI.ModelLike, + ::Type{T}, + d::ExprMap, + expr::MOI.ScalarNonlinearFunction, +) where {T} + operator = expr.head + operands = expr.args + if operator == :* + remaining = count(_is_expr, operands) + if _degree(expr) > 2 && remaining > 1 + new_operands = Vector{Any}(undef, length(operands)) + for i in sortperm(operands; by = _degree, rev = true) + op = operands[i] + if remaining > 1 && _is_expr(op) + new_operands[i] = _substitute_expr(model, T, d, op) + remaining -= 1 + else + new_operands[i] = _to_reformulated_polynomial!(model, T, d, op) + end + end + return MOI.ScalarNonlinearFunction(:*, new_operands) + end + elseif operator == :^ && length(operands) == 2 + base, exponent = operands + deg = _degree(expr) + if deg < 3 || (deg == 3 && length(_var_index_set(expr)) <= 2) + base = _to_reformulated_polynomial!(model, T, d, base) + else + base = _substitute_expr(model, T, d, base) + end + return MOI.ScalarNonlinearFunction(:^, Any[base, exponent]) + end + return MOI.ScalarNonlinearFunction( + operator, + Any[_to_reformulated_polynomial!(model, T, d, op) for op in operands], + ) +end + +_degree(::Number) = 0 +_degree(::MOI.VariableIndex) = 1 +_degree(expr) = _degree(convert(MOI.ScalarNonlinearFunction, expr)) + +function _degree(expr::MOI.ScalarNonlinearFunction) + operator = expr.head + operands = expr.args + + if operator == :+ || operator == :- + return maximum(_degree, operands) + elseif operator == :* + return sum(_degree, operands) + elseif operator == :/ && length(operands) == 2 + return _degree(operands[1]) + elseif operator == :^ && length(operands) == 2 + if _degree(operands[2]) != 0 + error("Cannot infer the degree of $expr: the exponent is not constant.") + end + return _degree(operands[1]) * Int(_evaluate_const_expr(operands[2])) + end + return error("Unknown expression $expr, cannot infer degree") +end + +_evaluate_const_expr(expr::MOI.VariableIndex) = error("Cannot evaluate $expr, expr needs to be constant") +_evaluate_const_expr(expr::Number) = expr +_evaluate_const_expr(expr) = _evaluate_const_expr(convert(MOI.ScalarNonlinearFunction, expr)) + +function _evaluate_const_expr(expr::MOI.ScalarNonlinearFunction) + operator = expr.head + operands = expr.args + if operator == :+ + return sum(_evaluate_const_expr, operands) + elseif operator == :- && length(operands) == 2 + return _evaluate_const_expr(operands[1]) - _evaluate_const_expr(operands[2]) + elseif operator == :- && length(operands) == 1 + return -_evaluate_const_expr(operands[1]) + elseif operator == :* + return prod(_evaluate_const_expr, operands) + elseif operator == :/ && length(operands) == 2 + return _evaluate_const_expr(operands[1]) / _evaluate_const_expr(operands[2]) + elseif operator == :^ && length(operands) == 2 + return _evaluate_const_expr(operands[1])^_evaluate_const_expr(operands[2]) + end + return error("Cannot evaluate $expr") +end + +_var_index_set(::Number, ::Set{Int}) = nothing +_var_index_set(expr::MOI.VariableIndex, vars::Set{Int}) = push!(vars, expr.value) +function _var_index_set(expr, vars::Set{Int}) + return _var_index_set(convert(MOI.ScalarNonlinearFunction, expr), vars) +end + +function _var_index_set(expr)::Set{Int} + ret = Set{Int}() + return _var_index_set(convert(MOI.ScalarNonlinearFunction, expr), ret) +end + +function _var_index_set(expr::MOI.ScalarNonlinearFunction, vars::Set{Int}) + for operand in expr.args + _var_index_set(operand, vars) + end + return vars +end + +""" + _substitute_expr(model, T, d::ExprMap, expr) + +If `expr` is a sum or a difference, returns an auxiliary variable equal to +its reformulation (reusing the one created for an identical expression of +the same function); otherwise returns the reformulation of `expr`. +""" +function _substitute_expr( + model::MOI.ModelLike, + ::Type{T}, + d::ExprMap, + expr, +) where {T} + func = convert(MOI.ScalarNonlinearFunction, expr) + func = _to_reformulated_polynomial!(model, T, d, func) + if !(func.head in (:+, :-)) + return func + end + key = _expr_key(func) + v = get(d, key, nothing) + if v === nothing + lb, ub = IA.bounds(_get_nonlinfunc_bounds(model, T, func)) + v, _ = MOI.add_constrained_variable(model, MOI.Interval(lb, ub)) + MOI.Utilities.normalize_and_add_constraint( + model, + MOI.ScalarNonlinearFunction(:-, Any[func, v]), + MOI.EqualTo(zero(T)), + ) + d[key] = v + end + return v +end + +function _substitute_expr( + ::MOI.ModelLike, + ::Type{T}, + ::ExprMap, + expr::Union{MOI.VariableIndex,Number}, +) where {T} + return expr +end + +""" + _is_expr(expr) + +Returns whether `expr` contains a sum or a difference of several terms. +""" +function _is_expr(expr::MOI.ScalarNonlinearFunction)::Bool + return (expr.head in (:+, :-) && length(expr.args) > 1) || any(_is_expr, expr.args) +end +_is_expr(::Union{MOI.VariableIndex,Number}) = false +_is_expr(expr) = _is_expr(convert(MOI.ScalarNonlinearFunction, expr)) + +""" + _get_nonlinfunc_bounds(model, T, expr) + +Interval containing the values of `expr` given the bounds of the variables +in `model`, used as bounds of the auxiliary variable replacing `expr`. A +wrong (too narrow) interval cuts feasible solutions, so every operation must +return a valid enclosure. +""" +function _get_nonlinfunc_bounds(model::MOI.ModelLike, ::Type{T}, expr) where {T} + return _get_nonlinfunc_bounds( + model, + T, + convert(MOI.ScalarNonlinearFunction, expr), + ) +end + +function _get_nonlinfunc_bounds(::MOI.ModelLike, ::Type{T}, expr::Number) where {T} + return IA.interval(float(expr)) +end + +function _get_nonlinfunc_bounds( + model::MOI.ModelLike, + ::Type{T}, + expr::MOI.VariableIndex, +) where {T} + F = float(T) + lb, ub = MOI.Utilities.get_bounds(model, T, expr) + return IA.interval( + lb == typemin(lb) ? typemin(F) : float(lb), + ub == typemax(ub) ? typemax(F) : float(ub), + ) +end + +function _get_nonlinfunc_bounds( + model::MOI.ModelLike, + ::Type{T}, + expr::MOI.ScalarNonlinearFunction, +) where {T} + operator = expr.head + operands = expr.args + bounds(op) = _get_nonlinfunc_bounds(model, T, op) + if operator == :+ + return sum(bounds, operands) + elseif operator == :- && length(operands) == 1 + return -bounds(operands[1]) + elseif operator == :- && length(operands) == 2 + return bounds(operands[1]) - bounds(operands[2]) + elseif operator == :* + return prod(bounds, operands) + elseif operator == :/ && length(operands) == 2 + return bounds(operands[1]) / bounds(operands[2]) + elseif operator == :^ && length(operands) == 2 + base, exponent = operands + if _degree(exponent) == 0 && isinteger(_evaluate_const_expr(exponent)) + return bounds(base)^Int(_evaluate_const_expr(exponent)) + end + return bounds(base)^bounds(exponent) + end + return error("Unknown expression $expr, cannot compute its bounds") +end diff --git a/src/QCQP/MOI_wrapper.jl b/src/QCQP/MOI_wrapper.jl index 35b666d..54bdbbc 100644 --- a/src/QCQP/MOI_wrapper.jl +++ b/src/QCQP/MOI_wrapper.jl @@ -28,6 +28,22 @@ end Optimizer(model::MOI.ModelLike) = Optimizer{Float64}(model) +function optimizer(model::MOI.ModelLike; preprocessing=false) + optimizer = Optimizer(model) + if preprocessing + bridge = PolyJuMP.Bridges.Constraint.ToPolynomialBridge{Float64} + optimizer = MOI.Bridges.Constraint.SingleBridgeOptimizer{bridge}(optimizer) + bridge = PolyJuMP.Bridges.Objective.ToPolynomialBridge{Float64} + optimizer = MOI.Bridges.Objective.SingleBridgeOptimizer{bridge}(optimizer) + + bridge = PolyJuMP.Bridges.Constraint.ToReformulatedPolynomialBridge{Float64} + optimizer = MOI.Bridges.Constraint.SingleBridgeOptimizer{bridge}(optimizer) + bridge = PolyJuMP.Bridges.Objective.ToReformulatedPolynomialBridge{Float64} + optimizer = MOI.Bridges.Objective.SingleBridgeOptimizer{bridge}(optimizer) + end + return optimizer +end + function MOI.get( model::Optimizer{T}, attr::MOI.Bridges.ListOfNonstandardBridges, diff --git a/src/nl_to_polynomial.jl b/src/nl_to_polynomial.jl index 964c267..f69b07f 100644 --- a/src/nl_to_polynomial.jl +++ b/src/nl_to_polynomial.jl @@ -90,8 +90,9 @@ function _to_polynomial!( return a - b elseif _is_operator(expr, :^) && length(operands) == 2 a, b = _to_polynomial!.(Ref(d), T, operands) - if !(b isa Integer && round(Int, b) == b) - b = round(Int, b) + b_int = convert(Int, b) + if !(b isa Integer && b_int == b) + b = b_int end return a^b elseif _is_operator(expr, :/) && length(operands) == 2 diff --git a/test/qcqp.jl b/test/qcqp.jl index f190061..5d695a3 100644 --- a/test/qcqp.jl +++ b/test/qcqp.jl @@ -393,6 +393,39 @@ function test_constraint_bridging_cost_polynomial(x, y, T) return end +function test_reformulated_polynomial_bridge(_, _, T) + # The bounds of the auxiliary variables are floating point + T <: AbstractFloat || return + inner = MOI.Utilities.Model{T}() + B = PolyJuMP.Bridges.Constraint.ToReformulatedPolynomialBridge{T} + model = MOI.Bridges.Constraint.SingleBridgeOptimizer{B}(inner) + a, _ = MOI.add_constrained_variable(model, MOI.Interval{T}(0, 1)) + b, _ = MOI.add_constrained_variable(model, MOI.Interval{T}(0, 1)) + sum_ab = MOI.ScalarNonlinearFunction(:+, Any[a, b]) + pow(f, n) = MOI.ScalarNonlinearFunction(:^, Any[f, n]) + # Number of auxiliary variables created when bridging `f` + function num_aux(f) + n = MOI.get(inner, MOI.NumberOfVariables()) + MOI.add_constraint(model, f, MOI.LessThan(one(T))) + return MOI.get(inner, MOI.NumberOfVariables()) - n + end + # Degree 3 with 2 variables: expanding is cheap, nothing is substituted + @test num_aux(pow(sum_ab, 3)) == 0 + # Degree 4: the base is replaced by `y ∈ [0, 2]` with `a + b - y = 0` + @test num_aux(pow(sum_ab, 4)) == 1 + y = last(MOI.get(inner, MOI.ListOfVariableIndices())) + @test MOI.Utilities.get_bounds(inner, T, y) == (0, 2) + F = MOI.ScalarNonlinearFunction + @test MOI.get(inner, MOI.NumberOfConstraints{F,MOI.EqualTo{T}}()) == 1 + # Product of two sums: only one of them is substituted + mul = MOI.ScalarNonlinearFunction(:*, Any[sum_ab, sum_ab, a]) + @test num_aux(mul) == 1 + # The same sum appearing twice in a function reuses its auxiliary variable + add = MOI.ScalarNonlinearFunction(:+, Any[pow(sum_ab, 4), pow(sum_ab, 5)]) + @test num_aux(add) == 1 + return +end + function runtests(x, y) for name in names(@__MODULE__; all = true) if startswith("$name", "test_")