From 3c1d94142f2a086f5e919e0603d5ce7bc48816a6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 17 Jun 2026 15:45:51 +0200 Subject: [PATCH 01/17] Use Image bridge for constraints with domain too --- src/Bridges/Constraint/image.jl | 104 +++--- .../sos_polynomial_in_semialgebraic_set.jl | 330 +++++++++--------- 2 files changed, 207 insertions(+), 227 deletions(-) diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index cb47bf4c9..b77d2d2e7 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -69,7 +69,9 @@ struct ImageBridge{ constraints::Vector{MOI.ConstraintIndex{F}} zero_constraint::Union{Nothing,MOI.ConstraintIndex{G,MOI.Zeros}} set::SOS.WeightedSOSCone{M} - first::Vector{Union{Nothing,Int}} + # `first[t] = (gram_basis_idx, k_local)` identifies which entry in which + # gram basis anchors the `t`-th basis monomial. + first::Vector{Union{Nothing,Tuple{Int,Int}}} end function MOI.Bridges.Constraint.bridge_constraint( @@ -80,27 +82,33 @@ function MOI.Bridges.Constraint.bridge_constraint( ) where {T,F,G,M} @assert MOI.output_dimension(g) == length(set.basis) scalars = MOI.Utilities.scalarize(g) - # `found[mono] = k` records the index in the current `f` where the - # monomial `mono` has been anchored, with factor `factor[mono]`. - found = Dict{eltype(set.basis),Int}() - factor = Dict{eltype(set.basis),T}() - # `pending[t]` lists `(k, factor)` contributions to the `t`-th basis - # monomial that is not yet anchored; the contributions will be absorbed - # when the monomial gets anchored, or otherwise turned into a zero - # constraint. - pending = Dict{Int,Vector{Tuple{Int,T}}}() - first = Union{Nothing,Int}[nothing for _ in eachindex(scalars)] + # `found[mono] = (i_basis, k, factor)` records that the monomial `mono` + # is anchored by the `k`-th entry of the `i_basis`-th gram matrix, with + # the recorded factor. + found = Dict{eltype(set.basis),Tuple{Int,Int,T}}() + # `pending[t]` lists `(i_basis, k, factor)` contributions to the `t`-th + # basis monomial that is not yet anchored; the contributions will be + # absorbed when the monomial gets anchored, or otherwise turned into a + # zero constraint. + pending = Dict{Int,Vector{Tuple{Int,Int,T}}}() + first = Union{Nothing,Tuple{Int,Int}}[nothing for _ in eachindex(scalars)] variables = MOI.VariableIndex[] - constraints = MOI.ConstraintIndex{F}[] - # `fs[i]` keeps the function passed to `MOI.add_constraint` for the - # `i`-th gram basis. It is used below to build the residual zero - # constraint for unanchored basis monomials. + # The gram matrices are built into `fs` first, then handed to MOI at + # the end: when an entry from a later gram basis slacks against the + # anchor of an earlier basis, the earlier `f` still has to be mutable. fs = F[] - for (gram_basis, weight) in zip(set.gram_bases, set.weights) + for gram_basis in set.gram_bases cone = SOS.matrix_cone(M, length(gram_basis)) - f = MOI.Utilities.zero_with_output_dimension(F, MOI.dimension(cone)) + push!( + fs, + MOI.Utilities.zero_with_output_dimension(F, MOI.dimension(cone)), + ) + end + for (i_basis, (gram_basis, weight)) in + enumerate(zip(set.gram_bases, set.weights)) weight_basis = SA.basis(weight) weight_coeffs = SA.coeffs(weight) + f = fs[i_basis] k = 0 for j in eachindex(gram_basis) for i in 1:j @@ -115,25 +123,18 @@ function MOI.Bridges.Constraint.bridge_constraint( mono = mono_w * SA.star(gram_basis[i]) * gram_basis[j] f_kw = w_coef * diag_factor if haskey(found, mono) - k_a = found[mono] - f_a = factor[mono] + i_a, k_a, f_a = found[mono] + f_a_func = fs[i_a] if entry_anchored - # Entry `k` is anchored by a previous weight term: - # subtract its current contribution from the - # anchor of `mono`. f_k = MOI.Utilities.eachscalar(f)[k] MOI.Utilities.operate_output_index!( -, T, k_a, - f, + f_a_func, (f_kw / f_a) * f_k, ) else - # Entry `k` is free so far: introduce a slack - # variable representing its value (or extend the - # one already introduced by an earlier weight - # term) and adjust the anchor. if slack_var === nothing slack_var = MOI.add_variable(model) push!(variables, slack_var) @@ -149,17 +150,15 @@ function MOI.Bridges.Constraint.bridge_constraint( +, T, k_a, - f, + f_a_func, (f_kw / f_a) * slack_var, ) end elseif mono in set.basis t = set.basis[mono] if !entry_anchored && slack_var === nothing - # Anchor `mono` at entry `k` with factor `f_kw`. - found[mono] = k - factor[mono] = f_kw - first[t] = k + found[mono] = (i_basis, k, f_kw) + first[t] = (i_basis, k) MOI.Utilities.operate_output_index!( +, T, @@ -168,8 +167,9 @@ function MOI.Bridges.Constraint.bridge_constraint( inv(f_kw) * scalars[t], ) if haskey(pending, t) - for (k_p, f_p) in pending[t] - f_kp = MOI.Utilities.eachscalar(f)[k_p] + for (i_p, k_p, f_p) in pending[t] + f_kp = + MOI.Utilities.eachscalar(fs[i_p])[k_p] MOI.Utilities.operate_output_index!( -, T, @@ -182,37 +182,32 @@ function MOI.Bridges.Constraint.bridge_constraint( end entry_anchored = true else - # Entry `k` is busy (already anchored or slack): - # record the contribution for when `mono` gets - # anchored or for a zero constraint. - push!(get!(pending, t, Tuple{Int,T}[]), (k, f_kw)) + push!( + get!(pending, t, Tuple{Int,Int,T}[]), + (i_basis, k, f_kw), + ) end end end end end - push!(fs, f) - push!(constraints, MOI.add_constraint(model, f, cone)) end - # Build zero constraints for unanchored basis monomials. For each such - # monomial, the polynomial coefficient must equal the sum of the pending - # contributions; if no contributions exist at all, the coefficient must - # simply be zero. + constraints = MOI.ConstraintIndex{F}[] + for (i_basis, gram_basis) in enumerate(set.gram_bases) + cone = SOS.matrix_cone(M, length(gram_basis)) + push!(constraints, MOI.add_constraint(model, fs[i_basis], cone)) + end + # Build zero constraints for unanchored basis monomials. z = findall(isnothing, first) zero_terms = if isempty(z) empty(scalars) else - # The per-gram-basis `f`s above were built with a local `k`. With a - # single gram basis (the typical case) the global pending index is - # the local one. Supporting multiple gram bases would require - # tracking the originating basis for each pending contribution. - @assert length(fs) == 1 - f_local = MOI.Utilities.eachscalar(fs[1]) map(z) do t term = scalars[t] if haskey(pending, t) - for (k_p, f_p) in pending[t] - term = MA.operate!!(-, term, f_p * f_local[k_p]) + for (i_p, k_p, f_p) in pending[t] + f_kp = MOI.Utilities.eachscalar(fs[i_p])[k_p] + term = MA.operate!!(-, term, f_p * f_kp) end end return term @@ -373,11 +368,12 @@ function MOI.get( z_idx += 1 return z[z_idx] else + i_basis, k = bridge.first[i] f = MOI.Utilities.filter_variables( !Base.Fix2(in, bridge.variables), - funcs[1][bridge.first[i]], # FIXME + funcs[i_basis][k], ) - if !MOI.Utilities.is_diagonal_vectorized_index(bridge.first[i]) + if !MOI.Utilities.is_diagonal_vectorized_index(k) f = T(2) * f end return f diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 51ff1da12..8b557dae0 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -1,50 +1,31 @@ -function lagrangian_multiplier( - model::MOI.ModelLike, - certificate, - index, - preprocessed, - T::Type, -) - basis = Certificate.multiplier_basis(certificate, index, preprocessed) - MCT = SOS.matrix_cone_type(typeof(certificate)) - return SOS.add_gram_matrix(model, MCT, basis, T)..., basis -end - struct SOSPolynomialInSemialgebraicSetBridge{ T, F<:MOI.AbstractVectorFunction, DT<:SemialgebraicSets.AbstractSemialgebraicSet, CT<:Certificate.AbstractIdealCertificate, - B<:Union{Vector{<:MB.SA.ExplicitBasis},MB.SA.ExplicitBasis}, - UMCT<:Union{ - Vector{<:MOI.ConstraintIndex{MOI.VectorOfVariables}}, - MOI.ConstraintIndex{MOI.VectorOfVariables}, - }, - UMST, - MCT, BT<:MB.SubBasis{MB.Monomial}, + M, + NB<:SA.ExplicitBasis, + GB<:SA.ExplicitBasis, + W<:SA.AlgebraElement, } <: MOI.Bridges.Constraint.AbstractBridge - lagrangian_bases::Vector{B} - lagrangian_variables::Vector{ - Union{Vector{MOI.VariableIndex},Vector{Vector{MOI.VariableIndex}}}, - } - lagrangian_constraints::Vector{UMCT} - constraint::MOI.ConstraintIndex{F,SOS.SOSPolynomialSet{DT,BT,CT}} + constraint::MOI.ConstraintIndex{F,SOS.WeightedSOSCone{M,NB,GB,W}} basis::BT + # `sigma_0_indices` records the range of `gram_bases` that correspond to + # σ_0 (possibly multiple if the certificate uses sparsity). The + # remaining gram bases (after this range) are the Lagrangian multipliers + # σ_1, σ_2, … in order. + sigma_0_indices::Union{Int,UnitRange{Int}} end function MOI.Bridges.Constraint.bridge_constraint( - ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}}, + ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}}, model::MOI.ModelLike, f::MOI.AbstractVectorFunction, set::SOS.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}, -) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT} +) where {T,F,DT,CT,BT,M,NB,GB,W} @assert MOI.output_dimension(f) == length(set.basis) - # MOI does not modify the coefficients of the functions so we can modify `p`. - # without altering `f`. - # The monomials may be copied by MA however so we need to copy it. - # TODO remove `collect` when `DynamicPolynomials.MonomialVector` can be used as keys - p = MB.algebra_element( + poly = MB.algebra_element( SA.SparseCoefficients( copy(collect(set.basis.keys)), MOI.Utilities.scalarize(f), @@ -52,60 +33,76 @@ function MOI.Bridges.Constraint.bridge_constraint( ), MB.implicit_basis(set.basis), ) - λ_bases = B[] - λ_variables = - Union{Vector{MOI.VariableIndex},Vector{Vector{MOI.VariableIndex}}}[] - λ_constraints = UMCT[] - preprocessed = - Certificate.preprocessed_domain(set.certificate, set.domain, p) + ideal_cert = Certificate.ideal_certificate(set.certificate) + # Reduce by the algebraic ideal part `V` of the basic semialgebraic + # domain so the rest of the certificate operates modulo it. + domain_V = MP.similar(set.domain.V, T) + poly_reduced = Certificate.reduced_polynomial(ideal_cert, poly, domain_V) + sigma_0_basis = Certificate.gram_basis( + ideal_cert, + Certificate.with_variables(poly_reduced, set.domain), + ) implicit_basis = MB.implicit_basis(set.basis) - cache = zero(MOI.ScalarAffineFunction{T}, MB.algebra(implicit_basis)) + # Build `[1, g_1, …]` weights and `[σ_0_basis, σ_1_basis, …]` gram bases + # so the whole Putinar decomposition becomes a single `WeightedSOSCone`. + preprocessed = + Certificate.preprocessed_domain(set.certificate, set.domain, poly) + # Use the same `sparse_coefficients`-based construction for the unit + # weight as for the polynomial g_i weights so they all share the same + # concrete element type and fit in a homogeneous `Vector{W}`. + some_mono = first(MB.keys_as_monomials(set.basis)) + unit_poly = MP.polynomial(MP.term(one(T), MP.constant_monomial(some_mono))) + gram_bases_raw = Any[sigma_0_basis] + weights_raw = W[ + MB.algebra_element( + MB.sparse_coefficients(unit_poly), + implicit_basis, + ), + ] for index in Certificate.preorder_indices(set.certificate, preprocessed) - λ, λ_variable, λ_constraint, λ_basis = lagrangian_multiplier( - model, - set.certificate, - index, - preprocessed, - T, + push!( + gram_bases_raw, + Certificate.multiplier_basis(set.certificate, index, preprocessed), ) - push!(λ_variables, λ_variable) - push!(λ_constraints, λ_constraint) - push!(λ_bases, λ_basis) - # As `*(::MOI.ScalarAffineFunction{T}, ::S)` is only defined if `S == T`, we - # need to call `similar`. This is critical since `T` is - # `Float64` when used with JuMP and the coefficient type is often `Int` if - # `set.domain.V` is `FullSpace` or `FixedPolynomialSet`. g = Certificate.generator(set.certificate, index, preprocessed) - MA.operate_to!(cache, +, λ) - # TODO replace with `MA.sub_mul` when it works. - p = MA.operate!( - SA.UnsafeAddMul(*), - p, - cache, + push!( + weights_raw, MB.algebra_element( - MB.sparse_coefficients(-one(T) * similar(g, T)), + MB.sparse_coefficients(one(T) * similar(g, T)), implicit_basis, ), ) end - MA.operate!(SA.canonical, SA.coeffs(p)) - new_set = SOS.SOSPolynomialSet( - set.domain.V, - MB.explicit_basis(p), - Certificate.ideal_certificate(set.certificate), + # `gram_basis` returns either a single basis or a `Vector{basis}` when + # sparsity is in play. Reuse `_flatten` from the SOSPolynomial bridge to + # turn the latter into a uniform `Vector{GB}` (replicating the matching + # weight for each sparsity block). The σ_0 block is the first entry, so + # `sigma_0_indices` is what `_flatten` reports for that single entry. + sigma_0_gram_bases, _, sigma_0_indices = _flatten( + identity.([sigma_0_basis]), + [weights_raw[1]], ) + gram_bases, weights, _ = _flatten( + identity.(gram_bases_raw), + weights_raw, + ) + new_basis = Certificate.zero_basis( + ideal_cert, + MB.explicit_basis(poly_reduced), + domain_V, + gram_bases, + weights, + ) + new_coeffs = SA.coeffs(poly_reduced, new_basis) constraint = MOI.add_constraint( model, - MOI.Utilities.vectorize(SA.values(SA.coeffs(p))), - new_set, + MOI.Utilities.vectorize(SA.values(new_coeffs)), + SOS.WeightedSOSCone{M}(new_basis, gram_bases, weights), ) - - return SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}( - λ_bases, - λ_variables, - λ_constraints, + return SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}( constraint, set.basis, + sigma_0_indices, ) end @@ -117,14 +114,14 @@ function MOI.supports_constraint( return true end function MOI.Bridges.added_constrained_variable_types( - ::Type{<:SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT}}, -) where {T,F,DT,CT} - return constrained_variable_types(SOS.matrix_cone_type(CT)) + ::Type{<:SOSPolynomialInSemialgebraicSetBridge}, +) + return Tuple{Type}[] end function MOI.Bridges.added_constraint_types( - ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}}, -) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT} - return [(F, SOS.SOSPolynomialSet{DT,BT,CT})] + ::Type{SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}}, +) where {T,F,DT,CT,BT,M,NB,GB,W} + return Tuple{Type,Type}[(F, SOS.WeightedSOSCone{M,NB,GB,W})] end function MOI.Bridges.Constraint.concrete_bridge_type( ::Type{<:SOSPolynomialInSemialgebraicSetBridge{T}}, @@ -137,89 +134,77 @@ function MOI.Bridges.Constraint.concrete_bridge_type( }, }, ) where {T,S,PS,AT,CT,BT<:MB.MonomialIndexedBasis{MB.Monomial}} - - # promotes VectorOfVariables into VectorAffineFunction, it should be enough - # for most use cases G = MOI.Utilities.promote_operation(-, T, F, MOI.VectorOfVariables) - MCT = SOS.matrix_cone_type(CT) - MT = MP.monomial_type(BT) - B = Certificate.multiplier_basis_type(CT, MT) - UMCT = union_constraint_types(MCT) - UMST = union_set_types(MCT) IC = Certificate.ideal_certificate(CT) - return SOSPolynomialInSemialgebraicSetBridge{T,G,AT,IC,B,UMCT,UMST,MCT,BT} + M = SOS.matrix_cone_type(IC) + # The multiplier gram basis of the inner ideal certificate is used both + # for σ_0 and (via `multiplier_basis_type`) for each Lagrangian σ_i. + GB = + SOSPolynomialBridgeType_helper_gram_basis_type(IC, MP.monomial_type(BT)) + # All weights (the σ_0 unit weight and each g_i) are built with + # `MB.sparse_coefficients` so they share a single concrete type with + # `Vector{Int}` keys and `Vector{T}` values. + A = MA.promote_operation( + MB.algebra, + MA.promote_operation(MB.implicit_basis, BT), + ) + C = SA.SparseCoefficients{ + Vector{Int}, + T, + Vector{Vector{Int}}, + Vector{T}, + MP.Graded{MP.LexOrder}, + } + W = SA.AlgebraElement{T,A,C} + NB = MA.promote_operation( + Certificate.zero_basis, + IC, + BT, + SemialgebraicSets.similar_type(AT, T), + Vector{GB}, + Vector{W}, + ) + return SOSPolynomialInSemialgebraicSetBridge{T,G,AT,IC,BT,M,NB,GB,W} +end + +# A small helper to keep the `concrete_bridge_type` readable. +function SOSPolynomialBridgeType_helper_gram_basis_type(::Type{IC}, ::Type) where {IC} + return _eltype(MA.promote_operation(SOS.Certificate.gram_basis, IC)) end -# Attributes, Bridge acting as an model +# Attributes, Bridge acting as a model function MOI.get( - bridge::SOSPolynomialInSemialgebraicSetBridge, + ::SOSPolynomialInSemialgebraicSetBridge, ::MOI.NumberOfVariables, ) - return mapreduce(_num_variables, +, bridge.lagrangian_variables, init = 0) + return 0 end function MOI.get( - bridge::SOSPolynomialInSemialgebraicSetBridge, + ::SOSPolynomialInSemialgebraicSetBridge, ::MOI.ListOfVariableIndices, ) - return Iterators.flatten([ - _list_variables(Qi) for Qi in bridge.lagrangian_variables - ]) -end -function MOI.get( - bridge::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST}, - ::MOI.NumberOfConstraints{MOI.VectorOfVariables,S}, -) where {T,F,DT,CT,B,UMCT,UMST,S<:UMST} - return mapreduce( - cQ -> - _num_constraints(cQ, MOI.ConstraintIndex{MOI.VectorOfVariables,S}), - +, - bridge.lagrangian_constraints, - init = 0, - ) + return MOI.VariableIndex[] end function MOI.get( - ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST}, - ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,S}, -) where {T,F,DT,CT,B,UMCT,UMST,S<:UMST} - C = MOI.ConstraintIndex{MOI.VectorOfVariables,S} - return Iterators.flatten([ - _list_constraints(cQ, C) for cQ in bridge.lagrangian_constraints - ]) -end -function MOI.get( - ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}, - ::MOI.NumberOfConstraints{F,SOS.SOSPolynomialSet{DT,BT,CT}}, -) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT} + ::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}, + ::MOI.NumberOfConstraints{F,SOS.WeightedSOSCone{M,NB,GB,W}}, +) where {T,F,DT,CT,BT,M,NB,GB,W} return 1 end function MOI.get( - b::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,B,UMCT,UMST,MCT,BT}, - ::MOI.ListOfConstraintIndices{F,SOS.SOSPolynomialSet{DT,BT,CT}}, -) where {T,F,DT,CT,B,UMCT,UMST,MCT,BT} + b::SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}, + ::MOI.ListOfConstraintIndices{F,SOS.WeightedSOSCone{M,NB,GB,W}}, +) where {T,F,DT,CT,BT,M,NB,GB,W} return [b.constraint] end # Indices -function _delete_variables(model, Q::Vector{MOI.VariableIndex}) - if !isempty(Q) - # FIXME Since there is not variables in the list, we cannot - # identify the `EmptyBridge` to delete - MOI.delete(model, Q) - end -end -function _delete_variables(model, Qs::Vector{Vector{MOI.VariableIndex}}) - for Q in Qs - _delete_variables(model, Q) - end -end function MOI.delete( model::MOI.ModelLike, bridge::SOSPolynomialInSemialgebraicSetBridge, ) MOI.delete(model, bridge.constraint) - for variables in bridge.lagrangian_variables - _delete_variables(model, variables) - end + return end # Attributes, Bridge acting as a constraint @@ -251,57 +236,56 @@ function MOI.get( attr::PolyJuMP.MomentsAttribute, bridge::SOSPolynomialInSemialgebraicSetBridge, ) - return MOI.get(model, attr, bridge.constraint) + # We can't forward `MomentsAttribute` to `bridge.constraint` because + # nothing in the `WeightedSOSCone` chain knows about + # `MomentsAttribute`. Build it directly from the dual and the new basis + # stored on the `WeightedSOSCone` set. + set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint) + return MultivariateMoments.moment_vector( + MOI.get(model, MOI.ConstraintDual(attr.result_index), bridge.constraint), + set.basis, + ) end function MOI.get( model::MOI.ModelLike, - attr::Union{ - SOS.CertificateBasis, - SOS.GramMatrixAttribute, - SOS.MomentMatrixAttribute, - }, + attr::SOS.CertificateBasis, bridge::SOSPolynomialInSemialgebraicSetBridge, ) return MOI.get(model, attr, bridge.constraint) end -function _gram( - f::Function, - Q::Vector{MOI.VariableIndex}, - gram_basis, - T::Type, - MCT, +function MOI.get( + model::MOI.ModelLike, + attr::Union{SOS.GramMatrixAttribute,SOS.MomentMatrixAttribute}, + bridge::SOSPolynomialInSemialgebraicSetBridge, ) - return SOS.build_gram_matrix(convert(Vector{T}, f(Q)), gram_basis, MCT, T) + # `multiplier_index = 0` means "σ_0"; remap it to whichever indices the + # σ_0 block(s) occupy in the underlying `WeightedSOSCone`'s + # `gram_bases`. + SOS.check_multiplier_index_bounds(attr, 0:0) + return _get(model, attr, bridge.constraint, bridge.sigma_0_indices) end -function _gram( - f::Function, - Qs::Vector{Vector{MOI.VariableIndex}}, - gram_bases, - T::Type, - MCT, -) - return SOS.build_gram_matrix(gram_bases, MCT, T) do i - return convert(Vector{T}, f(Qs[i])) - end -end function MOI.get( model::MOI.ModelLike, attr::SOS.LagrangianMultipliers, - bridge::SOSPolynomialInSemialgebraicSetBridge{T,F,CT,B,UMCT,UMST,MCT,M}, -) where {T,F,CT,B,UMCT,UMST,MCT,M} - @assert eachindex(bridge.lagrangian_variables) == - eachindex(bridge.lagrangian_bases) - return map( - i -> _gram( - Q -> MOI.get(model, MOI.VariablePrimal(attr.result_index), Q), - bridge.lagrangian_variables[i], - bridge.lagrangian_bases[i], - T, - MCT, - ), - eachindex(bridge.lagrangian_variables), - ) + bridge::SOSPolynomialInSemialgebraicSetBridge, +) + # σ_0 occupies the `sigma_0_indices` slot(s); the remaining gram bases + # are σ_1, σ_2, … (one Lagrangian multiplier per inequality). + set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint) + n_total = length(set.gram_bases) + sigma_0_max = bridge.sigma_0_indices isa Int ? bridge.sigma_0_indices : + last(bridge.sigma_0_indices) + return map((sigma_0_max + 1):n_total) do i + return MOI.get( + model, + SOS.GramMatrixAttribute(; + multiplier_index = i, + result_index = attr.result_index, + ), + bridge.constraint, + ) + end end From df3d173f6ec2b35cbce2bacf726284e9c988a413 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 08:51:23 +0200 Subject: [PATCH 02/17] Fix --- src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 8b557dae0..3144f0ac5 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -153,7 +153,7 @@ function MOI.Bridges.Constraint.concrete_bridge_type( T, Vector{Vector{Int}}, Vector{T}, - MP.Graded{MP.LexOrder}, + MP.ordering(BT), } W = SA.AlgebraElement{T,A,C} NB = MA.promote_operation( From 056819ca76ab3eaf9fd6d8df20cc23ecd8cef56e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 08:51:40 +0200 Subject: [PATCH 03/17] Add examples for ImageBridge --- src/Bridges/Constraint/image.jl | 80 +++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index b77d2d2e7..21eff4cc8 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -39,6 +39,86 @@ The gram matrix is therefore: \\end{bmatrix} ``` +### Non-unit weight + +When the weight `w` of the gram basis is not one, each entry of `b * b'` +contributes to a polynomial coefficient with an additional factor `w`. +Concretely, every `(i, j)` entry of the gram matrix gets divided by `w` +(and by the usual extra `2` for off-diagonal entries), while the slack +adjustment ratio `(factor / factor_anchor)` is unchanged since `w` cancels. + +Take the same polynomial `p` and gram basis as above, but now write the +SOS constraint with weight `w = 1 + x` and `q = q_0 + q_1 x`, where each `q_i` +is a coefficient of the gram matrix. +For brevity, denote the polynomial coefficients +`p = a_0 + a_1 x + a_2 x^2 + a_3 x^3` with gram basis `b = [1, x]`. +The product `b * b'` is +```math +\\begin{bmatrix} +1 & x\\\\ +x & x^2 +\\end{bmatrix} +``` +and `(1 + x) * b * b'` distributes one extra power of `x`: +the `(1, 1)` entry contributes both to the `1` and the `x` coefficient of `p`, +the `(1, 2)` (and `(2, 1)`) entry contributes to `x` and `x^2`, and the +`(2, 2)` entry contributes to `x^2` and `x^3`. +Anchoring each monomial of `p` to the first gram entry that produces it +(while accumulating the contributions of every earlier entry) gives: +- `q_{1,1} = a_0`; +- `q_{1,2}` is anchored on the `x` coefficient and absorbs `q_{1,1}`'s + contribution, yielding `q_{1,2} = (a_1 - a_0) / 2`; +- `q_{2,2}` is anchored on the `x^2` coefficient and absorbs the contribution + of `q_{1,2}`, yielding `q_{2,2} = a_2 - a_1 + a_0`. +There is no entry left to anchor the `x^3` coefficient, so its equation +`a_3 = q_{2,2}` becomes the zero constraint +`a_0 - a_1 + a_2 - a_3 = 0`. +The gram matrix is therefore: +```math +\\begin{bmatrix} +a_0 & \\frac{a_1 - a_0}{2}\\\\ +\\frac{a_1 - a_0}{2} & a_0 - a_1 + a_2 +\\end{bmatrix} +``` +with the additional zero constraint `a_0 - a_1 + a_2 - a_3 = 0`. + +### Multiple weights and gram bases + +The same logic extends to a `WeightedSOSCone` with several gram bases +`b_i` and weights `w_i`: the bridge walks through every `(i, j, k)` entry +of every gram matrix and treats each weight monomial as an additional +contribution to the same anchoring/slack bookkeeping. Anchors and slack +variables can therefore span across different gram matrices. + +Consider for instance Putinar-style decomposition +`p = σ_0 + x * σ_1` with +`p = a_0 + a_1 x + a_2 x^2`, +weights `w_0 = 1`, `w_1 = x`, +and gram bases `b_0 = [1, x]` (a 2 × 2 gram matrix `Q`) +and `b_1 = [1]` (a 1 × 1 gram matrix `R`). + +Walking through `b_0`'s entries first (weight `1`): +- `Q_{1,1}` anchors `1`, so `Q_{1,1} = a_0`; +- `Q_{1,2}` anchors `x` and `Q_{1,2} = a_1 / 2`; +- `Q_{2,2}` anchors `x^2` and `Q_{2,2} = a_2`. + +Then `R_{1,1}` with weight `x` produces the monomial +`x · 1 · 1 = x`, which is already anchored by `Q_{1,2}` (off-diagonal, +factor `2`). Introducing a slack variable `λ` to balance, with adjustment +ratio `1 / 2`: +`R_{1,1} = -λ` and `Q_{1,2}` is updated to `(a_1 + λ) / 2`. +The gram matrices are therefore +```math +Q = \\begin{bmatrix} +a_0 & \\frac{a_1 + \\lambda}{2}\\\\ +\\frac{a_1 + \\lambda}{2} & a_2 +\\end{bmatrix}, +\\qquad +R = \\begin{bmatrix}-\\lambda\\end{bmatrix}. +``` +The slack `λ` is constrained only by `Q ⪰ 0` and `R ⪰ 0`, +i.e. by the PSD constraints on the gram matrices. + ## Source node `ImageBridge` supports: From a7a7af81a9aa3cc4e694d4b614b847cda42d7ada Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 09:37:23 +0200 Subject: [PATCH 04/17] Add bridge section in docs --- docs/make.jl | 2 + docs/src/bridges.md | 212 ++++++++++++++++++ docs/src/index.md | 2 +- docs/src/reference/bridges.md | 56 +++++ docs/src/reference/standard_form.md | 25 +-- src/Bridges/Constraint/diagonally_dominant.jl | 28 +++ src/Bridges/Constraint/empty.jl | 18 ++ src/Bridges/Constraint/psd2x2.jl | 25 +++ src/Bridges/Constraint/sos_polynomial.jl | 26 +++ .../sos_polynomial_in_semialgebraic_set.jl | 30 +++ src/Bridges/Variable/copositive_inner.jl | 27 +++ src/Bridges/Variable/kernel.jl | 39 ++++ src/Bridges/Variable/lowrank.jl | 39 ++++ src/Bridges/Variable/psd2x2.jl | 22 ++ .../Variable/scaled_diagonally_dominant.jl | 34 ++- 15 files changed, 550 insertions(+), 35 deletions(-) create mode 100644 docs/src/bridges.md create mode 100644 docs/src/reference/bridges.md diff --git a/docs/make.jl b/docs/make.jl index 2d510ca71..7dfca361a 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -63,9 +63,11 @@ makedocs( "Sum-of-Squares Programming" => "sumofsquares.md", "Variables" => "variables.md", "Constraints" => "constraints.md", + "Bridging mechanism" => "bridges.md", "API Reference" => [ "reference/standard_form.md", "reference/constraints.md", + "reference/bridges.md", "reference/certificate.md", "reference/internal.md", ], diff --git a/docs/src/bridges.md b/docs/src/bridges.md new file mode 100644 index 000000000..a18229ada --- /dev/null +++ b/docs/src/bridges.md @@ -0,0 +1,212 @@ +# Bridging mechanism + +[Bridges](https://jump.dev/MathOptInterface.jl/stable/submodules/Bridges/overview/) +are MathOptInterface's mechanism for taking a constraint or constrained +variable expressed in one set and rewriting it in terms of other sets, until +something that the underlying solver supports natively is reached. + +In SumOfSquares, the user-facing object — a JuMP `@constraint` such as +```julia +@constraint(model, p in SOSCone(), domain = (@set x^2 + y^2 == 1 && x >= 0)) +``` +is several bridge hops away from the cone that the solver will actually see. +This page walks through those hops, so that you can predict — and tune — the +chain of reformulations between the SOS constraint that you wrote and the +PSD (or rotated-second-order-cone, or rank-1) constraint that hits the +solver. + +The high-level shape of the chain is: + +```text +SOSPolynomialSet ────► WeightedSOSCone ────► (PSD-like cone) ────► solver + ▲ ▲ + constraint variable + bridges bridges +``` + +The first hop is always a *constraint bridge*. The second hop, which produces +the gram-matrix PSD constraints, can be either a *variable bridge* (e.g. +[`SumOfSquares.Bridges.Variable.KernelBridge`](@ref)) or a *constraint +bridge* (e.g. [`SumOfSquares.Bridges.Constraint.ImageBridge`](@ref)). The +choice is driven by which form the underlying solver supports more directly, +and the bridge graph picks the cheapest path. + +## From a JuMP constraint to a `WeightedSOSCone` + +[`SumOfSquares.SOSPolynomialSet`](@ref) is the MOI set that JuMP creates +when you ask for an SOS constraint over a domain. The domain is carried as +a parameter: + + * `SOSPolynomialSet{<:SemialgebraicSets.AbstractAlgebraicSet}` — equality + constraints only (or no `domain` at all); + * `SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}` — at least + one inequality constraint. + +Two bridges turn this into a [`SumOfSquares.WeightedSOSCone`](@ref): + + * [`SumOfSquares.Bridges.Constraint.SOSPolynomialBridge`](@ref) handles the + algebraic case. It reduces the polynomial modulo the algebraic ideal, + asks the certificate for a single gram basis, and produces a + `WeightedSOSCone` with weight ``1``. + * [`SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge`](@ref) + handles the basic-semialgebraic case (Putinar-style). It collects one + multiplier basis per inequality ``g_i`` and produces a single + `WeightedSOSCone` with weights ``[1, g_1, g_2, \ldots]`` and gram bases + ``[\sigma_0, \sigma_1, \sigma_2, \ldots]``. The σ_0 basis comes from the + inner ideal certificate. + +After this hop, every SOS constraint is a single `WeightedSOSCone` +constraint, regardless of whether the user supplied a domain. + +## From a `WeightedSOSCone` to PSD + +There are two natural ways to lower a `WeightedSOSCone`: + + * **Variable side** — add the gram matrices ``Q_i`` as PSD-constrained + variables and express each polynomial coefficient as an affine function + of the ``Q_i``. This is what + [`SumOfSquares.Bridges.Variable.KernelBridge`](@ref) does (and what + [`SumOfSquares.Bridges.Variable.LowRankBridge`](@ref) does for the + rank-1 case, see below). + + * **Constraint side** — pick anchor entries of each ``Q_i`` so that they + are *fixed by* the polynomial coefficients (no new variables), introduce + slack variables only when two gram entries produce the same monomial, + and finally add the matrices as PSD constraints. This is what + [`SumOfSquares.Bridges.Constraint.ImageBridge`](@ref) does; the + docstring contains a worked example. + +Mathematically, if ``Σ`` is the SOS cone and ``S`` the PSD cone, then +``Σ = A(S)`` for a linear map ``A``. The variable bridge gives the user +``A(q)`` directly; the constraint bridge gives the user the preimage +``q \in A^{-1}(p) = A^\dagger p + \ker(A)`` instead. + +Which one is cheaper depends on the solver: + +| Solver style | Wins | +| -------------------------------------------------------------------------------------------- | ---- | +| Native PSD-constrained variables (e.g. Mosek's `barvar`, Hypatia) | `KernelBridge` | +| PSD as a constraint only (e.g. Clarabel, SCS) | `ImageBridge` | +| Specialised low-rank PSD format (e.g. Loraine, Burer–Monteiro) | `LowRankBridge` + LowRankOpt bridges | + +The bridge graph uses `MOI.Bridges.bridging_cost` to pick. For the SOS +pipeline, `ImageBridge` carries an explicit `bridging_cost(::Type{<:ImageBridge}) = 2.0` +so that `KernelBridge` wins for solvers that support constrained PSD +variables (and vice versa for solvers that only support PSD as a +constraint). + +## How Hypatia gets specialised cones + +Hypatia natively supports `MOI.PositiveSemidefiniteConeTriangle` *and* a +range of more specialised cones (`HypatiaSOSCone`, etc.). Because the SOS +pipeline lands in `WeightedSOSCone` after the first hop, +[`SumOfSquares.Bridges.Variable.KernelBridge`](@ref) wins (Hypatia accepts +PSD as constrained variables) and the chain stops there: + +```text +SOSPolynomialSet + └─SOSPolynomialBridge──► WeightedSOSCone + └─KernelBridge──► PositiveSemidefiniteConeTriangle (Hypatia native) +``` + +When the cone is structurally special (Lagrange basis with rank-1 factors, +see next section), the chain takes a different exit through +`LowRankBridge` and Hypatia still ends up receiving a cone it understands +natively — the bridge graph is what makes this transparent. + +## How LowRankOpt exploits rank-1 structure + +When the gram basis is a [`MultivariateBases.LagrangeBasis`](@ref) — the +polynomial is sampled at a fixed set of points — the relation between +gram-matrix entries and polynomial coefficients factorises through a list +of rank-1 outer products, one per Lagrange node: + +```math +p_j = \sum_i w_i(x_j) \langle u_{i,j} u_{i,j}^\top, Q_i \rangle +``` + +where ``u_{i,j}`` is a column of the basis transformation. The +[LowRankOpt](https://github.com/blegat/LowRankOpt.jl) package provides MOI +sets that carry this rank-1 structure all the way down to the solver: + + * `LRO.SetDotProducts{W,S,V}` — the set of vectors + ``(\langle a_1, x \rangle, \ldots, \langle a_m, x \rangle)`` for ``x`` + in some PSD-like set `S`. The factors `a_k` are exposed as + `LRO.TriangleVectorization{LRO.Factorization{...}}`, i.e. rank-1 + matrices given by their factor and weight rather than by their full + matrix. + * `LRO.LinearCombinationInSet{W,S,V}` — the dual side: vectors ``y`` + such that ``\sum_k y_k a_k`` lies in ``S``. + +[`SumOfSquares.Bridges.Variable.LowRankBridge`](@ref) reformulates a +`WeightedSOSCone` whose gram basis is a `LagrangeBasis` directly into +`LRO.SetDotProducts{WITHOUT_SET}` wrapping the appropriate +`SOS.matrix_cone`, with each Lagrange node giving rise to one rank-1 +factor. So the chain looks like: + +```text +SOSPolynomialSet + └─SOSPolynomial(InSemialgebraicSet)Bridge──► WeightedSOSCone (Lagrange) + └─LowRankBridge──► SetDotProducts{WITHOUT_SET, PSD, Factorization} +``` + +From there, what happens next is solver-dependent. + +### Solvers that natively support the rank-1 format + +* **Loraine.jl** — implements an interior-point method whose Hessian + computations exploit the rank-1 structure of the constraint matrices. + It supports `LRO.SetDotProducts` directly, so the chain stops at the + output of `LowRankBridge`. + +* **LowRankOpt.BurerMonteiro** — a Burer–Monteiro-style nonlinear-program + back-end. It also consumes `LRO.SetDotProducts` directly. + +### Solvers that don't + +If the solver only knows ordinary PSD, the LowRankOpt bridge layer takes +over and unfolds the structure step by step. The most relevant LowRankOpt +bridges, in roughly the order they fire, are: + + * `LowRankOpt.Bridges.Variable.AppendSetBridge` — converts between + `SetDotProducts{WITHOUT_SET}` and `SetDotProducts{WITH_SET}` (the + `WITH_SET` form additionally carries the underlying PSD slice as + explicit variables). + * `LowRankOpt.Bridges.Variable.DotProductsBridge` — drops the rank-1 + factorisation and rewrites the dot products as explicit affine + combinations of PSD variables. This is the bridge that "loses" the + rank-1 information for solvers that can't use it. + * `LowRankOpt.Bridges.Variable.ToPositiveBridge` and + `LowRankOpt.Bridges.Variable.ToRankOneBridge` — convert between + different concrete `LRO.Factorization` types (e.g. rank-`r` factors + expressed as `r` rank-1 factors). + * `LowRankOpt.Bridges.Constraint.LinearCombinationBridge` and + `LowRankOpt.Bridges.Constraint.AppendZeroBridge` — the constraint-side + duals of the variable bridges above, used when the underlying solver + accepts the cones only as constraints. + +The net effect is that the SOS chain stays the same up to the +`LRO.SetDotProducts` level, and a fallback to plain +`MOI.PositiveSemidefiniteConeTriangle` only happens at the very end if no +solver in the stack can do better. + +## Choosing a different bridge + +The bridge graph is driven entirely by the *types* the solver supports and +by `MOI.Bridges.bridging_cost`. There are two practical ways to nudge the +chain: + + * Use `JuMP.remove_bridge(model, SomeBridge)` to disable a bridge — for + example, `JuMP.remove_bridge(model, SumOfSquares.Bridges.Variable.KernelBridge{Float64})` + forces the chain to go through `ImageBridge` even for solvers that + accept PSD-constrained variables. This is occasionally useful for + debugging or for benchmarking. + * Define a new bridge with a lower cost than the existing one. This is + how solver packages opt into specialised cones: by providing a + `MOI.supports_constrained_variable` (or `MOI.supports_constraint`) for + the special set and a bridge with a cost that makes the bridge graph + prefer the specialised route over the generic one. + +In both cases, you can inspect the resulting chain with +`MOI.Bridges.print_active_bridges(JuMP.backend(model))` to confirm the +exact path that a given constraint takes. diff --git a/docs/src/index.md b/docs/src/index.md index e7f955c0d..fa1bcb303 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -11,6 +11,6 @@ Starting with v0.8, we also have experimental support for any algebra implementi ## Contents ```@contents -Pages = ["sumofsquares.md", "variables.md", "constraints.md"] +Pages = ["sumofsquares.md", "variables.md", "constraints.md", "bridges.md"] Depth = 2 ``` diff --git a/docs/src/reference/bridges.md b/docs/src/reference/bridges.md new file mode 100644 index 000000000..81e16d9b6 --- /dev/null +++ b/docs/src/reference/bridges.md @@ -0,0 +1,56 @@ +# Bridges + +[Bridges](https://jump.dev/MathOptInterface.jl/stable/submodules/Bridges/overview/) +are the mechanism by which a Sum-of-Squares constraint is rewritten until it +hits a set that the underlying solver supports natively. See the +[Bridging mechanism](@ref) section of the manual for a high-level walk-through +of the bridge chain. + +## Utilities + +Bridges are automatically added by JuMP using the following PolyJuMP utilities: +```@docs +SumOfSquares.PolyJuMP.bridgeable +SumOfSquares.PolyJuMP.bridges +``` + +## Bridges for polynomial optimization + +```@docs +PolyJuMP.ScalarPolynomialFunction +PolyJuMP.Bridges.Objective.ToPolynomialBridge +PolyJuMP.Bridges.Constraint.ToPolynomialBridge +``` + +## Constraint bridges + +The bridges that act on the constraint side of a Sum-of-Squares constraint. + +```@docs +SumOfSquares.Bridges.Constraint.SOSPolynomialBridge +SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge +SumOfSquares.Bridges.Constraint.ImageBridge +SumOfSquares.Bridges.Constraint.PositiveSemidefinite2x2Bridge +SumOfSquares.Bridges.Constraint.DiagonallyDominantBridge +SumOfSquares.Bridges.Constraint.EmptyBridge +``` + +## Variable bridges + +The bridges that act on the variable side — they back a `WeightedSOSCone` or +one of its PSD-approximation cones with constrained variables. + +```@docs +SumOfSquares.Bridges.Variable.KernelBridge +SumOfSquares.Bridges.Variable.LowRankBridge +SumOfSquares.Bridges.Variable.PositiveSemidefinite2x2Bridge +SumOfSquares.Bridges.Variable.ScaledDiagonallyDominantBridge +SumOfSquares.Bridges.Variable.CopositiveInnerBridge +``` + +## SAGE bridges + +```@docs +SumOfSquares.PolyJuMP.SAGE.SignomialsBridge +SumOfSquares.PolyJuMP.SAGE.AGEBridge +``` diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index 1da2a687f..2d7cd6dbf 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -37,8 +37,6 @@ SAGE cones: ```@docs SumOfSquares.PolyJuMP.SAGE.Polynomials SumOfSquares.PolyJuMP.SAGE.Signomials -SumOfSquares.PolyJuMP.SAGE.SignomialsBridge -SumOfSquares.PolyJuMP.SAGE.AGEBridge ``` ## MOI Sets @@ -64,25 +62,4 @@ SumOfSquares.ScaledDiagonallyDominantConeTriangle ## Bridges -Bridges are automatically added using the following utilities: -```@docs -SumOfSquares.PolyJuMP.bridgeable -SumOfSquares.PolyJuMP.bridges -``` - -Bridges for polynomial optimization -```@docs -PolyJuMP.ScalarPolynomialFunction -PolyJuMP.Bridges.Objective.ToPolynomialBridge -PolyJuMP.Bridges.Constraint.ToPolynomialBridge -``` - -Sum-of-Squares bridges -```@docs -SumOfSquares.Bridges.Constraint.ImageBridge -``` - -Bridges for PSD cone approximations -```@docs -SumOfSquares.Bridges.Variable.ScaledDiagonallyDominantBridge -``` +See the dedicated [Bridges](@ref) reference page. diff --git a/src/Bridges/Constraint/diagonally_dominant.jl b/src/Bridges/Constraint/diagonally_dominant.jl index 5dcabc661..1c5fde127 100644 --- a/src/Bridges/Constraint/diagonally_dominant.jl +++ b/src/Bridges/Constraint/diagonally_dominant.jl @@ -1,3 +1,31 @@ +""" + DiagonallyDominantBridge{T,F,G} <: Bridges.Constraint.AbstractBridge + +`DiagonallyDominantBridge` implements a reformulation from +[`SumOfSquares.DiagonallyDominantConeTriangle`](@ref) into a set of linear +constraints. + +A diagonally dominant matrix ``Q`` satisfies ``Q_{jj} \\ge \\sum_{i \\ne j} |Q_{ij}|`` +for every row ``j``. The bridge introduces auxiliary variables ``t_{ij}`` +for each off-diagonal entry with the constraints ``t_{ij} \\ge +Q_{ij}`` and +``t_{ij} \\ge -Q_{ij}`` (so that ``t_{ij} \\ge |Q_{ij}|``), and then the row +dominance inequality ``Q_{jj} - \\sum_{i \\ne j} t_{ij} \\ge 0``. + +## Source node + +`DiagonallyDominantBridge` supports: + + * `G` in [`SumOfSquares.DiagonallyDominantConeTriangle`](@ref) + +## Target nodes + +`DiagonallyDominantBridge` creates: + + * `F` in `MOI.GreaterThan{T}` (the row dominance inequalities) + * `MOI.ScalarAffineFunction{T}` in `MOI.GreaterThan{T}` (the + ``t_{ij} \\ge \\pm Q_{ij}`` inequalities), when `F` is not already + `MOI.ScalarAffineFunction{T}` +""" struct DiagonallyDominantBridge{T,F,G} <: MOI.Bridges.Constraint.AbstractBridge # |Qij| variables abs_vars::Vector{MOI.VariableIndex} diff --git a/src/Bridges/Constraint/empty.jl b/src/Bridges/Constraint/empty.jl index 09b95409c..d6a0e68b2 100644 --- a/src/Bridges/Constraint/empty.jl +++ b/src/Bridges/Constraint/empty.jl @@ -1,3 +1,21 @@ +""" + EmptyBridge{T,F<:MOI.AbstractVectorFunction} <: Bridges.Constraint.AbstractBridge + +`EmptyBridge` is a no-op bridge for [`SumOfSquares.EmptyCone`](@ref). It is +used as the inner-most node when an SOS decomposition produces no gram +entries (e.g. when the basis is empty), so that the bridge graph still has +a target node to terminate at. + +## Source node + +`EmptyBridge` supports: + + * `F` in [`SumOfSquares.EmptyCone`](@ref) + +## Target nodes + +`EmptyBridge` creates no constraints and no variables. +""" struct EmptyBridge{T,F<:MOI.AbstractVectorFunction} <: MOI.Bridges.Constraint.AbstractBridge end diff --git a/src/Bridges/Constraint/psd2x2.jl b/src/Bridges/Constraint/psd2x2.jl index 3dee91b45..772917814 100644 --- a/src/Bridges/Constraint/psd2x2.jl +++ b/src/Bridges/Constraint/psd2x2.jl @@ -1,6 +1,31 @@ # PSD constraints on 2x2 matrices are SOC representable. # [Q11 Q12] is PSD iff Q11, Q22 ≥ 0 and Q11*Q22 ≥ Q12 ^2 # [Q12 Q22] <=> 2*Q11*Q22 ≥ (√2*Q12)^2 + +""" + PositiveSemidefinite2x2Bridge{T,F} <: Bridges.Constraint.AbstractBridge + +`PositiveSemidefinite2x2Bridge` implements a reformulation from +[`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref) into a rotated +second-order cone constraint. + +A symmetric ``2 \\times 2`` matrix is PSD iff +``Q_{11}, Q_{22} \\ge 0`` and ``Q_{11} Q_{22} \\ge Q_{12}^2``, which +rewrites as ``2 Q_{11} Q_{22} \\ge (\\sqrt{2} Q_{12})^2``, exactly the +definition of `MOI.RotatedSecondOrderCone(3)`. + +## Source node + +`PositiveSemidefinite2x2Bridge` supports: + + * `G` in [`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref) + +## Target nodes + +`PositiveSemidefinite2x2Bridge` creates: + + * `F` in `MOI.RotatedSecondOrderCone` (of dimension 3) +""" struct PositiveSemidefinite2x2Bridge{T,F} <: MOI.Bridges.Constraint.AbstractBridge rsoc::MOI.ConstraintIndex{F,MOI.RotatedSecondOrderCone} diff --git a/src/Bridges/Constraint/sos_polynomial.jl b/src/Bridges/Constraint/sos_polynomial.jl index 9fa0e6244..19c741954 100644 --- a/src/Bridges/Constraint/sos_polynomial.jl +++ b/src/Bridges/Constraint/sos_polynomial.jl @@ -1,3 +1,29 @@ +""" + SOSPolynomialBridge{T,F,DT,M,BT,B,G,CT,W} <: Bridges.Constraint.SetMapBridge + +`SOSPolynomialBridge` implements a reformulation from +[`SumOfSquares.SOSPolynomialSet{<:AbstractAlgebraicSet}`](@ref) into +[`SumOfSquares.WeightedSOSCone`](@ref). + +The polynomial coefficients are first reduced modulo the algebraic ideal +``\\mathcal{V}`` of the source domain, and the gram basis is computed by +the inner ideal `certificate`. The resulting decomposition +``p = b^\\top Q b`` is encoded as a `WeightedSOSCone` with a single gram +basis and unit weight (one for each sparsity block when the certificate +exploits sparsity, after `_flatten`). + +## Source node + +`SOSPolynomialBridge` supports: + + * `F` in [`SumOfSquares.SOSPolynomialSet{<:SemialgebraicSets.AbstractAlgebraicSet}`](@ref) + +## Target nodes + +`SOSPolynomialBridge` creates: + + * `F` in [`SumOfSquares.WeightedSOSCone`](@ref) +""" struct SOSPolynomialBridge{ T, F<:MOI.AbstractVectorFunction, diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 3144f0ac5..0965ec7f2 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -1,3 +1,33 @@ +""" + SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W} <: Bridges.Constraint.AbstractBridge + +`SOSPolynomialInSemialgebraicSetBridge` implements a reformulation from +[`SumOfSquares.SOSPolynomialSet{<:BasicSemialgebraicSet}`](@ref) into a +single [`SumOfSquares.WeightedSOSCone`](@ref) carrying the whole +Putinar-style decomposition +``p = \\sigma_0 + \\sum_i g_i \\sigma_i``. + +For each inequality ``g_i`` of the basic semialgebraic domain, the bridge +collects a multiplier basis ``b_i`` and produces a `WeightedSOSCone` with +gram bases ``[b_0, b_1, \\dots]`` and weights ``[1, g_1, \\dots]``. The +polynomial is reduced modulo the algebraic part ``\\mathcal{V}`` before +being passed on. Whichever bridge handles `WeightedSOSCone` downstream +(`Variable.KernelBridge`, `Variable.LowRankBridge` or `Constraint.ImageBridge`) +ultimately produces the PSD constraints; this bridge adds no variables of +its own. + +## Source node + +`SOSPolynomialInSemialgebraicSetBridge` supports: + + * `F` in [`SumOfSquares.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}`](@ref) + +## Target nodes + +`SOSPolynomialInSemialgebraicSetBridge` creates: + + * `F` in [`SumOfSquares.WeightedSOSCone`](@ref) +""" struct SOSPolynomialInSemialgebraicSetBridge{ T, F<:MOI.AbstractVectorFunction, diff --git a/src/Bridges/Variable/copositive_inner.jl b/src/Bridges/Variable/copositive_inner.jl index b73803450..c781bfbb6 100644 --- a/src/Bridges/Variable/copositive_inner.jl +++ b/src/Bridges/Variable/copositive_inner.jl @@ -1,3 +1,30 @@ +""" + CopositiveInnerBridge{T,S} <: Bridges.Variable.AbstractBridge + +`CopositiveInnerBridge` implements a reformulation from +[`SumOfSquares.CopositiveInnerCone`](@ref) into a PSD-like inner cone plus +a nonnegativity constraint on the strictly off-diagonal entries. + +The cone of inner-copositive matrices contains symmetric matrices that can +be written as the sum of a matrix in `set.psd_inner` (e.g. PSD) and a +matrix that is zero on the diagonal and has nonnegative off-diagonal +entries. The bridge introduces constrained variables for each summand. + +## Source node + +`CopositiveInnerBridge` supports: + + * [`SumOfSquares.CopositiveInnerCone{S}`](@ref) + +## Target nodes + +`CopositiveInnerBridge` creates: + + * `MOI.VectorOfVariables` in `S` (the inner PSD-like cone, e.g. + `MOI.PositiveSemidefiniteConeTriangle`) + * `MOI.VectorOfVariables` in `MOI.Nonnegatives` (one per strictly + upper-triangular entry) +""" struct CopositiveInnerBridge{T,S} <: MOI.Bridges.Variable.AbstractBridge matrix_variables::Vector{MOI.VariableIndex} matrix_constraint::MOI.ConstraintIndex{MOI.VectorOfVariables,S} diff --git a/src/Bridges/Variable/kernel.jl b/src/Bridges/Variable/kernel.jl index 58c6aa93d..7c602e0ef 100644 --- a/src/Bridges/Variable/kernel.jl +++ b/src/Bridges/Variable/kernel.jl @@ -1,3 +1,42 @@ +""" + KernelBridge{T,M} <: Bridges.Variable.AbstractBridge + +`KernelBridge` implements a reformulation from +[`SumOfSquares.WeightedSOSCone`](@ref) into one PSD-style constrained +variable per gram basis. + +For each gram basis ``b_i`` and weight ``w_i`` of the source cone, +`KernelBridge` adds the corresponding gram matrix ``Q_i`` as constrained +variables in `matrix_cone(M, length(b_i))` (which downstream bridges may +unfold further, e.g. into `MOI.PositiveSemidefiniteConeTriangle`). It then +sets the bridged polynomial-coefficient functions to the affine combination +``\\sum_i w_i \\cdot b_i^\\top Q_i b_i`` expressed in the source basis. + +This is one of the two main back-ends for `WeightedSOSCone`. The other is +[`SumOfSquares.Bridges.Constraint.ImageBridge`](@ref). `KernelBridge` is +preferred for solvers that natively support PSD as constrained variables +(e.g. Mosek's `barvar` matrix variables), while `ImageBridge` is preferred +for solvers that only support PSD as a constraint (e.g. Clarabel). The +choice is driven by `MOI.Bridges.bridging_cost`. + +`KernelBridge` does not handle [`MultivariateBases.LagrangeBasis`](@ref) — +those bases dispatch to [`SumOfSquares.Bridges.Variable.LowRankBridge`](@ref) +instead. + +## Source node + +`KernelBridge` supports: + + * [`SumOfSquares.WeightedSOSCone{M,B}`](@ref) when + `!(B <: MB.LagrangeBasis)` + +## Target nodes + +`KernelBridge` creates: + + * `MOI.VectorOfVariables` in `SOS.matrix_cone(M, length(b_i))` for each + gram basis ``b_i`` of the cone +""" struct KernelBridge{T,M} <: MOI.Bridges.Variable.AbstractBridge affine::Vector{MOI.ScalarAffineFunction{T}} variables::Vector{Vector{MOI.VariableIndex}} diff --git a/src/Bridges/Variable/lowrank.jl b/src/Bridges/Variable/lowrank.jl index c55b85dcd..d1cca4723 100644 --- a/src/Bridges/Variable/lowrank.jl +++ b/src/Bridges/Variable/lowrank.jl @@ -1,3 +1,42 @@ +""" + LowRankBridge{T,M} <: Bridges.Variable.AbstractBridge + +`LowRankBridge` implements a reformulation from +[`SumOfSquares.WeightedSOSCone`](@ref) with a Lagrange basis into +`LowRankOpt.SetDotProducts` over a PSD-style cone. + +When the gram basis is sampled at a fixed set of points (a Lagrange basis), +the polynomial-coefficient ⇄ gram-matrix linear map factorises through a +list of rank-1 outer products, one per Lagrange node. This bridge exposes +that structure by adding constrained variables in +`LowRankOpt.SetDotProducts{WITHOUT_SET}` over +`SOS.matrix_cone(M, length(b_i))`, with the `LowRankOpt.Factorization` of +each node carrying the basis-transformation row and the weight value. +Solvers that understand `LowRankOpt.SetDotProducts` (e.g. +[Loraine.jl](https://github.com/kocvara/Loraine.jl) and +[BurerMonteiro](https://github.com/blegat/LowRankOpt.jl) backends of +`LowRankOpt`) can then exploit the rank-1 structure directly; otherwise, +the [`LowRankOpt` bridges](https://github.com/blegat/LowRankOpt.jl) unfold +the cone into a plain PSD constraint. + +For non-Lagrange bases, the structure is no longer rank-1 and the bridge +falls back to [`SumOfSquares.Bridges.Variable.KernelBridge`](@ref). + +## Source node + +`LowRankBridge` supports: + + * [`SumOfSquares.WeightedSOSCone{M,B}`](@ref) when + `B <: MB.LagrangeBasis` + +## Target nodes + +`LowRankBridge` creates: + + * `MOI.VectorOfVariables` in + `LowRankOpt.SetDotProducts{WITHOUT_SET, S, LowRankOpt.TriangleVectorization{T,LowRankOpt.Factorization{T,Vector{T},Array{T,0}}}}` + wrapping `SOS.matrix_cone(M, length(b_i))` for each gram basis ``b_i`` +""" struct LowRankBridge{T,M} <: MOI.Bridges.Variable.AbstractBridge affine::Vector{MOI.ScalarAffineFunction{T}} variables::Vector{Vector{MOI.VariableIndex}} diff --git a/src/Bridges/Variable/psd2x2.jl b/src/Bridges/Variable/psd2x2.jl index 6e62f111d..c963f7ecc 100644 --- a/src/Bridges/Variable/psd2x2.jl +++ b/src/Bridges/Variable/psd2x2.jl @@ -1,3 +1,25 @@ +""" + PositiveSemidefinite2x2Bridge{T} <: Bridges.Variable.AbstractBridge + +`PositiveSemidefinite2x2Bridge` is the variable counterpart of +[`SumOfSquares.Bridges.Constraint.PositiveSemidefinite2x2Bridge`](@ref): it +reformulates constrained variables in +[`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref) into +constrained variables in `MOI.RotatedSecondOrderCone(3)`, using +``Q_{11} Q_{22} \\ge Q_{12}^2 \\Leftrightarrow 2 Q_{11} Q_{22} \\ge (\\sqrt 2 Q_{12})^2``. + +## Source node + +`PositiveSemidefinite2x2Bridge` supports: + + * [`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref) + +## Target nodes + +`PositiveSemidefinite2x2Bridge` creates: + + * `MOI.VectorOfVariables` in `MOI.RotatedSecondOrderCone` (of dimension 3) +""" struct PositiveSemidefinite2x2Bridge{T} <: MOI.Bridges.Variable.AbstractBridge variables::Vector{MOI.VariableIndex} rsoc::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.RotatedSecondOrderCone} diff --git a/src/Bridges/Variable/scaled_diagonally_dominant.jl b/src/Bridges/Variable/scaled_diagonally_dominant.jl index 55ca58210..0dc3017a1 100644 --- a/src/Bridges/Variable/scaled_diagonally_dominant.jl +++ b/src/Bridges/Variable/scaled_diagonally_dominant.jl @@ -1,15 +1,29 @@ """ - struct ScaledDiagonallyDominantBridge{T} <: MOI.Bridges.Variable.AbstractBridge - side_dimension::Int - variables::Vector{Vector{MOI.VariableIndex}} - constraints::Vector{MOI.ConstraintIndex{ - MOI.VectorOfVariables, SOS.PositiveSemidefinite2x2ConeTriangle}} - end + ScaledDiagonallyDominantBridge{T} <: Bridges.Variable.AbstractBridge + +`ScaledDiagonallyDominantBridge` implements a reformulation from +[`SumOfSquares.ScaledDiagonallyDominantConeTriangle`](@ref) into a sum of +``2 \\times 2`` PSD blocks. + +A matrix is scaled diagonally dominant iff it is the sum of PSD matrices +``M^{(i,j)}`` that are zero except on entries ``(i, i)``, ``(i, j)`` and +``(j, j)`` [Ahmadi2017; Lemma 9](@cite). Each such block ``M^{(i,j)}`` is +encoded by adding constrained variables in +[`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref). + +## Source node + +`ScaledDiagonallyDominantBridge` supports: + + * [`SumOfSquares.ScaledDiagonallyDominantConeTriangle`](@ref) + +## Target nodes + +`ScaledDiagonallyDominantBridge` creates: -A matrix is SDD iff it is the sum of psd matrices Mij that are zero except -for entries ii, ij and jj [Ahmadi2017; Lemma 9](@cite). This bridge substitute the -constrained variables in [`SOS.ScaledDiagonallyDominantConeTriangle`](@ref) -into a sum of constrained variables in [`SOS.PositiveSemidefinite2x2ConeTriangle`](@ref). + * `MOI.VectorOfVariables` in + [`SumOfSquares.PositiveSemidefinite2x2ConeTriangle`](@ref), one per + strictly upper-triangular position ``(i, j)`` """ struct ScaledDiagonallyDominantBridge{T} <: MOI.Bridges.Variable.AbstractBridge side_dimension::Int From 44475644cad85fa270d3fdbb58f5013c010af3ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 10:07:07 +0200 Subject: [PATCH 05/17] test and example for augmented path --- src/Bridges/Constraint/image.jl | 40 ++++++++++++++++++++ test/Bridges/Constraint/image.jl | 65 ++++++++++++++++++++++++++++++++ 2 files changed, 105 insertions(+) diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index 21eff4cc8..5026bff05 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -119,6 +119,46 @@ R = \\begin{bmatrix}-\\lambda\\end{bmatrix}. The slack `λ` is constrained only by `Q ⪰ 0` and `R ⪰ 0`, i.e. by the PSD constraints on the gram matrices. +### Sub-optimal greedy: slacks that could be avoided + +The order in which the contributions of one entry are processed matters. +The bridge walks each entry's weight monomials in their stored order and +commits to either anchoring the entry to the first unanchored target it +sees, or to a slack variable if that first target is already anchored. +That decision can introduce a slack variable that a different order would +have avoided. + +Consider weight `w = 1 + x` with gram basis `[1, x, x^2]` and basis +`{1, x, x^2, x^3, x^4, x^5}`. The ``6`` PSD entries are +```text +k = 1: (1, 1) contributes to 1 and x +k = 2: (1, x) contributes to x and x^2 +k = 3: (x, x) contributes to x^2 and x^3 +k = 4: (1, x^2) contributes to x^2 and x^3 +k = 5: (x, x^2) contributes to x^3 and x^4 +k = 6: (x^2, x^2) contributes to x^4 and x^5 +``` +The greedy walks weight monomials in the order `[1, x]`, so for entry +``k = 4`` it sees the contribution to `x^2` first. By then `x^2` is +already anchored at ``k = 3``, so the bridge introduces a slack `λ`, +records the second contribution (`x^3`, currently unanchored) into +`pending`, and the trail of consequences cascades into 5 anchors and 1 +slack — with `x^5` left to a zero constraint. + +A bipartite matching between gram entries (left) and basis monomials +(right) — an edge `(k, t)` for each weight-monomial contribution from +``k`` to a target ``t`` in the basis — admits a perfect matching of size +6 here: +``1 \\to 1,\\ 2 \\to x,\\ 3 \\to x^2,\\ 4 \\to x^3,\\ 5 \\to x^4,\\ 6 \\to x^5``. +Anchoring along that matching would use 6 anchors and 0 slack. Even just +swapping the two contributions of entry ``k = 4`` (and similarly for +``k = 5`` and ``k = 6``) — i.e. preferring an unanchored target whenever +one is available, and falling back to subtracting from the existing +anchor for the other contribution — recovers the same outcome. The +underlying problem is unweighted bipartite max-cardinality matching, for +which Hopcroft–Karp or simple augmenting paths gives the optimum; the +weighted Hungarian assignment is not needed. + ## Source node `ImageBridge` supports: diff --git a/test/Bridges/Constraint/image.jl b/test/Bridges/Constraint/image.jl index 7e090aa3b..7a7f3c648 100644 --- a/test/Bridges/Constraint/image.jl +++ b/test/Bridges/Constraint/image.jl @@ -163,6 +163,71 @@ function test_runtests_polynomial_weight() return end +# Same polynomial weight as `test_runtests_polynomial_weight` but with a +# 3-element gram basis `[1, x, x^2]`. Here the greedy contribution order +# inside an entry is sub-optimal: for entry `(1, 3)` the algorithm sees +# `x^2` first (already anchored by entry `(2, 2)`) and so introduces a slack +# `λ`, even though processing `(1, 3)`'s second contribution (`x^3`, +# currently unanchored) first would have allowed it to anchor a fresh +# monomial and use the `entry_anchored = true` subtract branch for `x^2`. +# Bipartite max-cardinality matching between gram entries and basis +# monomials would yield 6 anchors and zero slacks here; the current greedy +# yields 5 anchors and 1 slack, with `x^5` becoming a zero constraint. +function test_runtests_polynomial_weight_avoidable_slack() + T = Float64 + @polyvar x + MOI.Bridges.runtests( + SumOfSquares.Bridges.Constraint.ImageBridge, + model -> begin + MOI.add_constraint( + model, + MOI.VectorAffineFunction{T}( + MOI.VectorAffineTerm{T}[], + T[1, 3, 5, 5, 3, 1], + ), + SumOfSquares.WeightedSOSCone{ + MOI.PositiveSemidefiniteConeTriangle, + }( + MB.SubBasis{MB.Monomial}([ + x^0, + x, + x^2, + x^3, + x^4, + x^5, + ]), + [MB.SubBasis{MB.Monomial}([x^0, x, x^2])], + [MB.algebra_element(T(1) * x^0 + T(1) * x)], + ), + ) + end, + model -> begin + λ = MOI.add_variable(model) + MOI.add_constraint( + model, + MOI.Utilities.operate( + vcat, + T, + T(1), + T(1), + T(3) + T(2) * λ, + T(-1) * λ, + T(1), + T(1), + ), + MOI.PositiveSemidefiniteConeTriangle(3), + ) + MOI.add_constraint( + model, + MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[0]), + MOI.Zeros(1), + ) + end; + cannot_unbridge = true, + ) + return +end + # Regression test for an oversight in # `MOI.Bridges.added_constraint_types(ImageBridge)`: the bridge can produce a # constraint in any of the four sets returned by `SOS.matrix_cone(M, k)` for From f8fe7914c528f24d2e37efb3959a20037b900d12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 10:11:24 +0200 Subject: [PATCH 06/17] Fix format --- .../sos_polynomial_in_semialgebraic_set.jl | 37 +++++++++---------- test/Bridges/Constraint/image.jl | 9 +---- 2 files changed, 19 insertions(+), 27 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 0965ec7f2..c5efa15f4 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -83,12 +83,8 @@ function MOI.Bridges.Constraint.bridge_constraint( some_mono = first(MB.keys_as_monomials(set.basis)) unit_poly = MP.polynomial(MP.term(one(T), MP.constant_monomial(some_mono))) gram_bases_raw = Any[sigma_0_basis] - weights_raw = W[ - MB.algebra_element( - MB.sparse_coefficients(unit_poly), - implicit_basis, - ), - ] + weights_raw = + W[MB.algebra_element(MB.sparse_coefficients(unit_poly), implicit_basis),] for index in Certificate.preorder_indices(set.certificate, preprocessed) push!( gram_bases_raw, @@ -108,14 +104,9 @@ function MOI.Bridges.Constraint.bridge_constraint( # turn the latter into a uniform `Vector{GB}` (replicating the matching # weight for each sparsity block). The σ_0 block is the first entry, so # `sigma_0_indices` is what `_flatten` reports for that single entry. - sigma_0_gram_bases, _, sigma_0_indices = _flatten( - identity.([sigma_0_basis]), - [weights_raw[1]], - ) - gram_bases, weights, _ = _flatten( - identity.(gram_bases_raw), - weights_raw, - ) + sigma_0_gram_bases, _, sigma_0_indices = + _flatten(identity.([sigma_0_basis]), [weights_raw[1]]) + gram_bases, weights, _ = _flatten(identity.(gram_bases_raw), weights_raw) new_basis = Certificate.zero_basis( ideal_cert, MB.explicit_basis(poly_reduced), @@ -198,7 +189,10 @@ function MOI.Bridges.Constraint.concrete_bridge_type( end # A small helper to keep the `concrete_bridge_type` readable. -function SOSPolynomialBridgeType_helper_gram_basis_type(::Type{IC}, ::Type) where {IC} +function SOSPolynomialBridgeType_helper_gram_basis_type( + ::Type{IC}, + ::Type, +) where {IC} return _eltype(MA.promote_operation(SOS.Certificate.gram_basis, IC)) end @@ -272,7 +266,11 @@ function MOI.get( # stored on the `WeightedSOSCone` set. set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint) return MultivariateMoments.moment_vector( - MOI.get(model, MOI.ConstraintDual(attr.result_index), bridge.constraint), + MOI.get( + model, + MOI.ConstraintDual(attr.result_index), + bridge.constraint, + ), set.basis, ) end @@ -306,9 +304,10 @@ function MOI.get( # are σ_1, σ_2, … (one Lagrangian multiplier per inequality). set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint) n_total = length(set.gram_bases) - sigma_0_max = bridge.sigma_0_indices isa Int ? bridge.sigma_0_indices : - last(bridge.sigma_0_indices) - return map((sigma_0_max + 1):n_total) do i + sigma_0_max = + bridge.sigma_0_indices isa Int ? bridge.sigma_0_indices : + last(bridge.sigma_0_indices) + return map((sigma_0_max+1):n_total) do i return MOI.get( model, SOS.GramMatrixAttribute(; diff --git a/test/Bridges/Constraint/image.jl b/test/Bridges/Constraint/image.jl index 7a7f3c648..952f59279 100644 --- a/test/Bridges/Constraint/image.jl +++ b/test/Bridges/Constraint/image.jl @@ -188,14 +188,7 @@ function test_runtests_polynomial_weight_avoidable_slack() SumOfSquares.WeightedSOSCone{ MOI.PositiveSemidefiniteConeTriangle, }( - MB.SubBasis{MB.Monomial}([ - x^0, - x, - x^2, - x^3, - x^4, - x^5, - ]), + MB.SubBasis{MB.Monomial}([x^0, x, x^2, x^3, x^4, x^5]), [MB.SubBasis{MB.Monomial}([x^0, x, x^2])], [MB.algebra_element(T(1) * x^0 + T(1) * x)], ), From d671fef4e4843564250492e93ad8fcea2d271c1d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 11:21:53 +0200 Subject: [PATCH 07/17] fixes --- .../sos_polynomial_in_semialgebraic_set.jl | 79 ++++++++++++------- 1 file changed, 50 insertions(+), 29 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index c5efa15f4..23ba16fdc 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -41,11 +41,11 @@ struct SOSPolynomialInSemialgebraicSetBridge{ } <: MOI.Bridges.Constraint.AbstractBridge constraint::MOI.ConstraintIndex{F,SOS.WeightedSOSCone{M,NB,GB,W}} basis::BT - # `sigma_0_indices` records the range of `gram_bases` that correspond to - # σ_0 (possibly multiple if the certificate uses sparsity). The - # remaining gram bases (after this range) are the Lagrangian multipliers - # σ_1, σ_2, … in order. - sigma_0_indices::Union{Int,UnitRange{Int}} + # `multiplier_indices[1]` is the range of `gram_bases` that corresponds + # to σ_0 (possibly multiple if the certificate uses sparsity). + # `multiplier_indices[1+i]` is the range of `gram_bases` of the + # Lagrangian multiplier σ_i. + multiplier_indices::Vector{Union{Int,UnitRange{Int}}} end function MOI.Bridges.Constraint.bridge_constraint( @@ -82,9 +82,13 @@ function MOI.Bridges.Constraint.bridge_constraint( # concrete element type and fit in a homogeneous `Vector{W}`. some_mono = first(MB.keys_as_monomials(set.basis)) unit_poly = MP.polynomial(MP.term(one(T), MP.constant_monomial(some_mono))) - gram_bases_raw = Any[sigma_0_basis] - weights_raw = - W[MB.algebra_element(MB.sparse_coefficients(unit_poly), implicit_basis),] + # Append the Lagrangian multiplier bases (σ_1, σ_2, …) first and σ_0 + # last so that, downstream, `Variable.KernelBridge` allocates the σ_i + # gram variables before the σ_0 ones — matching the variable ordering + # of the original two-bridge pipeline (and therefore the order assumed + # by the cached `Mock/` tests). + gram_bases_raw = Any[] + weights_raw = W[] for index in Certificate.preorder_indices(set.certificate, preprocessed) push!( gram_bases_raw, @@ -99,13 +103,35 @@ function MOI.Bridges.Constraint.bridge_constraint( ), ) end + push!(gram_bases_raw, sigma_0_basis) + push!( + weights_raw, + MB.algebra_element(MB.sparse_coefficients(unit_poly), implicit_basis), + ) # `gram_basis` returns either a single basis or a `Vector{basis}` when # sparsity is in play. Reuse `_flatten` from the SOSPolynomial bridge to # turn the latter into a uniform `Vector{GB}` (replicating the matching - # weight for each sparsity block). The σ_0 block is the first entry, so - # `sigma_0_indices` is what `_flatten` reports for that single entry. - sigma_0_gram_bases, _, sigma_0_indices = - _flatten(identity.([sigma_0_basis]), [weights_raw[1]]) + # weight for each sparsity block). At the same time, record where each + # multiplier σ_i lives in the flattened `gram_bases` so that we can + # rebuild block-diagonal gram and moment matrices per multiplier. + # `multiplier_indices` is `[σ_0, σ_1, …]` in that order, even though the + # underlying `gram_bases` lay them out as `[σ_1, …, σ_0]`. + raw_ranges = Union{Int,UnitRange{Int}}[] + cur = 0 + for raw in gram_bases_raw + n = raw isa AbstractVector ? length(raw) : 1 + if n == 1 + cur += 1 + push!(raw_ranges, cur) + else + push!(raw_ranges, (cur+1):(cur+n)) + cur += n + end + end + # σ_0 is the last entry of `gram_bases_raw`; bring it back to position 1 + # in `multiplier_indices` without splatting the `UnitRange` element. + multiplier_indices = Union{Int,UnitRange{Int}}[raw_ranges[end]] + append!(multiplier_indices, @view raw_ranges[1:(end-1)]) gram_bases, weights, _ = _flatten(identity.(gram_bases_raw), weights_raw) new_basis = Certificate.zero_basis( ideal_cert, @@ -123,7 +149,7 @@ function MOI.Bridges.Constraint.bridge_constraint( return SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}( constraint, set.basis, - sigma_0_indices, + multiplier_indices, ) end @@ -292,7 +318,7 @@ function MOI.get( # σ_0 block(s) occupy in the underlying `WeightedSOSCone`'s # `gram_bases`. SOS.check_multiplier_index_bounds(attr, 0:0) - return _get(model, attr, bridge.constraint, bridge.sigma_0_indices) + return _get(model, attr, bridge.constraint, bridge.multiplier_indices[1]) end function MOI.get( @@ -300,21 +326,16 @@ function MOI.get( attr::SOS.LagrangianMultipliers, bridge::SOSPolynomialInSemialgebraicSetBridge, ) - # σ_0 occupies the `sigma_0_indices` slot(s); the remaining gram bases - # are σ_1, σ_2, … (one Lagrangian multiplier per inequality). - set = MOI.get(model, MOI.ConstraintSet(), bridge.constraint) - n_total = length(set.gram_bases) - sigma_0_max = - bridge.sigma_0_indices isa Int ? bridge.sigma_0_indices : - last(bridge.sigma_0_indices) - return map((sigma_0_max+1):n_total) do i - return MOI.get( - model, - SOS.GramMatrixAttribute(; - multiplier_index = i, - result_index = attr.result_index, - ), - bridge.constraint, + # `multiplier_indices[1]` is σ_0; the rest are σ_1, σ_2, … in order. + # Each multiplier is returned as a `BlockDiagonalGramMatrix`, even when + # its `gram_bases` range has a single entry, so that the result type is + # uniform across multipliers. + return map(2:length(bridge.multiplier_indices)) do i + idx = bridge.multiplier_indices[i] + range = idx isa Int ? (idx:idx) : idx + gram_attr = SOS.GramMatrixAttribute(; + result_index = attr.result_index, ) + return _get(model, gram_attr, bridge.constraint, range) end end From 96f1773bb007784e77cbabc26fe14237f6e9c2dc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 11:29:39 +0200 Subject: [PATCH 08/17] fix format --- src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 23ba16fdc..603570119 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -333,9 +333,7 @@ function MOI.get( return map(2:length(bridge.multiplier_indices)) do i idx = bridge.multiplier_indices[i] range = idx isa Int ? (idx:idx) : idx - gram_attr = SOS.GramMatrixAttribute(; - result_index = attr.result_index, - ) + gram_attr = SOS.GramMatrixAttribute(; result_index = attr.result_index) return _get(model, gram_attr, bridge.constraint, range) end end From 1ce98e817a335e9e9c458ad1b37d86a68ef67de3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 14:29:37 +0200 Subject: [PATCH 09/17] Fixes --- .../sos_polynomial_in_semialgebraic_set.jl | 30 ++++--- test/Bridges/lazy.jl | 88 +++++++++++++++++++ 2 files changed, 107 insertions(+), 11 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 603570119..cb24e1097 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -119,13 +119,18 @@ function MOI.Bridges.Constraint.bridge_constraint( raw_ranges = Union{Int,UnitRange{Int}}[] cur = 0 for raw in gram_bases_raw - n = raw isa AbstractVector ? length(raw) : 1 - if n == 1 - cur += 1 - push!(raw_ranges, cur) - else + # Preserve whether the multiplier's basis was a `Vector` (sparsity is + # in play, even with a single block) or a single basis: that controls + # whether the corresponding gram/moment attribute comes back as a + # `BlockDiagonalGramMatrix` (`UnitRange` index) or a plain + # `GramMatrix` (single `Int` index). + if raw isa AbstractVector + n = length(raw) push!(raw_ranges, (cur+1):(cur+n)) cur += n + else + cur += 1 + push!(raw_ranges, cur) end end # σ_0 is the last entry of `gram_bases_raw`; bring it back to position 1 @@ -327,13 +332,16 @@ function MOI.get( bridge::SOSPolynomialInSemialgebraicSetBridge, ) # `multiplier_indices[1]` is σ_0; the rest are σ_1, σ_2, … in order. - # Each multiplier is returned as a `BlockDiagonalGramMatrix`, even when - # its `gram_bases` range has a single entry, so that the result type is - # uniform across multipliers. + # `_get` returns a `BlockDiagonalGramMatrix` when given a `UnitRange` + # and a plain `GramMatrix` when given an `Int`, mirroring the original + # bridge's behaviour: sparsity → block-diagonal, no sparsity → plain. return map(2:length(bridge.multiplier_indices)) do i - idx = bridge.multiplier_indices[i] - range = idx isa Int ? (idx:idx) : idx gram_attr = SOS.GramMatrixAttribute(; result_index = attr.result_index) - return _get(model, gram_attr, bridge.constraint, range) + return _get( + model, + gram_attr, + bridge.constraint, + bridge.multiplier_indices[i], + ) end end diff --git a/test/Bridges/lazy.jl b/test/Bridges/lazy.jl index 9acc59cae..3bccb5a96 100644 --- a/test/Bridges/lazy.jl +++ b/test/Bridges/lazy.jl @@ -104,6 +104,94 @@ function test_dual_clarabel_uses_kernel_bridge() return end +# With a `domain` keyword, the constraint enters the bridge graph as a +# `SOSPolynomialSet{<:BasicSemialgebraicSet}`, and is first reformulated by +# `SOSPolynomialInSemialgebraicSetBridge` into a single `WeightedSOSCone`. +# From there, Clarabel — which natively supports +# `MOI.PositiveSemidefiniteConeTriangle` as a constraint — routes through +# `VectorSlackBridge` and `Variable.KernelBridge` rather than through +# `Constraint.ImageBridge`, because the bump in `ImageBridge`'s cost shifts +# the balance once the immediate PSD lands on the variable side. +function test_clarabel_with_domain_uses_kernel_bridge() + T = Float64 + @polyvar x y + model = Model(Clarabel.Optimizer) + set_silent(model) + @variable(model, α) + @objective(model, Max, α) + S = @set x >= 0 + cref = @constraint( + model, + 10 - x^2 - α * y in SOSCone(), + domain = S, + maxdegree = 4, + ) + optimize!(model) + backend = JuMP.backend(model) + # Walk down the chain: SOSPolynomialSet → InSemialgebraicSetBridge → + # WeightedSOSCone (via VectorSlackBridge) → KernelBridge. + lbo = backend.optimizer + cis = MOI.ConstraintIndex[] + for (F, S) in MOI.get(lbo, MOI.ListOfConstraintTypesPresent()) + if S <: SumOfSquares.SOSPolynomialSet + append!(cis, MOI.get(lbo, MOI.ListOfConstraintIndices{F,S}())) + end + end + @test length(cis) == 1 + outer = MOI.Bridges.bridge(lbo, only(cis)) + @test outer isa + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge + slack = MOI.Bridges.bridge(backend.optimizer, outer.constraint) + @test slack isa MOI.Bridges.Constraint.VectorSlackBridge + inner = MOI.Bridges.bridge(backend.optimizer, slack.slack_in_set) + @test inner isa SumOfSquares.Bridges.Variable.KernelBridge + return +end + +# With `Dualization.dual_optimizer`, the role of +# `PositiveSemidefiniteConeTriangle` flips between primal/dual. The chain +# is otherwise the same as `test_clarabel_with_domain_uses_kernel_bridge`: +# `SOSPolynomialInSemialgebraicSetBridge` → `VectorSlackBridge` → +# `Variable.KernelBridge`. `Constraint.ImageBridge` is not added to the +# graph at all here, because the dualization wrapper hides the +# `PositiveSemidefiniteConeTriangle` one level deeper than the cost +# heuristic in `image.jl` looks (a follow-up could tune the bridging cost +# so that `ImageBridge` wins on the dual side, like it does for +# `Variable.KernelBridge` on the primal side). +function test_dual_clarabel_with_domain_uses_kernel_bridge() + T = Float64 + @polyvar x y + model = Model(Dualization.dual_optimizer(Clarabel.Optimizer)) + set_silent(model) + @variable(model, α) + @objective(model, Max, α) + S = @set x >= 0 + cref = @constraint( + model, + 10 - x^2 - α * y in SOSCone(), + domain = S, + maxdegree = 4, + ) + optimize!(model) + backend = JuMP.backend(model) + lbo = backend.optimizer + cis = MOI.ConstraintIndex[] + for (F, S) in MOI.get(lbo, MOI.ListOfConstraintTypesPresent()) + if S <: SumOfSquares.SOSPolynomialSet + append!(cis, MOI.get(lbo, MOI.ListOfConstraintIndices{F,S}())) + end + end + @test length(cis) == 1 + outer = MOI.Bridges.bridge(lbo, only(cis)) + @test outer isa + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge + slack = MOI.Bridges.bridge(lbo, outer.constraint) + @test slack isa MOI.Bridges.Constraint.VectorSlackBridge + inner = MOI.Bridges.bridge(lbo, slack.slack_in_set) + @test inner isa SumOfSquares.Bridges.Variable.KernelBridge + return +end + end TestLazy.runtests() From 6dbc16f82af1a424733cbab7039b7e7bb3dc3527 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 15:18:58 +0200 Subject: [PATCH 10/17] Fix --- .../sos_polynomial_in_semialgebraic_set.jl | 68 ++++++++++++++----- 1 file changed, 51 insertions(+), 17 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index cb24e1097..415f5b254 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -68,13 +68,7 @@ function MOI.Bridges.Constraint.bridge_constraint( # domain so the rest of the certificate operates modulo it. domain_V = MP.similar(set.domain.V, T) poly_reduced = Certificate.reduced_polynomial(ideal_cert, poly, domain_V) - sigma_0_basis = Certificate.gram_basis( - ideal_cert, - Certificate.with_variables(poly_reduced, set.domain), - ) implicit_basis = MB.implicit_basis(set.basis) - # Build `[1, g_1, …]` weights and `[σ_0_basis, σ_1_basis, …]` gram bases - # so the whole Putinar decomposition becomes a single `WeightedSOSCone`. preprocessed = Certificate.preprocessed_domain(set.certificate, set.domain, poly) # Use the same `sparse_coefficients`-based construction for the unit @@ -82,23 +76,63 @@ function MOI.Bridges.Constraint.bridge_constraint( # concrete element type and fit in a homogeneous `Vector{W}`. some_mono = first(MB.keys_as_monomials(set.basis)) unit_poly = MP.polynomial(MP.term(one(T), MP.constant_monomial(some_mono))) - # Append the Lagrangian multiplier bases (σ_1, σ_2, …) first and σ_0 - # last so that, downstream, `Variable.KernelBridge` allocates the σ_i - # gram variables before the σ_0 ones — matching the variable ordering - # of the original two-bridge pipeline (and therefore the order assumed - # by the cached `Mock/` tests). - gram_bases_raw = Any[] - weights_raw = W[] + # Collect the multiplier bases (σ_1, σ_2, …) and the inequality + # polynomials g_i. + multiplier_bases = Any[] + g_polynomials = [] for index in Certificate.preorder_indices(set.certificate, preprocessed) push!( - gram_bases_raw, + multiplier_bases, Certificate.multiplier_basis(set.certificate, index, preprocessed), ) - g = Certificate.generator(set.certificate, index, preprocessed) + push!( + g_polynomials, + Certificate.generator(set.certificate, index, preprocessed), + ) + end + # σ_0's gram basis must absorb both `p` and the monomials produced by + # `g_i * b_i_j * b_i_k`. We mirror the original two-bridge pipeline by + # symbolically augmenting `poly_reduced` with the unit-coefficient + # versions of those products before asking the ideal certificate for + # the gram basis — otherwise the Newton polytope would be that of `p` + # alone, which is generally too small to host an SOS decomposition. + augmented = SOS.MA.copy(poly_reduced) + for (mb, g_i) in zip(multiplier_bases, g_polynomials) + g_alg = MB.algebra_element( + MB.sparse_coefficients(one(T) * similar(g_i, T)), + implicit_basis, + ) + block_iter = mb isa AbstractVector ? mb : (mb,) + for block in block_iter + for b1 in block, b2 in block + term = MB.algebra_element(SA.star(b1) * b2) + MA.operate!( + SA.UnsafeAddMul(*), + augmented, + term, + g_alg, + ) + end + end + end + MA.operate!(SA.canonical, SA.coeffs(augmented)) + sigma_0_basis = Certificate.gram_basis( + ideal_cert, + Certificate.with_variables(augmented, set.domain), + ) + # Build `[g_1, …, 1]` weights and `[σ_1_basis, …, σ_0_basis]` gram + # bases. σ_0 goes last so that downstream `Variable.KernelBridge` + # allocates the σ_i gram variables before the σ_0 ones — matching the + # variable ordering of the original two-bridge pipeline (and therefore + # the order assumed by the cached `Mock/` tests). + gram_bases_raw = Any[] + weights_raw = W[] + for (mb, g_i) in zip(multiplier_bases, g_polynomials) + push!(gram_bases_raw, mb) push!( weights_raw, MB.algebra_element( - MB.sparse_coefficients(one(T) * similar(g, T)), + MB.sparse_coefficients(one(T) * similar(g_i, T)), implicit_basis, ), ) @@ -148,7 +182,7 @@ function MOI.Bridges.Constraint.bridge_constraint( new_coeffs = SA.coeffs(poly_reduced, new_basis) constraint = MOI.add_constraint( model, - MOI.Utilities.vectorize(SA.values(new_coeffs)), + MOI.Utilities.vectorize(new_coeffs), SOS.WeightedSOSCone{M}(new_basis, gram_bases, weights), ) return SOSPolynomialInSemialgebraicSetBridge{T,F,DT,CT,BT,M,NB,GB,W}( From 225a887b6c99d21001d6d565eb2246b9659c4498 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 15:28:29 +0200 Subject: [PATCH 11/17] Fix format --- .../Constraint/sos_polynomial_in_semialgebraic_set.jl | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 415f5b254..276cc1865 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -106,12 +106,7 @@ function MOI.Bridges.Constraint.bridge_constraint( for block in block_iter for b1 in block, b2 in block term = MB.algebra_element(SA.star(b1) * b2) - MA.operate!( - SA.UnsafeAddMul(*), - augmented, - term, - g_alg, - ) + MA.operate!(SA.UnsafeAddMul(*), augmented, term, g_alg) end end end From 4fb02a20d31564979b518f25ed418b69c3799189 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 16:03:55 +0200 Subject: [PATCH 12/17] Fixes --- docs/src/bridges.md | 2 +- docs/src/reference/standard_form.md | 2 +- src/Bridges/Variable/copositive_inner.jl | 4 ++-- src/Bridges/Variable/kernel.jl | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/bridges.md b/docs/src/bridges.md index a18229ada..8b52ddea7 100644 --- a/docs/src/bridges.md +++ b/docs/src/bridges.md @@ -116,7 +116,7 @@ natively — the bridge graph is what makes this transparent. ## How LowRankOpt exploits rank-1 structure -When the gram basis is a [`MultivariateBases.LagrangeBasis`](@ref) — the +When the gram basis is a `MultivariateBases.LagrangeBasis` — the polynomial is sampled at a fixed set of points — the relation between gram-matrix entries and polynomial coefficients factorises through a list of rank-1 outer products, one per Lagrange node: diff --git a/docs/src/reference/standard_form.md b/docs/src/reference/standard_form.md index 2d7cd6dbf..2bd127728 100644 --- a/docs/src/reference/standard_form.md +++ b/docs/src/reference/standard_form.md @@ -62,4 +62,4 @@ SumOfSquares.ScaledDiagonallyDominantConeTriangle ## Bridges -See the dedicated [Bridges](@ref) reference page. +See the dedicated bridges page in the API Reference section. diff --git a/src/Bridges/Variable/copositive_inner.jl b/src/Bridges/Variable/copositive_inner.jl index c781bfbb6..68685bae1 100644 --- a/src/Bridges/Variable/copositive_inner.jl +++ b/src/Bridges/Variable/copositive_inner.jl @@ -2,7 +2,7 @@ CopositiveInnerBridge{T,S} <: Bridges.Variable.AbstractBridge `CopositiveInnerBridge` implements a reformulation from -[`SumOfSquares.CopositiveInnerCone`](@ref) into a PSD-like inner cone plus +`SumOfSquares.CopositiveInnerCone` into a PSD-like inner cone plus a nonnegativity constraint on the strictly off-diagonal entries. The cone of inner-copositive matrices contains symmetric matrices that can @@ -14,7 +14,7 @@ entries. The bridge introduces constrained variables for each summand. `CopositiveInnerBridge` supports: - * [`SumOfSquares.CopositiveInnerCone{S}`](@ref) + * `SumOfSquares.CopositiveInnerCone{S}` ## Target nodes diff --git a/src/Bridges/Variable/kernel.jl b/src/Bridges/Variable/kernel.jl index 7c602e0ef..0f7a5b40a 100644 --- a/src/Bridges/Variable/kernel.jl +++ b/src/Bridges/Variable/kernel.jl @@ -19,7 +19,7 @@ preferred for solvers that natively support PSD as constrained variables for solvers that only support PSD as a constraint (e.g. Clarabel). The choice is driven by `MOI.Bridges.bridging_cost`. -`KernelBridge` does not handle [`MultivariateBases.LagrangeBasis`](@ref) — +`KernelBridge` does not handle `MultivariateBases.LagrangeBasis` — those bases dispatch to [`SumOfSquares.Bridges.Variable.LowRankBridge`](@ref) instead. From 40e9bec69ff9a3c774e0e75407b75d00aaaa9d8e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 19:39:51 +0200 Subject: [PATCH 13/17] Add tests --- .../sos_polynomial_in_semialgebraic_set.jl | 179 ++++++++++++++++++ 1 file changed, 179 insertions(+) create mode 100644 test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl diff --git a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl new file mode 100644 index 000000000..b58c35797 --- /dev/null +++ b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -0,0 +1,179 @@ +module TestConstraintSOSPolynomialInSemialgebraicSet + +using Test +import MathOptInterface as MOI +import MultivariateBases as MB +import MultivariatePolynomials as MP +import StarAlgebras as SA +using DynamicPolynomials +using SumOfSquares + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +# Build a `Putinar(Newton, Newton, maxdegree)` certificate, the default that +# `@constraint(model, ... in SOSCone(), domain = K, maxdegree = d)` produces +# under the hood. +function _newton_putinar(vars, maxdegree) + full_basis = MB.FullBasis{MB.Monomial}(vars) + newton = SumOfSquares.Certificate.Newton( + SOSCone(), + full_basis, + full_basis, + SumOfSquares.Certificate.NewtonFilter( + SumOfSquares.Certificate.NewtonDegreeBounds(tuple()), + ), + ) + return SumOfSquares.Certificate.Putinar(newton, newton, maxdegree) +end + +# Cover the core code path of the bridge: σ_0's gram basis is determined by +# `Newton` on the *augmented* polynomial `p + Σ g_i · b_i ⋆ b_i` (a recent +# bugfix; without the augmentation, Newton on `p` alone gives a different +# basis), and a single Lagrangian multiplier σ_1 is added for the +# inequality `x ≥ 0`. +# +# Input polynomial `p(x) = (1 - a) x + a x^2` with `basis = [x, x^2]`, +# `domain = {x ≥ 0}` and `maxdegree = 2`. +# - multiplier basis for `g_1 = x`: `[1]` +# - augmented poly monomials: `{x, x^2}` (Newton polytope ⇒ gram basis `[x]`) +# - σ_0 gram basis: `[x]` +# - σ_1 gram basis: `[1]` +# - WeightedSOSCone weights: `[g_1, 1] = [x, 1]` +# - `new_basis = [x, x^2]` (coefficients align with the input function). +function test_runtests() + T = Float64 + @polyvar x + basis = MB.SubBasis{MB.Monomial}([x, x^2]) + domain = (@set x >= 0) + cert = _newton_putinar([x], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + sigma_0_basis = MB.SubBasis{MB.Monomial}([x]) + sigma_1_basis = MB.SubBasis{MB.Monomial}([x^0]) + new_basis = MB.SubBasis{MB.Monomial}([x, x^2]) + implicit_basis = MB.implicit_basis(basis) + w_unit = MB.algebra_element( + MB.sparse_coefficients(MP.polynomial(MP.term(one(T), MP.constant_monomial(x)))), + implicit_basis, + ) + g1 = MP.polynomial(domain.p[1]) + w_g1 = MB.algebra_element( + MB.sparse_coefficients(one(T) * MP.similar(g1, T)), + implicit_basis, + ) + weighted = SumOfSquares.WeightedSOSCone{ + MOI.PositiveSemidefiniteConeTriangle, + }( + new_basis, + [sigma_1_basis, sigma_0_basis], + [w_g1, w_unit], + ) + MOI.Bridges.runtests( + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, + model -> begin + a = MOI.add_variable(model) + MOI.add_constraint( + model, + MOI.Utilities.operate( + vcat, + T, + T(1) - T(1) * a, + T(0) + T(1) * a, + ), + set, + ) + end, + model -> begin + a = MOI.add_variable(model) + MOI.add_constraint( + model, + MOI.Utilities.operate( + vcat, + T, + T(1) - T(1) * a, + T(0) + T(1) * a, + ), + weighted, + ) + end; + cannot_unbridge = true, + ) + return +end + +# Two inequalities to exercise the `multiplier_indices` machinery for σ_0 +# (last) and σ_1, σ_2 (first two `gram_bases`). The augmentation now adds +# monomials from *both* `g_i b_i ⋆ b_i` products to the polynomial used to +# compute σ_0's Newton polytope. +function test_runtests_two_inequalities() + T = Float64 + @polyvar x y + basis = MB.SubBasis{MB.Monomial}([x * y]) + domain = @set x >= 0 && y >= 0 + cert = _newton_putinar([x, y], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + # With the constant input function `1`, Newton's filter on the augmented + # polynomial gives empty multiplier gram bases — what matters for this + # test is that we exercise the multi-inequality `for index in preorder_indices` + # loop and the `multiplier_indices` bookkeeping. We have to use the + # same parent `FullBasis{Monomial}([x, y])` as the certificate's gram + # basis (the parent is part of the `SubBasis`'s identity). + sigma_basis = SA.SubBasis( + MB.FullBasis{MB.Monomial}([x, y]), + Vector{Vector{Int}}(), + ) + new_basis = MB.SubBasis{MB.Monomial}([x * y]) + implicit_basis = MB.implicit_basis(basis) + w_unit = MB.algebra_element( + MB.sparse_coefficients(MP.polynomial(MP.term(one(T), MP.constant_monomial(x * y)))), + implicit_basis, + ) + w_g1 = MB.algebra_element( + MB.sparse_coefficients(one(T) * MP.similar(MP.polynomial(domain.p[1]), T)), + implicit_basis, + ) + w_g2 = MB.algebra_element( + MB.sparse_coefficients(one(T) * MP.similar(MP.polynomial(domain.p[2]), T)), + implicit_basis, + ) + weighted = SumOfSquares.WeightedSOSCone{ + MOI.PositiveSemidefiniteConeTriangle, + }(new_basis, [sigma_basis, sigma_basis, sigma_basis], [w_g1, w_g2, w_unit]) + MOI.Bridges.runtests( + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, + model -> begin + MOI.add_constraint( + model, + MOI.VectorAffineFunction{T}( + MOI.VectorAffineTerm{T}[], + T[1], + ), + set, + ) + end, + model -> begin + MOI.add_constraint( + model, + MOI.VectorAffineFunction{T}( + MOI.VectorAffineTerm{T}[], + T[1], + ), + weighted, + ) + end; + cannot_unbridge = true, + ) + return +end + +end # module + +TestConstraintSOSPolynomialInSemialgebraicSet.runtests() From a7f6103422e2e925ec6936c28a30f83071da3d99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 19:39:59 +0200 Subject: [PATCH 14/17] Fix format --- .../sos_polynomial_in_semialgebraic_set.jl | 54 ++++++++++--------- 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index b58c35797..a15f96c7c 100644 --- a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -61,7 +61,9 @@ function test_runtests() new_basis = MB.SubBasis{MB.Monomial}([x, x^2]) implicit_basis = MB.implicit_basis(basis) w_unit = MB.algebra_element( - MB.sparse_coefficients(MP.polynomial(MP.term(one(T), MP.constant_monomial(x)))), + MB.sparse_coefficients( + MP.polynomial(MP.term(one(T), MP.constant_monomial(x))), + ), implicit_basis, ) g1 = MP.polynomial(domain.p[1]) @@ -69,13 +71,12 @@ function test_runtests() MB.sparse_coefficients(one(T) * MP.similar(g1, T)), implicit_basis, ) - weighted = SumOfSquares.WeightedSOSCone{ - MOI.PositiveSemidefiniteConeTriangle, - }( - new_basis, - [sigma_1_basis, sigma_0_basis], - [w_g1, w_unit], - ) + weighted = + SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( + new_basis, + [sigma_1_basis, sigma_0_basis], + [w_g1, w_unit], + ) MOI.Bridges.runtests( SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, model -> begin @@ -126,46 +127,47 @@ function test_runtests_two_inequalities() # loop and the `multiplier_indices` bookkeeping. We have to use the # same parent `FullBasis{Monomial}([x, y])` as the certificate's gram # basis (the parent is part of the `SubBasis`'s identity). - sigma_basis = SA.SubBasis( - MB.FullBasis{MB.Monomial}([x, y]), - Vector{Vector{Int}}(), - ) + sigma_basis = + SA.SubBasis(MB.FullBasis{MB.Monomial}([x, y]), Vector{Vector{Int}}()) new_basis = MB.SubBasis{MB.Monomial}([x * y]) implicit_basis = MB.implicit_basis(basis) w_unit = MB.algebra_element( - MB.sparse_coefficients(MP.polynomial(MP.term(one(T), MP.constant_monomial(x * y)))), + MB.sparse_coefficients( + MP.polynomial(MP.term(one(T), MP.constant_monomial(x * y))), + ), implicit_basis, ) w_g1 = MB.algebra_element( - MB.sparse_coefficients(one(T) * MP.similar(MP.polynomial(domain.p[1]), T)), + MB.sparse_coefficients( + one(T) * MP.similar(MP.polynomial(domain.p[1]), T), + ), implicit_basis, ) w_g2 = MB.algebra_element( - MB.sparse_coefficients(one(T) * MP.similar(MP.polynomial(domain.p[2]), T)), + MB.sparse_coefficients( + one(T) * MP.similar(MP.polynomial(domain.p[2]), T), + ), implicit_basis, ) - weighted = SumOfSquares.WeightedSOSCone{ - MOI.PositiveSemidefiniteConeTriangle, - }(new_basis, [sigma_basis, sigma_basis, sigma_basis], [w_g1, w_g2, w_unit]) + weighted = + SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( + new_basis, + [sigma_basis, sigma_basis, sigma_basis], + [w_g1, w_g2, w_unit], + ) MOI.Bridges.runtests( SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, model -> begin MOI.add_constraint( model, - MOI.VectorAffineFunction{T}( - MOI.VectorAffineTerm{T}[], - T[1], - ), + MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1]), set, ) end, model -> begin MOI.add_constraint( model, - MOI.VectorAffineFunction{T}( - MOI.VectorAffineTerm{T}[], - T[1], - ), + MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1]), weighted, ) end; From eb2e72e4eff5e76794cf5c8bd98a86cac626563c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 18 Jun 2026 22:21:19 +0200 Subject: [PATCH 15/17] add tests --- .../sos_polynomial_in_semialgebraic_set.jl | 125 ++++++++++++++++++ 1 file changed, 125 insertions(+) diff --git a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index a15f96c7c..95e4ad742 100644 --- a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -35,6 +35,31 @@ function _newton_putinar(vars, maxdegree) return SumOfSquares.Certificate.Putinar(newton, newton, maxdegree) end +# Build a `Sparsity.Preorder(Variable, Putinar(Newton, MaxDegree, maxdegree))` +# certificate, matching what `sparsity = Sparsity.Variable()` produces in +# JuMP. +function _sparse_putinar(vars, maxdegree) + full_basis = MB.FullBasis{MB.Monomial}(vars) + newton = SumOfSquares.Certificate.Newton( + SOSCone(), + full_basis, + full_basis, + SumOfSquares.Certificate.NewtonFilter( + SumOfSquares.Certificate.NewtonDegreeBounds(tuple()), + ), + ) + max_deg = SumOfSquares.Certificate.MaxDegree( + SOSCone(), + full_basis, + full_basis, + maxdegree, + ) + return SumOfSquares.Certificate.Sparsity.Preorder( + SumOfSquares.Certificate.Sparsity.Variable(), + SumOfSquares.Certificate.Putinar(newton, max_deg, maxdegree), + ) +end + # Cover the core code path of the bridge: σ_0's gram basis is determined by # `Newton` on the *augmented* polynomial `p + Σ g_i · b_i ⋆ b_i` (a recent # bugfix; without the augmentation, Newton on `p` alone gives a different @@ -176,6 +201,106 @@ function test_runtests_two_inequalities() return end +# Exercise the `if raw isa AbstractVector` branch of the `raw_ranges` +# bookkeeping: when the certificate carries sparsity (here +# `Sparsity.Preorder(Variable, Putinar(Newton, MaxDegree, …))`), both the σ_i +# multiplier bases and σ_0's basis come back as `Vector{SubBasis}` and the +# bridge collapses them into a `UnitRange` slot of `multiplier_indices`. For +# `p = x^2 + y^2` over `domain = {x ≥ 0}`, the σ_0 basis is split into two +# sparsity blocks `[1, y]` and `[1, x]` (one per disconnected variable +# cluster), so `multiplier_indices = [2:3, 1]`. +function test_runtests_sparsity() + T = Float64 + @polyvar x y + basis = MB.SubBasis{MB.Monomial}([y^2, x^2]) + domain = (@set x >= 0) + cert = _sparse_putinar([x, y], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + full = MB.FullBasis{MB.Monomial}([x, y]) + # Construct the basis and gram bases with the same parent + # `FullBasis{Monomial}([x, y])` as the certificate uses, so that + # `==` on `SubBasis` returns true. + new_basis = SA.SubBasis(full, [[0, 0], [0, 1], [1, 0], [0, 2], [2, 0]]) + sigma_1_basis = SA.SubBasis(full, [[0, 0]]) + sigma_0_block_y = SA.SubBasis(full, [[0, 0], [0, 1]]) + sigma_0_block_x = SA.SubBasis(full, [[0, 0], [1, 0]]) + implicit_basis = MB.implicit_basis(basis) + w_unit = MB.algebra_element( + MB.sparse_coefficients( + MP.polynomial(MP.term(one(T), MP.constant_monomial(x * y))), + ), + implicit_basis, + ) + # `domain.p[1] = x` only carries the variable `x`, so we lift it into + # the `[x, y]` variable system with `+ 0*y` to match the bridge's + # `Certificate.generator(...)` output exactly. + g1_full = MP.polynomial(domain.p[1] + 0 * y) + w_g1 = MB.algebra_element( + MB.sparse_coefficients(one(T) * MP.similar(g1_full, T)), + implicit_basis, + ) + weighted = + SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( + new_basis, + [sigma_1_basis, sigma_0_block_y, sigma_0_block_x], + [w_g1, w_unit, w_unit], + ) + func = MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1, 1]) + new_func = + MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[0, 0, 0, 1, 1]) + MOI.Bridges.runtests( + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge, + model -> begin + MOI.add_constraint(model, func, set) + end, + model -> begin + MOI.add_constraint(model, new_func, weighted) + end; + cannot_unbridge = true, + ) + return +end + +# Verify the `multiplier_indices` shape that the `GramMatrixAttribute`, +# `MomentMatrixAttribute` and `LagrangianMultipliers` getters dispatch on: +# - σ_0 at `multiplier_indices[1]` must be a `UnitRange` (the two sparsity +# blocks ⇒ `_get` ⇒ `BlockDiagonalGramMatrix`); +# - σ_1 at `multiplier_indices[2]` must also be a `UnitRange` (single +# block wrapped in a `Vector` ⇒ same `BlockDiagonalGramMatrix` path). +# The actual gram-matrix recovery on a real-solver path is covered by +# `test/Tests/rearrangement.jl` (run from `test/Mock/rearrangement.jl`). +function test_multiplier_indices_sparsity() + T = Float64 + @polyvar x y + basis = MB.SubBasis{MB.Monomial}([y^2, x^2]) + domain = (@set x >= 0) + cert = _sparse_putinar([x, y], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + # Build the bridged model with a mock optimizer that we can poke. + mock = MOI.Utilities.MockOptimizer( + MOI.Utilities.UniversalFallback(MOI.Utilities.Model{T}()), + ) + wrapped = MOI.Bridges.Constraint.SingleBridgeOptimizer{ + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge{ + T, + }, + }( + mock, + ) + func = MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1, 1]) + ci = MOI.add_constraint(wrapped, func, set) + bridge = MOI.Bridges.bridge(wrapped, ci) + # σ_0 occupies `multiplier_indices[1]` and is split into 2 sparsity + # blocks, hence the `UnitRange` `2:3`. σ_1 has a single sparsity block + # but the sparsity certificate still wraps its basis in a `Vector`, so + # it comes back as the (length-1) `UnitRange` `1:1` instead of a plain + # `Int` — that single-block-but-Vector path is precisely the case the + # `BlockDiagonalGramMatrix` getter has to handle. + @test bridge.multiplier_indices[1] == 2:3 + @test bridge.multiplier_indices[2] == 1:1 + return +end + end # module TestConstraintSOSPolynomialInSemialgebraicSet.runtests() From d72af98e64347ac5380dd505fd919ac5b7b5fef5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 19 Jun 2026 09:47:54 +0200 Subject: [PATCH 16/17] Fixes --- src/Bridges/Constraint/image.jl | 10 +- .../sos_polynomial_in_semialgebraic_set.jl | 4 + src/constraints.jl | 25 ++- test/Bridges/lazy.jl | 155 ++++++++++++------ 4 files changed, 133 insertions(+), 61 deletions(-) diff --git a/src/Bridges/Constraint/image.jl b/src/Bridges/Constraint/image.jl index 5026bff05..8def8cf98 100644 --- a/src/Bridges/Constraint/image.jl +++ b/src/Bridges/Constraint/image.jl @@ -354,9 +354,13 @@ end function MOI.supports_constraint( ::Type{ImageBridge{T}}, ::Type{<:MOI.AbstractVectorFunction}, - ::Type{<:SOS.WeightedSOSCone}, -) where {T} - return true + ::Type{<:SOS.WeightedSOSCone{M,B}}, +) where {T,M,B} + # `LagrangeBasis` dispatches to `Variable.LowRankBridge`; the image-form + # bridge cannot iterate over its samples (`bridge_constraint` walks the + # source basis monomial-by-monomial), so excluding it here both fixes the + # runtime error and lets the cost graph correctly favour `LowRankBridge`. + return !(B <: MB.LagrangeBasis) end function MOI.Bridges.added_constrained_variable_types(::Type{<:ImageBridge}) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 276cc1865..029239ca2 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -55,6 +55,10 @@ function MOI.Bridges.Constraint.bridge_constraint( set::SOS.SOSPolynomialSet{<:SemialgebraicSets.BasicSemialgebraicSet}, ) where {T,F,DT,CT,BT,M,NB,GB,W} @assert MOI.output_dimension(f) == length(set.basis) + # MOI does not modify the coefficients of the functions so we can modify `p`. + # without altering `f`. + # The monomials may be copied by MA however so we need to copy it. + # TODO remove `collect` when `DynamicPolynomials.MonomialVector` can be used as keys poly = MB.algebra_element( SA.SparseCoefficients( copy(collect(set.basis.keys)), diff --git a/src/constraints.jl b/src/constraints.jl index ad52be71d..129ce2059 100644 --- a/src/constraints.jl +++ b/src/constraints.jl @@ -334,12 +334,25 @@ end function PolyJuMP.bridges( F::Type{<:MOI.AbstractVectorFunction}, - ::Type{<:WeightedSOSCone}, -) # Needed so that `Variable.KernelBridge` is added as well - return Tuple{Type,Type}[( - MOI.Bridges.Constraint.VectorSlackBridge, - PolyJuMP.coefficient_type_or_float(F), - )] + ::Type{<:WeightedSOSCone{M,B}}, +) where {M,B} + T = PolyJuMP.coefficient_type_or_float(F) + # `VectorSlackBridge` is needed so that `Variable.KernelBridge` / + # `Variable.LowRankBridge` (registered by + # `PolyJuMP.bridges(::Type{<:WeightedSOSCone})`) are added as well via the + # variable-bridge chain. `Bridges.Constraint.ImageBridge` is the + # constraint-side counterpart; it must be listed here too so that JuMP + # recursively registers the bridges for ImageBridge's target types (e.g. + # `Constraint.PositiveSemidefinite2x2Bridge`) — otherwise + # `MOI.Bridges.bridging_cost` overestimates the constraint-side cost and + # picks the variable-side path even when `ImageBridge` should win. + # `ImageBridge` doesn't handle `LagrangeBasis`; mirror its + # `supports_constraint` check so JuMP only registers it where applicable. + bridges = Tuple{Type,Type}[(MOI.Bridges.Constraint.VectorSlackBridge, T)] + if !(B <: MB.LagrangeBasis) + push!(bridges, (Bridges.Constraint.ImageBridge, T)) + end + return bridges end function _bridge_coefficient_type( diff --git a/test/Bridges/lazy.jl b/test/Bridges/lazy.jl index 3bccb5a96..f4c5446c5 100644 --- a/test/Bridges/lazy.jl +++ b/test/Bridges/lazy.jl @@ -4,9 +4,27 @@ using Test import Clarabel using DynamicPolynomials import Dualization +using JuMP import MultivariateBases as MB +import MultivariatePolynomials as MP using SumOfSquares +# Build a `Putinar(Newton, Newton, maxdegree)` certificate, matching what +# `@constraint(model, ... in SOSCone(), domain = K, maxdegree = d)` produces +# under the hood. +function _newton_putinar(vars, maxdegree) + full_basis = MB.FullBasis{MB.Monomial}(vars) + newton = SumOfSquares.Certificate.Newton( + SOSCone(), + full_basis, + full_basis, + SumOfSquares.Certificate.NewtonFilter( + SumOfSquares.Certificate.NewtonDegreeBounds(tuple()), + ), + ) + return SumOfSquares.Certificate.Putinar(newton, newton, maxdegree) +end + function runtests() for name in names(@__MODULE__; all = true) if startswith("$(name)", "test_") @@ -67,11 +85,9 @@ end # Wrapping `Clarabel.Optimizer` in `Dualization.dual_optimizer` flips # Clarabel's `VAF`-in-`PositiveSemidefiniteConeTriangle` support into # `PositiveSemidefiniteConeTriangle` supported as constrained variables on -# the outer bridge graph. `Constraint.ImageBridge` would still win on cost -# (the cost bump in `image.jl` only switches when the *immediate* -# `PositiveSemidefiniteConeTriangle` is variable-side, which the dualization -# wrapper hides one level deeper), so we explicitly remove it to expose the -# variable-side route: `VectorSlackBridge` lifts the `VAF`-in-`WeightedSOSCone` +# the outer bridge graph. The cost graph propagates the inner solver's +# support through the dualization layer, so the variable-side route wins +# naturally: `VectorSlackBridge` lifts the `VAF`-in-`WeightedSOSCone` # constraint into constrained variables in `WeightedSOSCone`, which # `Variable.KernelBridge` then breaks into PSD blocks. function test_dual_clarabel_uses_kernel_bridge() @@ -83,10 +99,6 @@ function test_dual_clarabel_uses_kernel_bridge() with_bridge_type = T, ) SumOfSquares.Bridges.add_all_bridges(optimizer, T) - MOI.Bridges.remove_bridge( - optimizer, - SumOfSquares.Bridges.Constraint.ImageBridge{T}, - ) set = SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( MB.SubBasis{MB.Monomial}([y^4, x * y^3, x^2 * y^2, x^3 * y, x^4]), [MB.SubBasis{MB.Monomial}([y^2, x * y, x^2])], @@ -104,33 +116,87 @@ function test_dual_clarabel_uses_kernel_bridge() return end -# With a `domain` keyword, the constraint enters the bridge graph as a -# `SOSPolynomialSet{<:BasicSemialgebraicSet}`, and is first reformulated by -# `SOSPolynomialInSemialgebraicSetBridge` into a single `WeightedSOSCone`. -# From there, Clarabel — which natively supports +# With a `domain` keyword, the constraint enters the bridge graph as +# `SOSPolynomialSet{<:BasicSemialgebraicSet}`, which +# `SOSPolynomialInSemialgebraicSetBridge` reformulates into a single +# `WeightedSOSCone`. From there, Clarabel — which natively supports # `MOI.PositiveSemidefiniteConeTriangle` as a constraint — routes through -# `VectorSlackBridge` and `Variable.KernelBridge` rather than through -# `Constraint.ImageBridge`, because the bump in `ImageBridge`'s cost shifts -# the balance once the immediate PSD lands on the variable side. -function test_clarabel_with_domain_uses_kernel_bridge() +# `Constraint.ImageBridge`, mirroring `test_clarabel_uses_image_bridge`. +function test_clarabel_with_domain_uses_image_bridge() T = Float64 - @polyvar x y + @polyvar x + optimizer = MOI.instantiate( + Clarabel.Optimizer; + with_cache_type = T, + with_bridge_type = T, + ) + SumOfSquares.Bridges.add_all_bridges(optimizer, T) + basis = MB.SubBasis{MB.Monomial}([x, x^2]) + domain = (@set x >= 0) + cert = _newton_putinar([x], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + F = MOI.VectorAffineFunction{T} + S = typeof(set) + # `SOSPolynomialInSemialgebraicSetBridge` (cost 1) + `ImageBridge` path + # (cost 6 for Clarabel; same as `test_clarabel_uses_image_bridge`). + @test MOI.Bridges.bridging_cost(optimizer, F, S) == 7.0 + func = MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1, 1]) + ci = MOI.add_constraint(optimizer, func, set) + outer = MOI.Bridges.bridge(optimizer, ci) + @test outer isa + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge + inner = MOI.Bridges.bridge(optimizer, outer.constraint) + @test inner isa SumOfSquares.Bridges.Constraint.ImageBridge + return +end + +# Wrapping `Clarabel.Optimizer` in `Dualization.dual_optimizer` flips +# Clarabel's `VAF`-in-`PositiveSemidefiniteConeTriangle` support into +# `PositiveSemidefiniteConeTriangle` supported as constrained variables on +# the outer bridge graph. As in `test_dual_clarabel_uses_kernel_bridge`, +# the cost graph propagates this through the layers and the variable-side +# route wins naturally: `SOSPolynomialInSemialgebraicSetBridge` → +# `VectorSlackBridge` → `Variable.KernelBridge`. +function test_dual_clarabel_with_domain_uses_kernel_bridge() + T = Float64 + @polyvar x + optimizer = MOI.instantiate( + Dualization.dual_optimizer(Clarabel.Optimizer); + with_cache_type = T, + with_bridge_type = T, + ) + SumOfSquares.Bridges.add_all_bridges(optimizer, T) + basis = MB.SubBasis{MB.Monomial}([x, x^2]) + domain = (@set x >= 0) + cert = _newton_putinar([x], 2) + set = SumOfSquares.SOSPolynomialSet(domain, basis, cert) + func = MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1, 1]) + ci = MOI.add_constraint(optimizer, func, set) + outer = MOI.Bridges.bridge(optimizer, ci) + @test outer isa + SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge + slack = MOI.Bridges.bridge(optimizer, outer.constraint) + @test slack isa MOI.Bridges.Constraint.VectorSlackBridge + inner = MOI.Bridges.bridge(optimizer, slack.slack_in_set) + @test inner isa SumOfSquares.Bridges.Variable.KernelBridge + return +end + +# Walk the bridge chain on a JuMP model directly (no `MOI.instantiate`, +# no explicit `add_all_bridges` call). `PolyJuMP.bridges` is responsible +# for transitively registering every bridge the chain might use; if any +# is missing, the lazy cost graph picks the wrong path. The JuMP-side +# chain should match `test_clarabel_with_domain_uses_image_bridge`. +function test_clarabel_with_domain_uses_image_bridge_jump() + @polyvar x model = Model(Clarabel.Optimizer) set_silent(model) @variable(model, α) @objective(model, Max, α) S = @set x >= 0 - cref = @constraint( - model, - 10 - x^2 - α * y in SOSCone(), - domain = S, - maxdegree = 4, - ) + @constraint(model, α - x^2 in SOSCone(), domain = S, maxdegree = 2) optimize!(model) - backend = JuMP.backend(model) - # Walk down the chain: SOSPolynomialSet → InSemialgebraicSetBridge → - # WeightedSOSCone (via VectorSlackBridge) → KernelBridge. - lbo = backend.optimizer + lbo = JuMP.backend(model).optimizer cis = MOI.ConstraintIndex[] for (F, S) in MOI.get(lbo, MOI.ListOfConstraintTypesPresent()) if S <: SumOfSquares.SOSPolynomialSet @@ -141,40 +207,25 @@ function test_clarabel_with_domain_uses_kernel_bridge() outer = MOI.Bridges.bridge(lbo, only(cis)) @test outer isa SumOfSquares.Bridges.Constraint.SOSPolynomialInSemialgebraicSetBridge - slack = MOI.Bridges.bridge(backend.optimizer, outer.constraint) - @test slack isa MOI.Bridges.Constraint.VectorSlackBridge - inner = MOI.Bridges.bridge(backend.optimizer, slack.slack_in_set) - @test inner isa SumOfSquares.Bridges.Variable.KernelBridge + inner = MOI.Bridges.bridge(lbo, outer.constraint) + @test inner isa SumOfSquares.Bridges.Constraint.ImageBridge return end -# With `Dualization.dual_optimizer`, the role of -# `PositiveSemidefiniteConeTriangle` flips between primal/dual. The chain -# is otherwise the same as `test_clarabel_with_domain_uses_kernel_bridge`: +# JuMP counterpart of `test_dual_clarabel_with_domain_uses_kernel_bridge`: +# wrap Clarabel in `Dualization.dual_optimizer` and verify the chain is # `SOSPolynomialInSemialgebraicSetBridge` → `VectorSlackBridge` → -# `Variable.KernelBridge`. `Constraint.ImageBridge` is not added to the -# graph at all here, because the dualization wrapper hides the -# `PositiveSemidefiniteConeTriangle` one level deeper than the cost -# heuristic in `image.jl` looks (a follow-up could tune the bridging cost -# so that `ImageBridge` wins on the dual side, like it does for -# `Variable.KernelBridge` on the primal side). -function test_dual_clarabel_with_domain_uses_kernel_bridge() - T = Float64 - @polyvar x y +# `Variable.KernelBridge` without any explicit bridge removal. +function test_dual_clarabel_with_domain_uses_kernel_bridge_jump() + @polyvar x model = Model(Dualization.dual_optimizer(Clarabel.Optimizer)) set_silent(model) @variable(model, α) @objective(model, Max, α) S = @set x >= 0 - cref = @constraint( - model, - 10 - x^2 - α * y in SOSCone(), - domain = S, - maxdegree = 4, - ) + @constraint(model, α - x^2 in SOSCone(), domain = S, maxdegree = 2) optimize!(model) - backend = JuMP.backend(model) - lbo = backend.optimizer + lbo = JuMP.backend(model).optimizer cis = MOI.ConstraintIndex[] for (F, S) in MOI.get(lbo, MOI.ListOfConstraintTypesPresent()) if S <: SumOfSquares.SOSPolynomialSet From 9d4ceb7e8aa617d2d8a189d797bad6ba9c370a33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 19 Jun 2026 13:50:20 +0200 Subject: [PATCH 17/17] Keep sigma0 from newton polytope --- .../sos_polynomial_in_semialgebraic_set.jl | 54 +++++-------------- src/Certificate/Sparsity/preorder.jl | 14 +++++ src/Certificate/preorder.jl | 48 +++++++++++++---- .../sos_polynomial_in_semialgebraic_set.jl | 4 +- 4 files changed, 68 insertions(+), 52 deletions(-) diff --git a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 029239ca2..08cc07964 100644 --- a/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/src/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -75,50 +75,18 @@ function MOI.Bridges.Constraint.bridge_constraint( implicit_basis = MB.implicit_basis(set.basis) preprocessed = Certificate.preprocessed_domain(set.certificate, set.domain, poly) + # σ_0's gram basis is obtained directly from the preorder certificate. + # `Putinar{<:Newton}` reuses the basis already computed by + # `half_newton_polytope` inside `preprocessed_domain` (which takes the + # inequality polynomials `g_i` into account); other certificates fall back + # to the ideal certificate's `gram_basis` on `poly_reduced`. + sigma_0_basis = + Certificate.gram_basis(set.certificate, preprocessed, poly_reduced) # Use the same `sparse_coefficients`-based construction for the unit # weight as for the polynomial g_i weights so they all share the same # concrete element type and fit in a homogeneous `Vector{W}`. some_mono = first(MB.keys_as_monomials(set.basis)) unit_poly = MP.polynomial(MP.term(one(T), MP.constant_monomial(some_mono))) - # Collect the multiplier bases (σ_1, σ_2, …) and the inequality - # polynomials g_i. - multiplier_bases = Any[] - g_polynomials = [] - for index in Certificate.preorder_indices(set.certificate, preprocessed) - push!( - multiplier_bases, - Certificate.multiplier_basis(set.certificate, index, preprocessed), - ) - push!( - g_polynomials, - Certificate.generator(set.certificate, index, preprocessed), - ) - end - # σ_0's gram basis must absorb both `p` and the monomials produced by - # `g_i * b_i_j * b_i_k`. We mirror the original two-bridge pipeline by - # symbolically augmenting `poly_reduced` with the unit-coefficient - # versions of those products before asking the ideal certificate for - # the gram basis — otherwise the Newton polytope would be that of `p` - # alone, which is generally too small to host an SOS decomposition. - augmented = SOS.MA.copy(poly_reduced) - for (mb, g_i) in zip(multiplier_bases, g_polynomials) - g_alg = MB.algebra_element( - MB.sparse_coefficients(one(T) * similar(g_i, T)), - implicit_basis, - ) - block_iter = mb isa AbstractVector ? mb : (mb,) - for block in block_iter - for b1 in block, b2 in block - term = MB.algebra_element(SA.star(b1) * b2) - MA.operate!(SA.UnsafeAddMul(*), augmented, term, g_alg) - end - end - end - MA.operate!(SA.canonical, SA.coeffs(augmented)) - sigma_0_basis = Certificate.gram_basis( - ideal_cert, - Certificate.with_variables(augmented, set.domain), - ) # Build `[g_1, …, 1]` weights and `[σ_1_basis, …, σ_0_basis]` gram # bases. σ_0 goes last so that downstream `Variable.KernelBridge` # allocates the σ_i gram variables before the σ_0 ones — matching the @@ -126,8 +94,12 @@ function MOI.Bridges.Constraint.bridge_constraint( # the order assumed by the cached `Mock/` tests). gram_bases_raw = Any[] weights_raw = W[] - for (mb, g_i) in zip(multiplier_bases, g_polynomials) - push!(gram_bases_raw, mb) + for index in Certificate.preorder_indices(set.certificate, preprocessed) + push!( + gram_bases_raw, + Certificate.multiplier_basis(set.certificate, index, preprocessed), + ) + g_i = Certificate.generator(set.certificate, index, preprocessed) push!( weights_raw, MB.algebra_element( diff --git a/src/Certificate/Sparsity/preorder.jl b/src/Certificate/Sparsity/preorder.jl index add48e88d..7fc4bddc6 100644 --- a/src/Certificate/Sparsity/preorder.jl +++ b/src/Certificate/Sparsity/preorder.jl @@ -85,6 +85,20 @@ function SumOfSquares.Certificate.generator( ) end +function SumOfSquares.Certificate.gram_basis( + certificate::Preorder, + preprocessed::Domain, + poly, +) + # σ_0's basis under sparsity is governed by the `Sparsity.Ideal` wrapper of + # the inner ideal certificate, which produces a `Vector` of bases (one per + # sparsity block). We forward `poly` so it can drive the sparsity pattern. + return SumOfSquares.Certificate.gram_basis( + SumOfSquares.Certificate.ideal_certificate(certificate), + SumOfSquares.Certificate.with_variables(poly, preprocessed.domain), + ) +end + function SumOfSquares.Certificate.ideal_certificate(certificate::Preorder) return Ideal( certificate.sparsity, diff --git a/src/Certificate/preorder.jl b/src/Certificate/preorder.jl index 5064b1b5b..fbb72d95b 100644 --- a/src/Certificate/preorder.jl +++ b/src/Certificate/preorder.jl @@ -50,6 +50,10 @@ _algebra_element(a::SA.AlgebraElement) = a struct WithFixedBases{S,B} inner::S bases::Vector{B} + # `σ_0`'s gram basis as computed by `half_newton_polytope` from `p` and + # the inequality polynomials `gs` (so it already accounts for the + # multiplier-times-g terms that show up in the SOS decomposition). + ideal_basis::B end # TODO temporary workaround because SS doesn't support AlgebraElement yet @@ -88,16 +92,14 @@ function with_fixed_basis( newton::AbstractNewtonPolytopeApproximation, ) v = with_variables(domain, p) - return WithFixedBases( - v.inner, - half_newton_polytope( - _algebra_element(p), - SemialgebraicSets.inequalities(v.inner), - v.variables, - maxdegree, - newton, - )[2], + ideal_basis, multiplier_bases = half_newton_polytope( + _algebra_element(p), + SemialgebraicSets.inequalities(v.inner), + v.variables, + maxdegree, + newton, ) + return WithFixedBases(v.inner, multiplier_bases, ideal_basis) end function preprocessed_domain( @@ -164,6 +166,34 @@ end ideal_certificate(certificate::Putinar) = certificate.ideal_certificate ideal_certificate(::Type{<:Putinar{MC,IC}}) where {MC,IC} = IC +""" + gram_basis( + certificate::AbstractPreorderCertificate, + preprocessed, + poly, + ) + +Return the gram basis of `σ_0`, the ideal multiplier in the Putinar-style +decomposition `p = σ_0 + Σ g_i σ_i`. + +Whenever the certificate already pre-computes `σ_0`'s basis as a side-effect +of [`preprocessed_domain`](@ref) (e.g. `Putinar{<:Newton}`, which gets it +from [`half_newton_polytope`](@ref)), the preprocessed value is returned +directly and `poly` is ignored. Otherwise (`Putinar{<:MaxDegree}` and +similar), the fallback computes the basis by calling +[`gram_basis`](@ref) on the inner [`ideal_certificate`](@ref) with `poly`. +""" +function gram_basis(::Putinar{<:Newton}, preprocessed::WithFixedBases, _poly) + return preprocessed.ideal_basis +end + +function gram_basis(certificate::Putinar, preprocessed::WithVariables, poly) + return gram_basis( + ideal_certificate(certificate), + with_variables(poly, preprocessed), + ) +end + function SumOfSquares.matrix_cone_type(::Type{<:Putinar{MC}}) where {MC} return SumOfSquares.matrix_cone_type(MC) end diff --git a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl index 95e4ad742..6933310da 100644 --- a/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl +++ b/test/Bridges/Constraint/sos_polynomial_in_semialgebraic_set.jl @@ -207,7 +207,7 @@ end # multiplier bases and σ_0's basis come back as `Vector{SubBasis}` and the # bridge collapses them into a `UnitRange` slot of `multiplier_indices`. For # `p = x^2 + y^2` over `domain = {x ≥ 0}`, the σ_0 basis is split into two -# sparsity blocks `[1, y]` and `[1, x]` (one per disconnected variable +# sparsity blocks `[1, x]` and `[1, y]` (one per disconnected variable # cluster), so `multiplier_indices = [2:3, 1]`. function test_runtests_sparsity() T = Float64 @@ -242,7 +242,7 @@ function test_runtests_sparsity() weighted = SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( new_basis, - [sigma_1_basis, sigma_0_block_y, sigma_0_block_x], + [sigma_1_basis, sigma_0_block_x, sigma_0_block_y], [w_g1, w_unit, w_unit], ) func = MOI.VectorAffineFunction{T}(MOI.VectorAffineTerm{T}[], T[1, 1])