# Implementation of Mahalanobis distances using Eigen

I'm trying to learn C++ with Eigen. Here's my attempt at computing Mahalanobis distances of a set of points x with respect to a sub-matrix xs.

The aim of the project is to turn an R code describing a statistical procedure in C++ (and in the process to learn a bit about numerical computing in c++). Some of these functions will have to be called literally 10,000 of times so the efficiency part is a bit important.

Is this an efficient implementation of what I'm trying to do in Eigen (the code does what I want)? Can it be improved? What are the beginner's mistake that I'm doing? (in the application I have in mind, x is an imported matrix, so 30 and 3 are only known at runtime)

#include <iostream>
#include <Eigen/Dense>
#include <fstream>
#include <Eigen/Cholesky>

using namespace Eigen;
using namespace std;

VectorXd mah(MatrixXd& x, MatrixXd& xs){
int ps = xs.cols(), ns=xs.rows();
RowVectorXd xs_mean=xs.colwise().sum()/ns;
MatrixXd xs_cen = (xs.rowwise()-xs_mean);
MatrixXd x_cen = (x.rowwise()-xs_mean);
MatrixXd w = xs_cen.transpose()*xs_cen/(ns-1);
MatrixXd b = MatrixXd::Identity(ps,ps);
w.ldlt().solveInPlace(b);
x_cen = (x_cen*b).cwiseProduct(x_cen);
return x_cen.rowwise().sum();
}

MatrixXd sub(MatrixXd&  x, VectorXi& s){
int p = x.cols();
MatrixXd xs(p+1,p);
for(int i=0;i<s.size();i++){
xs.row(i)=x.row(s(i));
}
return xs;
}

int main(){
MatrixXd x = MatrixXd::Random(30,3);
VectorXi s = (VectorXd::Random(p+1)*n).cast<int>();
MatrixXd xs = sub(x,s);
ofstream file("test1.txt");
if (file.is_open()){
file <<  x << '\n';
}
cout << "m" << endl <<  mah(x,xs) << endl;
cout << "s" << endl <<  s << endl;
return 0;
}


### Efficiency

Seem fine ! It's a good idea to compute all input samples at once.

If covariance matrix remains small (dim <= 4), you may use inverse().

I was wondering whether the inverse is needed at all with your approach. I've tried to replace :

MatrixXd b = MatrixXd::Identity(ps,ps);
w.ldlt().solveInPlace(b);
x_cen = (x_cen*b).cwiseProduct(x_cen);


by

MatrixXd b = w.ldlt().solve(x_cen.transpose());
x_cen = b.transpose().cwiseProduct(x_cen);


But crude benchmarking showed +20% of process time ! You may want to check with real data, since it highly depends on matrix dimensions.

To trigger efficiency boost :

• Use floats instead of double (for SIMD)
• Enable the right flags (-fopenmp, -msse2 ...)

### General review

Post compiling code ! In test_Mahalahobis, n and p are not defined.

It is a better design to compute covariance matrix in a separate function. You may wish to obtain the distance in regard to a covariance matrix known once for all.

Documentation : if you don't do it for you, do it for the reviewers (and vice-versa). For instance, says that there is one sample by row. And legible code starts with good names :

• mah : meh !
• sub : subset is less ambiguous, and not very long to type.

### Cosmetic points

Replace

 RowVectorXd xs_mean=xs.colwise().sum()/ns;   //---
RowVectorXd xs_mean=xs.colwise().mean();     //+++ Explicit and shorter

x_cen = (x_cen*b).cwiseProduct(x_cen);       //---Deprecated
x_cen = (x_cen*b).array() * x_cen.array() ;  //+++Array operations are element wise

MatrixXd xs(p+1,p);        //--- Implicit assumption
MatrixXd xs(s.size(),p);   //+++

• thanks for all these comments. I would up-vote more if i could. – user189035 Feb 29 '12 at 22:35

### Is it effecient

Imposable to tell as we have no idea how stuff in Eigen is implemented.

There should be some logic in the ordering of includes. So you can easily find ones that have been included or spot missing ones (or know where to put knew ones).

#include <iostream>
#include <Eigen/Dense>
#include <fstream>
#include <Eigen/Cholesky>


This set seems to be completely random. Personally (though the exact orders does not matter and you can pick your own scheme) I put them in groups (so in this case there would be an Eigen group and a standard lib group). Then I order the groups most specific to most general. Within the groups I am less pedantic but usually it is alphabetical unless there is some good reason to use another scheme.

#include <Eigen/Cholesky>
#include <Eigen/Dense>

#include <iostream>
#include <fstream>     // I usually group all the stream stuff together
// Then all containers then all the others in the standard lib.
//


Avoid the using-directive. In non trivial programs it causes more problems than it is worth. For short programs the extra cost of prefixing objects (in terms of typing) is negligible so prefix your types and objects (don't be lazy). (I always use std::cout not cout). If you feel this is too much work then selectively bring types/objects into the current namespace with a using statement.

using namespace Eigen;
using namespace std;

1. Prefer to prefix types/objects with their namespace
2. If you must selectively bring types/objects into current namespace

using std::cout; using std::endl;

• If you must do this then scope the using to minimize polution.
3. Try to not use using namespace X;

Does x change in the following code? If not then pass by const reference to indicate that it is not modified.

VectorXd mah(MatrixXd& x){


One variable per line:

 int p = x.cols(), n=x.rows();


This is a nasty trick to temporarily comment out blocks.

/*  ofstream file("test1.txt");
if (file.is_open()){
file <<  m << '\n';
}
/**/


It is so easily broken that I would avoid it.
If you want to comment out a block then use #if 0

// To uncomment block just replace 0 with 1.
#if 0
ofstream file("test1.txt");
if (file.is_open()){
file <<  m << '\n';
}
#endif

• thank you! One question: to declar x as constant, i just do: 'const MatrixXd x = MatrixXd::Random(30,3);' or should i also modify something in the function? Best – user189035 Feb 29 '12 at 14:28
• No. In the method do this: VectorXd mah(MatrixXd const& x){ Here you are telling people that mah() will not alter x. – Martin York Feb 29 '12 at 14:31
• thanks it was not initially clear. By the way: for clarity i'm leaving the code above unchanged...but i'm following your guidelines in the version of the code on my computer. Thanks again for taking the time. – user189035 Feb 29 '12 at 14:34
• +1, but using namespace X; is a using-directive, not a using-declaration. – Anton Golov Feb 29 '12 at 15:45