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.

_images/ellipses_polygon.png

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 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());

  // quadratic cones
  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 [BenTalN01] 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);
  }
}