Minimum volume ellipsoid covering a finite set

% Section 8.4.1, Boyd & Vandenberghe "Convex Optimization"
% Almir Mutapcic - 10/05
% (a figure is generated)
%
% Given a finite set of points x_i in R^2, we find the minimum volume
% ellipsoid (described by matrix A and vector b) that covers all of
% the points by solving the optimization problem:
%
%           maximize     log det A
%           subject to   || A x_i + b || <= 1   for all i
%
% CVX cannot yet handle the logdet function, but this problem can be
% represented in an equivalent way as follows:
%
%           maximize     det(A)^(1/n)
%           subject to   || A x_i + b || <= 1   for all i
%
% The expression det(A)^(1/n) is SDP-representable, and is implemented
% by the MATLAB function det_rootn().

% Generate data
x = [ 0.55  0.0;
      0.25  0.35
     -0.2   0.2
     -0.25 -0.1
     -0.0  -0.3
      0.4  -0.2 ]';
[n,m] = size(x);

% Create and solve the model
cvx_begin
    variable A(n,n) symmetric
    variable b(n)
    maximize( det_rootn( A ) )
    subject to
        norms( A * x + b * ones( 1, m ), 2 ) <= 1;
cvx_end

% Plot the results
clf
noangles = 200;
angles   = linspace( 0, 2 * pi, noangles );
ellipse  = A \ [ cos(angles) - b(1) ; sin(angles) - b(2) ];
plot( x(1,:), x(2,:), 'ro', ellipse(1,:), ellipse(2,:), 'b-' );
axis off
 
Calling sedumi: 31 variables, 9 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 9, order n = 19, dim = 39, blocks = 9
nnz(A) = 43 + 0, nnz(ADA) = 57, nnz(L) = 33
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            9.57E+00 0.000
  1 :   4.42E-01 2.33E+00 0.000 0.2435 0.9000 0.9000   2.99  1  1  1.0E+00
  2 :   1.50E+00 7.83E-01 0.000 0.3359 0.9000 0.9000   0.76  1  1  4.4E-01
  3 :   2.51E+00 1.76E-01 0.000 0.2253 0.9000 0.9000   0.42  1  1  1.1E-01
  4 :   2.68E+00 1.74E-02 0.000 0.0985 0.9900 0.9900   0.91  1  1  1.2E-02
  5 :   2.68E+00 3.35E-03 0.000 0.1930 0.9000 0.9000   1.00  1  1  2.4E-03
  6 :   2.68E+00 1.39E-04 0.259 0.0415 0.9900 0.9902   1.00  1  1  1.4E-04
  7 :   2.68E+00 4.69E-06 0.000 0.0337 0.6744 0.9000   1.00  1  1  3.9E-05
  8 :   2.68E+00 9.94E-08 0.000 0.0212 0.9900 0.9900   1.00  1  1  8.2E-07
  9 :   2.68E+00 3.71E-09 0.184 0.0373 0.9900 0.9900   1.00  1  1  3.1E-08
 10 :   2.68E+00 7.60E-10 0.000 0.2047 0.9000 0.9000   1.00  1  1  6.3E-09

iter seconds digits       c*x               b*y
 10      0.0   Inf  2.6839853844e+00  2.6839853874e+00
|Ax-b| =   8.6e-09, [Ay-c]_+ =   0.0E+00, |x|=  2.6e+00, |y|=  6.1e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    5.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 1,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.1766.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2.68399