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;
  }


}