10.4 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. 10.5.

_images/ellipses_polygon.png

Fig. 10.5 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 Section 10.4.3 we show how to implement bounds involving the determinant of a PSD matrix.

10.4.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:

(1)\[\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 10.7 Fusion implementation of model (1). Click here to download.
        function [C, d] = lownerjohn_inner(A, b)
        
            import mosek.fusion.*;
        
            M = Model('lownerjohn_inner');
        
            [m, n] = size(A);  

            % Setup variables
            t = M.variable('t', 1, Domain.greaterThan(0.0));
            C = M.variable('C', Domain.inPSDCone(n));
            d = M.variable('d', n, Domain.unbounded());
        
            % (b-Ad, AC) generate cones 
            M.constraint('qc', Expr.hstack(Expr.sub(b, Expr.mul(A,d)), Expr.mul(A,C.transpose())), Domain.inQCone());
            
            % t <= det(C)^{1/n}
            model_utils.det_rootn(M, C, t);
            
            % Objective: Maximize t
            M.objective(ObjectiveSense.Maximize, t);
            
            M.solve();
            
            C = reshape(C.level(), n, n);
            d = reshape(d.level(), n, 1);
    
            M.dispose();
        end

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 Section 10.4.3.

10.4.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:

(2)\[\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 10.8 Fusion implementation of model (2). Click here to download.
        function [P, c] = lownerjohn_outer(x)
    
            import mosek.fusion.*;
    
            M = Model('lownerjohn_outer');
            
            [m, n] = size(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 Q 
            M.constraint('qc', ...
                         Expr.hstack(Expr.ones(m), ...
                                     Expr.sub(Expr.mul(x,P.transpose()), ...
                                              Var.reshape(Var.repeat(c,m), [m,n]))), ...
                         Domain.inQCone());

            % t <= det(P)^{1/n}
            model_utils.det_rootn(M, P, t);
        
            % Objective: Maximize t
            M.objective(ObjectiveSense.Maximize, t);
            M.solve();
    
            P = reshape(P.level(), n, n);
            c = reshape(c.level(), n, 1);
            
            M.dispose();
        end

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

(3)\[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

(4)\[\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

(5)\[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\); otherwise we need to round \(n\) up to the nearest power of \(2\) as in Listing 10.10. For example, \(t\leq (\lambda_1\lambda_2\lambda_3\lambda_4)^{1/4}\) is equivalent to

\[\lambda_1\lambda_2\geq y_1^2,\ \lambda_3\lambda_4\geq y_2^2,\ y_1y_2\geq t^2\]

while \(t\leq (\lambda_1\lambda_2\lambda_3)^{1/3}\) can be achieved by writing \(t\leq (t\lambda_1\lambda_2\lambda_3)^{1/4}\).

For further details and proofs see [BenTalN01] or [MOSEKApS12].

Listing 10.9 Approaching the determinant, see (4). Click here to download.
        function [] = det_rootn(M, X, t)

            import mosek.fusion.*;

            n = X.shape.dim(1);

            % Setup variables
            Y = M.variable(Domain.inPSDCone(2*n));
        
            % Setup Y = [X Z; Z^T diag(Z)] 
            Y11 = Y.slice([1,   1],   [n+1,   n+1]);
            Y21 = Y.slice([n+1, 1],   [2*n+1, n+1]);
            Y22 = Y.slice([n+1, n+1], [2*n+1, 2*n+1]);

            S = Matrix.sparse(n, n, 1:n, 1:n, ones(1,n));
            M.constraint( Expr.sub(Y21.diag(), Y22.diag()), Domain.equalsTo(0.0) );
            M.constraint( Expr.sub(X, Y11), Domain.equalsTo(0.0) );
            
            % t^n <= (Z11*Z22*...*Znn)
            model_utils.geometric_mean(M, Y22.diag(), t);
        end
Listing 10.10 Bounding the geometric mean, see (5) Click here to download.
        function geometric_mean(M, x, t)
            import mosek.fusion.*;
            
            n = x.size();
            l = ceil(log2(n));
            m = 2^l - n;
            
            if (m == 0)
                x0 = x;
            else    
                x0 = Variable.vstack(x, M.variable(m, Domain.greaterThan(0.0)));
            end
            
            z = x0;
            
            for i=1:l-1,
                xi = M.variable(2^(l-i), Domain.greaterThan(0.0));
                
                for k=1:2^(l-i),        
                    M.constraint(Variable.hstack( z.index(2*k-1),z.index(2*k),xi.index(k)),...
                                 Domain.inRotatedQCone());
                end
                z = xi;
            end
        
            t0 = M.variable(1, Domain.greaterThan(0.0));
            M.constraint(Var.vstack(z, t0), Domain.inRotatedQCone());
            
            M.constraint(Expr.sub(Expr.mul(2^(0.5*l),t),t0), Domain.equalsTo(0.0));
            for i=2^l-m+1:2^l
                M.constraint(Expr.sub(x0.index(i), t), Domain.equalsTo(0.0));
            end
            
            t0 = M.variable(1, Domain.greaterThan(0.0));
            M.constraint(Var.vstack(z, t0), Domain.inRotatedQCone());
            
            M.constraint(Expr.sub(Expr.mul(2^(0.5*l),t),t0), Domain.equalsTo(0.0));
            for i=2^l-m+1:2^l
                M.constraint(Expr.sub(x0.index(i), t), Domain.equalsTo(0.0));
            end
        
        end