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.
Header file
/**
Bojan Nikolic <bojan@bnikolic.co.uk>
Initial version 2010
This file is part of BNMin1 and is licensed under GNU General
Public License version 2
\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
This file is part of BNMin1 and is licensed under GNU General
Public License version 2
\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;
}
}