C++ implementation of Khachiyan algorithm for the minimum enclosing (or covering) ellipsoid

The problem in this case consists of computing the hyper-ellipsoid (potentially in a space of high-dimensionality) such that it has the minimum volume but still encloses all of the points in a supplied set. The ellipsoid is specified by it's centre and a matrix describing it's axes.

There are many practical uses for such an algorithm -- I developed this C++ implementation as a part of a set of minimisation/inference routines. The full files will released as part of forthcoming version of BNMin1, but in the mean time here is a fully working version containing all of the essential functionality.

This code is licensed to visitors of this site under GNU Public License Version 2.

```/**
Bojan Nikolic <bojan@bnikolic.co.uk>
Initial version 2010

\file ellipsoids.hxx

Computation and use of ellipsoids releated to sets of
points. References are to Todd and Yildirim, "On Khachiyan's
Algorithm for the Computation of Minimum Volume Enclosing
Ellipsoids", 2005

*/
#ifndef _BNMIN1_SETS_ELIIPSOIDS_HXX__
#define _BNMIN1_SETS_ELIIPSOIDS_HXX__

#include <set>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>

namespace Minim {

namespace ublas=boost::numeric::ublas;

/**  Lift to a higher dimensional space that so that points are
centrally symmetric.

The input data are arranged in columns of A

c.f. Eq 9
*/
void Lift(const ublas::matrix<double> &A,
ublas::matrix<double> &Ap);

/** Compute the result of operator Lambda

c.f. eq. 13

\param Lambdap The result is stored here
*/
void KaLambda(const ublas::matrix<double> &Ap,
const ublas::vector<double> &p,
ublas::matrix<double> &Lambdap);

/** Calculate the inverse of a result of lambda operator
*/
void InvertLP(const ublas::matrix<double> &Lambdap,
ublas::matrix<double> &LpInv);

/** Iterate the Khachiyan algorithm one step

\param Ap is the lifted data set

\param p is the current guess, and the iterated guess is stored
here

\returns Error estimate
*/
double KhachiyanIter(const ublas::matrix<double> &Ap,
ublas::vector<double> &p);

/** Invert from the dual problem back to real space

\param Q The ellipsiod matrix
\param c The centre of the ellipse
*/
void KaInvertDual(const ublas::matrix<double> &A,
const ublas::vector<double> &p,
ublas::matrix<double> &Q,
ublas::vector<double> &c
);

/** \param Compute minimum containing ellipsoid using the Khachiyan
algorithm

\param A Input point set as columns of the matrix

\param eps Tolerance of the estimated ellipsoid: it's volume
will be approx (1+eps) of the true smallest ellipsoid

\param maxiter Maximum number of iterations to make before
stopping the algorithm

\param Q The matrix describing the axes of the resulting
ellipsoid

\param c The centre of the resulting ellipsoid
*/
double KhachiyanAlgo(const ublas::matrix<double> &A,
double eps,
size_t maxiter,
ublas::matrix<double> &Q,
ublas::vector<double> &c);

/** \brief Convenience class to store the results of minimum
covering ellipsoid calculation
*/
struct KhachiyanEllipsoid
{
ublas::matrix<double> Q;
ublas::vector<double> c;
};

}

#endif
```

Implementation:

```/**
Bojan Nikolic <bojan@bnikolic.co.uk>
Initial version 2010

\file ellipsoids.cxx

Computation and use of ellipsoids releated to sets of points
*/

#include <boost/numeric/ublas/triangular.hpp>
#include <boost/numeric/ublas/banded.hpp>
#include <boost/numeric/ublas/lu.hpp>

#include <iostream>
#include <boost/numeric/ublas/io.hpp>

#include "ellipsoids.hxx"
#include "../mcpoint.hxx"
#include "../bnmin_main.hxx"

namespace Minim {

template<class T>
bool InvertMatrix(const ublas::matrix<T> &input,
ublas::matrix<T> &inverse)
{
using namespace boost::numeric::ublas;
typedef permutation_matrix<std::size_t> pmatrix;
matrix<T> A(input);
pmatrix pm(A.size1());
int res = lu_factorize(A,pm);
if( res != 0 ) return false;
inverse.assign(ublas::identity_matrix<T>(A.size1()));
lu_substitute(A, pm, inverse);
return true;
}

void InvertLP(const ublas::matrix<double> &Lambdap,
ublas::matrix<double> &LpInv)
{
bool res=InvertMatrix(Lambdap, LpInv);
if (not res)
{
// throw an error of your choice here!
// throw MatrixErr("Could not invert matrix: ",
//                 Lambdap);
}
}

void Lift(const ublas::matrix<double> &A,
ublas::matrix<double> &Ap)
{
Ap.resize(A.size1()+1,
A.size2());
ublas::matrix_range<ublas::matrix<double> >
sub(Ap,
ublas::range(0, A.size1()),
ublas::range(0, A.size2()));
sub.assign(A);
ublas::row(Ap, Ap.size1()-1)=ublas::scalar_vector<double>(A.size2(),1.0);

}

void genDiag(const ublas::vector<double> &p,
ublas::matrix<double> &res)
{
res.assign(ublas::zero_matrix<double>(p.size(),
p.size()));
for(size_t i=0; i<p.size(); ++i)
{
res(i,i)=p(i);
}
}

void KaLambda(const ublas::matrix<double> &Ap,
const ublas::vector<double> &p,
ublas::matrix<double> &Lambdap)
{

ublas::matrix<double> dp(p.size(), p.size());
genDiag(p, dp);

dp=ublas::prod(dp, ublas::trans(Ap));
Lambdap=ublas::prod(Ap,
dp);
}

double KhachiyanIter(const ublas::matrix<double> &Ap,
ublas::vector<double> &p)
{
/// Dimensionality of the problem
const size_t d=Ap.size1()-1;

ublas::matrix<double> Lp;
ublas::matrix<double> M;
KaLambda(Ap, p, Lp);
ublas::matrix<double> ILp(Lp.size1(), Lp.size2());
InvertLP(Lp, ILp);
M=ublas::prod(ILp, Ap);
M=ublas::prod(ublas::trans(Ap), M);

double maxval=0;
size_t maxi=0;
for(size_t i=0; i<M.size1(); ++i)
{
if (M(i,i) > maxval)
{
maxval=M(i,i);
maxi=i;
}
}
const double step_size=(maxval -d - 1)/((d+1)*(maxval-1));
ublas::vector<double> newp=p*(1-step_size);
newp(maxi) += step_size;

const double err= ublas::norm_2(newp-p);
p=newp;
return err;

}

void KaInvertDual(const ublas::matrix<double> &A,
const ublas::vector<double> &p,
ublas::matrix<double> &Q,
ublas::vector<double> &c
)
{
const size_t d=A.size1();
ublas::matrix<double> dp(p.size(), p.size());
genDiag(p, dp);

ublas::matrix<double> PN=ublas::prod(dp, ublas::trans(A));
PN=ublas::prod(A, PN);

ublas::vector<double> M2=ublas::prod(A, p);
ublas::matrix<double> M3=ublas::outer_prod(M2, M2);

ublas::matrix<double> invert(PN.size1(), PN.size2());
InvertLP(PN- M3, invert);

Q.assign( 1.0/d *invert);
c=ublas::prod(A, p);

}

double KhachiyanAlgo(const ublas::matrix<double> &A,
double eps,
size_t maxiter,
ublas::matrix<double> &Q,
ublas::vector<double> &c)
{
ublas::vector<double> p=ublas::scalar_vector<double>(A.size2(), 1.0)*(1.0/A.size2());

ublas::matrix<double> Ap;
Minim::Lift(A, Ap);

double ceps=eps*2;
for (size_t i=0;  i<maxiter && ceps>eps; ++i)
{
ceps=KhachiyanIter(Ap, p);
}

KaInvertDual(A, p, Q, c);

return ceps;

}

double KhachiyanAlgo(const std::set<MCPoint> &ss,
double eps,
size_t maxiter,
KhachiyanEllipsoid &res)
{
const size_t d=ss.begin()->p.size();
ublas::matrix<double> A(d,
ss.size());

size_t j=0;
for (std::set<MCPoint>::const_iterator i=ss.begin();
i != ss.end();
++i)
{
for(size_t k=0; k <d ;++k)
A(k,j)=i->p[k];
++j;
}

ublas::matrix<double> Q(d,d);
ublas::vector<double> c(d);

const double ceps=KhachiyanAlgo(A, eps, maxiter,
Q, c);
res.Q=Q;
res.c=c;
return ceps;
}

}
```