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:

  • [1] Ben-Tal, Aharon, and Arkadi Nemirovski. Lectures on modern convex optimization: analysis, algorithms, and engineering applications. Society for Industrial and Applied Mathematics, 2001.
  • [2] "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 [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
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*}
\\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
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[1] .+ Δx
y = ci[2] .+ Δ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[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)
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()
        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)