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 
#include 
#include 
#include 

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

Is it effecient

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

Comments on code:

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 
#include 
#include 
#include 

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 
#include 

#include 
#include      // 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;

• Prefer to prefix types/objects with their namespace

-
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.

• 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

StackExchange Code Review Q#9554, answer score: 4