Strange behaviours in porting code from Eigen2 to Eigen3

Question

I'm considering to use this library to perform spectral clustering in my research project.

But, to do so, I need to port it from Eigen2 to Eigen3 (which is what I use in my code).

There's a portion of code that is causing me some troubles.

This is for Eigen2:

double Evrot::evqual(Eigen::MatrixXd& X) {
    // take the square of all entries and find max of each row
    Eigen::MatrixXd X2 = X.cwise().pow(2);
    Eigen::VectorXd max_values = X2.rowwise().maxCoeff();

    // compute cost
    for (int i=0; i<mNumData; i++ ) {
        X2.row(i) = X2.row(i) / max_values[i];
    }
    double J = 1.0 - (X2.sum()/mNumData -1.0)/mNumDims;
    if( DEBUG )
        std::cout << "Computed quality = "<< J << std::endl;

    return J;
}

as explained here, Eigen3 replaces .cwise() with the slightly different .array() functionality.

So, I wrote:

double Evrot::evqual(Eigen::MatrixXd& X) {
    // take the square of all entries and find max of each row
    Eigen::MatrixXd X2 = X.array().pow(2);

    Eigen::VectorXd max_values = X2.rowwise().maxCoeff();

    // compute cost
    for (int i=0; i<mNumData; i++ ) {
        X2.row(i) = X2.row(i) / max_values[i];
    }
    double J = 1.0 - (X2.sum()/mNumData -1.0)/mNumDims;
    if( DEBUG )
        std::cout << "Computed quality = "<< J << std::endl;

    return J;
}

and I got no compiler errors.

But, if I give to the two programs the same input (and check that they're actually getting consistent inputs), in the first case I get numbers and in the second only NANs.

My idea is that this is caused by the fact that max_values is badly computed and then using this vector in a division causes all the NANs. But I have no clue on how to fix that.

Can, please, someone explain me how to solve this problem?

Thanks!

Answer 1

Have you checked when the values start to diverge ? Are you sure there is no empty rows or that X^2 do not underflow. Anyways, you should had a guard before dividing by max_values[i]. Moreover, to avoid underflow in squaring you could rewrite it like that:

VectorXd max_values = X.array().abs().rowwise().maxCoeff();
double sum = 0;
for (int i=0; i<mNumData; i++ ) {
  if(max_values[i]>0)
    sum += (X.row(i)/max_values[i]).squaredNorm();
}
double J = 1.0 - (sum/mNumData -1.0)/mNumDims;

This will work even if X.abs().maxCoeff()==1e-170 whereas your code will underflow and produces NaN. Of course, if you are in a such a case, maybe you should check your inputs first as you are already on dangerous side regarding numerical issues.