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

For further mathematical details, such as uniqueness of the two ellipsoids, consult . Our solution is a mix of conic quadratic 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

$\mathcal{P} = \{ x \in \real^n \mid Ax \leq b \}$

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:

$\mathcal{E} = \{ x \mid x = Cu + d,\ u\in\real^n,\ \| u \|_2 \leq 1 \}.$

Every non-degenerate ellipsoid has a parametrization such that $$C$$ is a positive definite symmetric $$n\times n$$ matrix. Now the volume of $$\mathcal{E}$$ is proportional to $$\mbox{det}(C)^{1/n}$$. The condition $$\mathcal{E}\subseteq\mathcal{P}$$ is equivalent to the inequality $$A(Cu+d)\leq b$$ for all $$u$$ with $$\|u\|_2\leq 1$$. After a short computation we obtain the formulation:

(11.22)$\begin{split}\begin{array}{lll} \maximize & t & \\ \st & t \leq \mbox{det}(C)^{1/n}, & \\ & (b-Ad)_i\geq \|(AC)_i\|_2, & i=1,\ldots,m,\\ & C \succeq 0, & \end{array}\end{split}$

where $$X_i$$ denotes the $$i$$-th row of the matrix $$X$$. This can easily be implemented using Fusion, where the sequence of conic inequalities can be realized at once by feeding in the matrices $$b-Ad$$ and $$AC$$.

Listing 11.11 Fusion implementation of model (11.22). Click here to download.
std::pair<std::shared_ptr<ndarray<double, 1>>, std::shared_ptr<ndarray<double, 1>>>
lownerjohn_inner
( std::shared_ptr<ndarray<double, 2>> A,
std::shared_ptr<ndarray<double, 1>> b)
{
Model::t M = new Model("lownerjohn_inner"); auto _M = finally([&]() { M->dispose(); });
int m = A->size(0);
int n = A->size(1);

// Setup variables
Variable::t t = M->variable("t", 1, Domain::greaterThan(0.0));
Variable::t C = det_rootn(M, t, n);
Variable::t d = M->variable("d", n, Domain::unbounded());

M->constraint(Expr::hstack(Expr::sub(b, Expr::mul(A, d)), Expr::mul(A, C)),
Domain::inQCone());

// Objective: Maximize t
M->objective(ObjectiveSense::Maximize, t);
M->solve();

return std::make_pair(C->level(), d->level());
}


The only black box is the method det_rootn which implements the constraint $$t\leq \mbox{det}(C)^{1/n}$$. It will be described in Sec. 11.6.3 (Bound on the Determinant Root).

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

$\mathcal{P} = \mbox{conv}\{ x_1, x_2, \ldots , x_m \} \subseteq \real^n,$

of the polytope as a convex hull of a set of points. We are looking for an ellipsoid given by a quadratic inequality

$\mathcal{E} = \{ x\in\real^n \mid \| Px-c \|_2 \leq 1 \},$

whose volume is proportional to $$\mbox{det}(P)^{-1/n}$$, so we are after maximizing $$\mbox{det}(P)^{1/n}$$. Again, there is always such a representation with a symmetric, positive definite matrix $$P$$. The inclusion conditions $$x_i\in\mathcal{E}$$ translate into a straightforward problem formulation:

(11.23)$\begin{split}\begin{array}{lll} \maximize & t &\\ \st & t \leq \mbox{det}(P)^{1/n}, &\\ & \|Px_i - c\|_2 \leq 1, &i=1,\ldots,m,\\ & P \succeq 0, & \end{array}\end{split}$

and then directly into Fusion code:

Listing 11.12 Fusion implementation of model (11.23). Click here to download.
std::pair<std::shared_ptr<ndarray<double, 1>>, std::shared_ptr<ndarray<double, 1>>>
lownerjohn_outer(std::shared_ptr<ndarray<double, 2>> x)
{
Model::t M = new Model("lownerjohn_outer");
int m = x->size(0);
int n = x->size(1);

// Setup variables
Variable::t t = M->variable("t", 1, Domain::greaterThan(0.0));
Variable::t P = det_rootn(M, t, n);
Variable::t c = M->variable("c", n, Domain::unbounded());

// (1, Px-c) \in Q
M->constraint(Expr::hstack(
Expr::ones(m), Expr::sub(Expr::mul(x, P),
Var::reshape(Var::repeat(c, m), new_array_ptr<int, 1>({m, n}))) ),
Domain::inQCone());

// Objective: Maximize t
M->objective(ObjectiveSense::Maximize, t);
M->solve();

return std::make_pair(P->level(), c->level());
}


## 11.6.3 Bound on the Determinant Root¶

It remains to show how to express the bounds on $$\mbox{det}(X)^{1/n}$$ for a symmetric positive definite $$n\times n$$ matrix $$X$$ using PSD and conic quadratic variables. We want to model the set

(11.24)$C = \lbrace (X, t) \in \PSD^n \times \real \mid t \leq \mbox{det}(X)^{1/n} \rbrace.$

A standard approach when working with the determinant of a PSD matrix is to consider a semidefinite cone

(11.25)$\begin{split}\left( {\begin{array}{cc}X & Z \\ Z^T & \mbox{Diag}(Z) \\ \end{array} } \right) \succeq 0\end{split}$

where $$Z$$ is a matrix of additional variables and where we intuitively identify $$\mbox{Diag}(Z)=\{\lambda_1,\ldots,\lambda_n\}$$ with the eigenvalues of $$X$$. With this in mind, we are left with expressing the constraint

(11.26)$t \leq (\lambda_1\cdot\ldots\cdot\lambda_n)^{1/n}.$

This is easy to implement recursively using rotated quadratic cones when $$n$$ is a power of $$2$$. In general it is convenient to express (11.26) as a composition of power cones, see or Modeling Cookbook.

Listing 11.13 Approaching the determinant, see (11.25). Click here to download.
Variable::t det_rootn(Model::t M, Variable::t t, int n)
{
// Setup variables
Variable::t Y = M->variable(Domain::inPSDCone(2 * n));

Variable::t X   = Y->slice(new_array_ptr<int,1>({0, 0}), new_array_ptr<int,1>({n, n}));
Variable::t Z   = Y->slice(new_array_ptr<int,1>({0, n}), new_array_ptr<int,1>({n, 2 * n}));
Variable::t DZ  = Y->slice(new_array_ptr<int,1>({n, n}), new_array_ptr<int,1>({2 * n, 2 * n}));

// Z is lower-triangular
std::shared_ptr<ndarray<int,2>> low_tri( new ndarray<int,2>( shape_t<2>(n*(n-1)/2, 2) ));
int k = 0;
for(int i = 0; i < n; i++)
for(int j = i+1; j < n; j++)
(*low_tri)(k,0) = i, (*low_tri)(k,1) = j, ++k;
M->constraint(Z->pick(low_tri), Domain::equalsTo(0.0));
// DZ = Diag(Z)
M->constraint(Expr::sub(DZ, Expr::mulElm(Z, Matrix::eye(n))), Domain::equalsTo(0.0));

// t^n <= (Z11*Z22*...*Znn)
geometric_mean(M, DZ->diag(), t);

// Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
return X;
}

Listing 11.14 Bounding the geometric mean, see (11.26). Click here to download.
void geometric_mean(Model::t M, Variable::t x, Variable::t t)
{
int n = (int) x->getSize();
if (n==1)
M->constraint(Expr::sub(t, x), Domain::lessThan(0.0));
else
{
Variable::t t2 = M->variable();
M->constraint(Var::hstack(t2, x->index(n-1), t), Domain::inPPowerCone(1-1.0/n));
geometric_mean(M, x->slice(0,n-1), t2);
}
}