11.6 Inner and outer Löwner-John Ellipsoids¶
In this section we show how to compute the Löwner-John inner and outer ellipsoidal approximations of a polytope. They are defined as, respectively, the largest volume ellipsoid contained inside the polytope and the smallest volume ellipsoid containing the polytope, as seen in Fig. 11.7.

Fig. 11.7 The inner and outer Löwner-John ellipse of a polygon.¶
For further mathematical details, such as uniqueness of the two ellipsoids, consult [BenTalN01]. Our solution is a mix of conic and semidefinite programming. Among other things, in Sec. 11.6.3 (Bound on the Determinant Root) we show how to implement bounds involving the determinant of a PSD matrix.
11.6.1 Inner Löwner-John Ellipsoids¶
Suppose we have a polytope given by an h-representation
and we wish to find the inscribed ellipsoid with maximal volume. It will be convenient to parametrize the ellipsoid as an affine transformation of the standard disk:
Every non-degenerate ellipsoid has a parametrization such that
where
def lownerjohn_inner(A, b):
with Model("lownerjohn_inner") as M:
M.setLogHandler(sys.stdout)
m, n = len(A), len(A[0])
# Setup variables
t = M.variable("t", 1, Domain.greaterThan(0.0))
C = det_rootn(M, t, n)
d = M.variable("d", n, Domain.unbounded())
# (b-Ad, AC) generate cones
M.constraint("qc", Expr.hstack(b - A @ d, A @ C), Domain.inQCone())
# Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)
M.writeTask("lj-inner.ptf")
M.solve()
C, d = C.level(), d.level()
return ([C[i:i + n] for i in range(0, n * n, n)], d)
The only black box is the method det_rootn
which implements the constraint
11.6.2 Outer Löwner-John Ellipsoids¶
To compute the outer ellipsoidal approximation to a polytope, let us now start with a v-representation
of the polytope as a convex hull of a set of points. We are looking for an ellipsoid given by a quadratic inequality
whose volume is proportional to
and then directly into Fusion code:
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 = det_rootn(M, t, n)
c = M.variable("c", Domain.unbounded().withShape(1,n))
# (1, Px-c) in cone
M.constraint("qc", Expr.hstack(Expr.ones(m), x @ P - Expr.repeat(c,m,0)), Domain.inQCone())
# Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t)
M.writeTask("lj-outer.ptf")
M.solve()
P, c = P.level(), c.level()
return ([P[i:i + n] for i in range(0, n * n, n)], c)
11.6.3 Bound on the Determinant Root¶
It remains to show how to express the bounds on
A standard approach when working with the determinant of a PSD matrix is to consider a semidefinite cone
where
but this is exactly the geometric mean cone Domain.inPGeoMeanCone
. We obtain the following model:
def det_rootn(M, t, n):
# Setup variables
Y = M.variable(Domain.inPSDCone(2 * n))
# Setup Y = [X, Z; Z^T , diag(Z)]
X = Y[0:n, 0:n]
Z = Y[0:n, n:2*n]
DZ = Y[n:2*n, n:2*n]
# Z is lower-triangular
M.constraint(Z.pick([[i,j] for i in range(n) for j in range(i+1,n)]), Domain.equalsTo(0.0))
# DZ = Diag(Z)
M.constraint(DZ == Expr.mulElm(Z, Matrix.eye(n)))
# (Z11*Z22*...*Znn) >= t^n
M.constraint(Expr.vstack(DZ.diag(), t), Domain.inPGeoMeanCone())
# Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
return X