Largely inspired from tools/examples/fusion/python/lownerjohn_ellipsoid.py
in Mosek
Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope.
function _ilog2(n, i)
if n <= (one(n) << i)
i
else
_ilog2(n, i+1)
end
end
function ilog2(n::Integer)
@assert n > zero(n)
_ilog2(n, zero(n))
end
@assert ilog2(1) == 0
@assert ilog2(2) == 1
@assert ilog2(3) == 2
@assert ilog2(4) == 2
@assert ilog2(5) == 3
@assert ilog2(6) == 3
@assert ilog2(7) == 3
@assert ilog2(8) == 3
using JuMP # Needs branch jump/moi
using MathOptInterface
const MOI = MathOptInterface
"""
geometric_mean(model, x, t)
Models the convex set
S = { (x, t) \in R^n x R | x >= 0, t <= (x1 * x2 * ... * xn)^(1/n) }
as the intersection of rotated quadratic cones and affine hyperplanes.
see [1, p. 105] or [2, p. 21]. This set can be interpreted as the hypograph of the
geometric mean of x.
We illustrate the modeling procedure using the following example.
Suppose we have
t <= (x1 * x2 * x3)^(1/3)
for some t >= 0, x >= 0. We rewrite it as
t^4 <= x1 * x2 * x3 * x4, x4 = t
which is equivalent to (see [1])
x11^2 <= 2*x1*x2, x12^2 <= 2*x3*x4,
x21^2 <= 2*x11*x12,
x21 = sqrt(8)*t, x4 = t.
"""
function geometric_mean(model, x, t)
function rec(x)
n = length(x)
if n > 1
y = @variable model [1:div(n, 2)]
for i in 1:length(y)
@constraint model [x[2*i-1], x[2*i], y[i]] in MOI.RotatedSecondOrderCone(3)
end
rec(y)
else
x[1]
end
end
n = length(x)
l = ilog2(n)
m = (1 << l) - n
# if size of x is not a power of 2 we pad it:
if m > 0
x_padding = @variable model [1:m]
# set the last m elements equal to t
@constraint model x_padding .== t
# The type of x could be Vector{AffExpr} in which case
# x = [x; x_padding] would change the type of x
y = [x; x_padding]
else
y = x
end
z = rec(y)
@constraint model sqrt(2)^l * t == z
end
"""
det_rootn(model, X, t)
Purpose: Models the hypograph of the n-th power of the
determinant of a positive definite matrix. See [1, p. 149] or [2, p. 42] for more details.
The convex set (a hypograph)
```math
C = \\{ (X, t) \\in S^n_+ \\times R \\mid t \\leq \\det(X)^{1/n} \\},
```
can be modeled as the intersection of a semidefinite cone
```math
\\begin{pmatrix}
X & Δ\\\\
Δ^\\top & \\mathrm{Diag}(Δ)
\\end{\pmatrix} \\geq 0
```
and a number of rotated quadratic cones and affine hyperplanes,
```math
t \\leq (Δ_{11} Δ_{22} \\cdots Δ_{nn})^{1/n}.
```
See [`geometric_mean`](@ref).
"""
function det_rootn(model, X, t)
n = LinAlg.checksquare(X)
# Setup variables
Δ = Matrix{JuMP.AffExpr}(n, n)
for j in 1:n, i in j:n
Δ[i, j] = @variable model
end
z = zero(JuMP.AffExpr)
for j in 1:n, i in 1:j-1
Δ[i, j] = z
end
d = diag(Δ)
D = diagm(d)
@SDconstraint model [X Δ
Δ' D] ⪰ 0
# t^n <= (Z11*Z22*...*Znn)
geometric_mean(model, d, t)
end
using MathOptInterfaceMosek
solver = MosekSolver(QUIET=true)
using CSDP
solver = CSDPSolver(printlevel=0)
"""
lownerjohn_inner(A, b)
The inner ellipsoidal approximation to a polytope
```math
S = \\{\\, x \\in R^n \\mid Ax \\leq b \\,\\}
```
maximizes the volume of the inscribed ellipsoid,
```math
\\{\\, x \\mid x = Cu + d, \\lVert u \\lVert_2 \\leq 1 \\,\\}.
```
The volume is proportional to det(C)^(1/n), so the
problem can be solved as
```math
\\begin{align*}
\\max \\quad & t\\\\
\\text{subject to} \\quad & t \\leq \\det(C)^{1/n}\\\\
& \\lVert Ca_i \\lVert_2 \\leq b_i - a_i^T d, i=1,\\ldots,m\\\\
& C \\text{ is PSD}
\\end{align*}
```
which is equivalent to a mixed conic quadratic and semidefinite
programming problem.
"""
function lownerjohn_inner(A, b)
m, n = size(A)
model = Model(solver=solver)
# Setup variables
@variable model t >= 0
@variable model C[1:n,1:n] PSD
@variable model d[1:n]
# (b-Ad, AC) generate cones
T = b - A*d
X = A * C
for i in 1:m
@constraint model [T[i]; X[i,:]] in MOI.SecondOrderCone(1+n)
end
# t <= det(C)^{1/n}
det_rootn(model, C, t)
# Objective: Maximize t
@objective model Max t
solve(model)
@show JuMP.primalstatus(model)
JuMP.resultvalue.(C), JuMP.resultvalue.(d)
end
using Plots
plotlyjs()
#Vertices of a pentagon in 2D
x = [0., 1., 5.5, 7., 7., 3.]
y = [0., 3., 4.5, 4., 1., -2.]
nVerts = length(x)
#The hyperplane representation of the same polytope
A = Matrix{Float64}(nVerts, 2)
b = Vector{Float64}(nVerts)
for i in 1:nVerts
j = i == 1 ? nVerts : i-1
A[i, 1] = y[j] - y[i]
A[i, 2] = x[i] - x[j]
b[i] = A[i,1] * x[i] + A[i,2] * y[i]
end
plot(x, y, seriestype=:shape)
#Po, co = lownerjohn_outer(p)
#Ci, di = lownerjohn_inner(A, b)
Ci, di = lownerjohn_inner(A, b)
Qi = inv(Ci*Ci)
ci = di
θ = linspace(0, 2π, 100)
function scaling(θ)
xy = [cos(θ), sin(θ)]
sqrt(xy' * Qi * xy)
end
scl = scaling.(θ)
Δx = cos.(θ) ./ scl
Δy = sin.(θ) ./ scl
x = ci[1] .+ Δx
y = ci[2] .+ Δy
plot!(x, y, seriestype=:shape, fillalpha=0.6)
"""
The outer ellipsoidal approximation to a polytope given
as the convex hull of a set of points
S = conv{ x1, x2, ... , xm }
minimizes the volume of the enclosing ellipsoid,
{ x | || P*x-c ||_2 <= 1 }
The volume is proportional to det(P)^{-1/n}, so the problem can
be solved as
maximize t
subject to t <= det(P)^(1/n)
|| P*xi - c ||_2 <= 1, i=1,...,m
P is PSD.
"""
def lownerjohn_outer(x):
with Model("lownerjohn_outer") as M:
M.setLogHandler(sys.stdout)
m, n = len(x), len(x[0])
# Setup variables
t = M.variable("t", 1, Domain.greaterThan(0.0))
P = M.variable("P", Domain.inPSDCone(n))
c = M.variable("c", n, Domain.unbounded())
# (1, Px-c) in cone
M.constraint("qc",
Expr.hstack(Expr.ones(m),
Expr.sub(Expr.mul(x, P),
Var.reshape(Var.repeat(c, m), [m, n])
)
),
Domain.inQCone())
# t <= det(P)^{1/n}
#det_rootn(M, P, t)
# Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)
M.solve()
P, c = P.level(), c.level()
return ([P[i:i + n] for i in range(0, n * n, n)], c)
#Visualization
try:
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as patches
#Polygon
fig = plt.figure()
ax = fig.add_subplot(111)
ax.add_patch(patches.Polygon(p, fill=False, color="red"))
#The inner ellipse
theta = np.linspace(0, 2 * np.pi, 100)co
x = Ci[0][0] * np.cos(theta) + Ci[0][1] * np.sin(theta) + di[0]
y = Ci[1][0] * np.cos(theta) + Ci[1][1] * np.sin(theta) + di[1]
ax.plot(x, y)
#The outer ellipse
x, y = np.meshgrid(np.arange(-1.0, 8.0, 0.025),
np.arange(-3.0, 6.5, 0.025))
ax.contour(x, y,
(Po[0][0] * x + Po[0][1] * y - co[0])**2 + (Po[1][0] * x + Po[1][1] * y - co[1])**2, [1])
ax.autoscale_view()
ax.xaxis.set_visible(False)
ax.yaxis.set_visible(False)
fig.savefig('ellipsoid.png')
except:
print("Inner:")
print(" C = ", Ci)
print(" d = ", di)
print("Outer:")
print(" P = ", Po)
print(" c = ", co)