Volume Objective
In this example, we compute the maximal (resp. minimal) volume ellipsoids and polynomial sublevel sets contained in (resp. containing) the square with vertices $(\pm 1, \pm 1)$. We start by defining the square with Polyhedra.
using Polyhedra
h = HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1)
p = polyhedron(h)Polyhedron DefaultPolyhedron{Float64, Polyhedra.Intersection{Float64, Vector{Float64}, Int64}, Polyhedra.Hull{Float64, Vector{Float64}, Int64}}:
4-element iterator of HalfSpace{Float64, Vector{Float64}}:
HalfSpace([1.0, 0.0], 1.0)
HalfSpace([-1.0, 0.0], 1.0)
HalfSpace([0.0, 1.0], 1.0)
HalfSpace([0.0, -1.0], 1.0)We need to pick an SDP solver, see here for a list of available ones.
using SetProg
using SumOfSquares
import Hypatia
sdp_solver = optimizer_with_attributes(Hypatia.Optimizer, MOI.Silent() => true)MathOptInterface.OptimizerWithAttributes(Hypatia.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute, Any}[MathOptInterface.Silent() => true])John ellipsoid
The maximal volume ellipsoid contained in a convex body is called its John ellipsoid. The John ellipsoid for our square can be computed as follows.
model = Model(sdp_solver)
@variable(model, john, Ellipsoid(symmetric=true, dimension=2))
@constraint(model, john ⊆ p)
@objective(model, Max, nth_root(volume(john)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
SetProg.Sets.print_support_function(value(john))solve_time(model) = 0.0016121864318847656
termination_status(model) = ALMOST_OPTIMAL
objective_value(model) = 0.9999999663952056
h(S, x) = x[2]^2 + x[1]^2Which cone does the solver receive ?
The objective nth_root(volume(john)) is reformulated by SetProg into a new variable t constrained by [t; vec(Q)] in MOI.RootDetConeTriangle(2) where Q is the matrix of the ellipsoid. We can see these reformulated constraints with list_of_constraint_types:
list_of_constraint_types(model)3-element Vector{Tuple{Type, Type}}:
(AffExpr, MathOptInterface.LessThan{Float64})
(Vector{VariableRef}, MathOptInterface.RootDetConeTriangle)
(Vector{VariableRef}, SumOfSquares.PositiveSemidefinite2x2ConeTriangle)The solver may not support every one of these constraints natively. In that case, JuMP transforms them using bridges. For instance, a solver only supporting positive semidefinite constraints would receive the RootDetConeTriangle reformulated with a GeometricMeanCone and a PositiveSemidefiniteConeTriangle, so with additional variables and constraints. Hypatia supports the root-determinant cone natively, and we can check that with print_active_bridges. It shows, for each constraint type of the model, which bridges are used and which constraints the solver receives in the end.
print_active_bridges(model) * Supported objective: MOI.VariableIndex
* Unsupported constraint: MOI.ScalarAffineFunction{Float64}-in-MOI.LessThan{Float64}
| bridged by:
| MOIB.Constraint.LessToGreaterBridge{Float64, MOI.ScalarAffineFunction{Float64}, MOI.ScalarAffineFunction{Float64}}
| may introduce:
| * Unsupported constraint: MOI.ScalarAffineFunction{Float64}-in-MOI.GreaterThan{Float64}
| | bridged by:
| | MOIB.Constraint.VectorizeBridge{Float64, MOI.VectorAffineFunction{Float64}, MOI.Nonnegatives, MOI.ScalarAffineFunction{Float64}}
| | may introduce:
| | * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Nonnegatives
* Unsupported constraint: MOI.VectorOfVariables-in-MOI.RootDetConeTriangle
| bridged by:
| MOIB.Constraint.SetDotScalingBridge{Float64, MOI.RootDetConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}
| may introduce:
| * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.RootDetConeTriangle}
* Unsupported constraint: MOI.VectorOfVariables-in-SumOfSquares.PositiveSemidefinite2x2ConeTriangle
| bridged by:
| SumOfSquares.Bridges.Constraint.PositiveSemidefinite2x2Bridge{Float64, MOI.VectorAffineFunction{Float64}}
| may introduce:
| * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.RotatedSecondOrderConeTo only look at the root-determinant cone, we give its function and set type:
print_active_bridges(model, Vector{VariableRef}, MOI.RootDetConeTriangle) * Unsupported constraint: MOI.VectorOfVariables-in-MOI.RootDetConeTriangle
| bridged by:
| MOIB.Constraint.SetDotScalingBridge{Float64, MOI.RootDetConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorOfVariables}
| may introduce:
| * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.RootDetConeTriangle}The only bridge is SetDotScalingBridge. It only rescales the off-diagonal entries into the MOI.Scaled{MOI.RootDetConeTriangle} set, which is the form in which Hypatia supports the cone natively. There is no GeometricMeanCone or PositiveSemidefiniteConeTriangle.
Log-determinant objective
Instead of the nth root of the determinant, we can also maximize its logarithm with log(volume(john)). The optimal ellipsoid is the same, but the objective value is now log(det(Q)) = log(1) = 0 instead of det(Q)^(1/2) = 1.
model = Model(sdp_solver)
@variable(model, john_log, Ellipsoid(symmetric=true, dimension=2))
@constraint(model, john_log ⊆ p)
@objective(model, Max, log(volume(john_log)))
optimize!(model)
@show termination_status(model)
@show objective_value(model)
SetProg.Sets.print_support_function(value(john_log))termination_status(model) = ALMOST_OPTIMAL
objective_value(model) = -1.9304391854351855e-8
h(S, x) = x[2]^2 + x[1]^2The objective is now reformulated into [t; 1; vec(Q)] in MOI.LogDetConeTriangle(2) which means t ≤ log(det(Q)). As the constant 1 is part of the function, it is a Vector{AffExpr} instead of a Vector{VariableRef}. Hypatia also supports this cone natively:
print_active_bridges(model, Vector{AffExpr}, MOI.LogDetConeTriangle) * Unsupported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.LogDetConeTriangle
| bridged by:
| MOIB.Constraint.SetDotScalingBridge{Float64, MOI.LogDetConeTriangle, MOI.VectorAffineFunction{Float64}, MOI.VectorAffineFunction{Float64}}
| may introduce:
| * Supported constraint: MOI.VectorAffineFunction{Float64}-in-MOI.Scaled{MOI.LogDetConeTriangle}Löwner ellipsoid
The minimal volume ellipsoid containing a convex body is called its Löwner ellipsoid. The Löwner ellipsoid for our square can be computed as follows.
model = Model(sdp_solver)
@variable(model, löwner, Ellipsoid(symmetric=true, dimension=2))
@constraint(model, p ⊆ löwner)
@objective(model, Min, nth_root(volume(löwner)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
löwner_value = value(löwner)SetProg.Sets.Ellipsoid{Float64}([0.49999998368834947 6.323710909172985e-18; 6.323710909172985e-18 0.4999999836865327])We can visualize the Löwner and John ellipsoids as follows.
using Plots
plot(ratio=:equal)
plot!(löwner_value)
plot!(p)
plot!(value(john))
Higher degree polynomials
Ellipsoids are the sublevel sets of positive definite quadratic forms. To allow for more sophisticated shapes, we instead look for sublevel sets of quartic forms. For this, we simply replace Ellipsoid(dimension=2) by PolySet(degree=4, dimension=2). Note that the quantities optimized are not exactly the volume anymore but provide a reasonable heuristic.
Maximal volume quartic sublevel set contained in the square
model = Model(sdp_solver)
@variable(model, quartic_inner, PolySet(degree=4, symmetric=true, convex=true))
@constraint(model, quartic_inner ⊆ p)
@objective(model, Max, nth_root(volume(quartic_inner)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
quartic_inner_value = value(quartic_inner)SetProg.Sets.Polar{Float64, SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}}(SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}(4, GramMatrix with row/column basis:
SubBasis{Monomial}([x[2]^2, x[1]*x[2], x[1]^2])
And entries in a 3×3 MultivariateMoments.SymMatrix{Float64}:
0.9999999853674684 -1.618967945420618e-11 -0.14542692691517345
-1.618967945420618e-11 3.119280865022582 -1.0583473835805435e-11
-0.14542692691517345 -1.0583473835805435e-11 0.999999985367468, [11.999999824409617 -9.713807672523707e-11 -9.713722415988084e-11 8.485281139359092; -9.713807672523707e-11 5.656854022384468 2.828426905409852 -6.34913988442545e-11; -9.713722415988084e-11 2.828426905409852 5.656854022384474 -6.349705117710575e-11; 8.485281139359092 -6.34913988442545e-11 -6.349705117710575e-11 11.999999824409617]))Minimal volume quartic sublevel set containing the square
model = Model(sdp_solver)
@variable(model, quartic_outer, PolySet(symmetric=true, degree=4, convex=true))
@constraint(model, p ⊆ quartic_outer)
@objective(model, Min, nth_root(volume(quartic_outer)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
quartic_outer_value = value(quartic_outer)SetProg.Sets.ConvexPolySet{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, SumOfSquares.GramMatrix{Float64, StarAlgebras.SubBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Int64, Vector{Int64}, StarAlgebras.MappedBasis{MultivariateBases.Polynomial{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}, Vector{Int64}}, Vector{Int64}, MultivariatePolynomials.ExponentsIterator{MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}, Nothing, Vector{Int64}}, MultivariateBases.Variables{MultivariateBases.Monomial, Vector{DynamicPolynomials.Variable{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, MultivariatePolynomials.Graded{MultivariatePolynomials.LexOrder}}}}, typeof(MultivariatePolynomials.exponents)}, Vector{Vector{Int64}}}, Float64, MultivariateMoments.SymMatrix{Float64}}, MultivariateMoments.SymMatrix{Float64}}(4, GramMatrix with row/column basis:
SubBasis{Monomial}([x[2]^2, x[1]*x[2], x[1]^2])
And entries in a 3×3 MultivariateMoments.SymMatrix{Float64}:
0.3017766934996899 -7.395550080495531e-13 -0.09365064570406897
-7.395550080495531e-13 0.5837479004421742 7.395714491831187e-13
-0.09365064570406897 7.395714491831187e-13 0.3017766935003648, [3.6213203219962793 -4.437330048297319e-12 -4.437081288240348e-12 1.4999999969172813; -4.437330048297319e-12 0.7928932180680723 0.08578643921886356 4.438415866830229e-12; -4.437081288240348e-12 0.08578643921886356 0.7928932180680737 4.437395306732085e-12; 1.4999999969172813 4.438415866830229e-12 4.437395306732085e-12 3.621320322004377])We can visualize the quartic sublevel sets as follows.
plot(ratio=:equal)
plot!(quartic_outer_value)
plot!(p)
plot!(quartic_inner_value)
plot!(value(john))
Inner sublevel sets of increasing degree
We can also explore how the inner sublevel set tightens as the degree grows, using the L1_heuristic as a tractable volume proxy.
function inner_L1(d)
model = Model(sdp_solver)
@variable(model, S, PolySet(symmetric=true, degree=d, convex=true))
@constraint(model, S ⊆ p)
@objective(model, Max, L1_heuristic(volume(S), [1.0, 1.0]))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
return value(S)
end
S2 = inner_L1(2)
S4 = inner_L1(4)
S6 = inner_L1(6)
S8 = inner_L1(8)
S10 = inner_L1(10)
plot(ratio=:equal)
plot!(quartic_outer_value)
plot!(p)
plot!(S10)
plot!(S8)
plot!(S6)
plot!(S4)
plot!(S2)
Non-homogeneous case
For non-symmetric bodies, the John/Löwner ellipsoids are not centered at the origin. We need to provide a point as a hint of an interior point. We illustrate on a horizontally-shifted square.
shift = 1.1
h_shift = HalfSpace([1, 0], 1.0 + shift) ∩ HalfSpace([-1, 0], 1.0 - shift) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1)
p_shift = polyhedron(h_shift)
model = Model(sdp_solver)
@variable(model, john_shift, Ellipsoid(point=SetProg.InteriorPoint([shift, 0.0])))
@constraint(model, john_shift ⊆ p_shift)
@objective(model, Max, nth_root(volume(john_shift)))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
plot(ratio=:equal)
plot!(p_shift)
plot!(value(john_shift))
Likewise for the quartic sublevel set.
model = Model(sdp_solver)
@variable(model, quartic_shift, PolySet(convex=true, degree=2, point=SetProg.InteriorPoint([shift, 0.0])))
@constraint(model, quartic_shift ⊆ p_shift)
@objective(model, Max, L1_heuristic(volume(quartic_shift), [1.0, 1.0]))
optimize!(model)
@show solve_time(model)
@show termination_status(model)
@show objective_value(model)
plot(ratio=:equal)
plot!(p_shift)
plot!(value(quartic_shift))
This page was generated using Literate.jl.