Largely inspired from tools/examples/fusion/python/lownerjohn_ellipsoid.py in Mosek

### Purpose:¶

Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope.

### References:¶

•  Ben-Tal, Aharon, and Arkadi Nemirovski. Lectures on modern convex optimization: analysis, algorithms, and engineering applications. Society for Industrial and Applied Mathematics, 2001.
•  "MOSEK modeling manual", 2013
In [ ]:
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

In [ ]:
using JuMP # Needs branch jump/moi
using MathOptInterface
const MOI = MathOptInterface

In [ ]:
"""
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 )

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
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

# set the last m elements equal to t

# The type of x could be Vector{AffExpr} in which case
# x = [x; x_padding] would change the type of x
else
y = x
end

z = rec(y)
@constraint model sqrt(2)^l * t == z
end

In [ ]:
"""
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

In [ ]:
using MathOptInterfaceMosek
solver = MosekSolver(QUIET=true)

In [ ]:
using CSDP
solver = CSDPSolver(printlevel=0)

In [ ]:
"""
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*}
\\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]

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

In [ ]:
using Plots
plotlyjs()

In [ ]:
#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)

In [ ]:
Ci, di = lownerjohn_inner(A, b)
Qi = inv(Ci*Ci)
ci = di

In [ ]:
θ = linspace(0, 2π, 100)
function scaling(θ)
xy = [cos(θ), sin(θ)]
sqrt(xy' * Qi * xy)
end
scl = scaling.(θ)
Δx = cos.(θ) ./ scl
Δy = sin.(θ) ./ scl
x = ci .+ Δx
y = ci .+ Δy
plot!(x, y, seriestype=:shape, fillalpha=0.6)

In [ ]:
"""
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)

# 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)

In [ ]:
#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()
#The inner ellipse
theta = np.linspace(0, 2 * np.pi, 100)co
x = Ci * np.cos(theta) + Ci * np.sin(theta) + di
y = Ci * np.cos(theta) + Ci * np.sin(theta) + di
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 * x + Po * y - co)**2 + (Po * x + Po * y - co)**2, )
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)