C++ matrix computation efficiency

Question

We are trying to optimize our C++ code and the we have the following matrix computation (using the Eigen library)

#include<Eigen/Dense>

int main(){

   MatrixXd P = MatrixXd::Random(30,30); // a  random double 30 x 30 matrix P
   MatrixXd M = MatrixXd::Random(30,30); // a  random double 30 x 30 matrix M
   Matrix<double, 30, 30> I; 
   I.setIdentity(); // I is an 30 x 30 identity matirx

   P = (I-M)*P

   return 0;

   }

Where they are all nxn matrices and I is the identity matrix. We found the rewriting the above matrix computation

   P= (I- M)*P

as

   P = P-M*P

results in ~4-8x speed up in a Linux ubuntu system using a gcc 6.2 compiler. I realize the fact that the compiler might not know anything about identity matrix and the fact I*P = P, but still cannot wrap my head around what makes the efficiency improve so much. Anyone knows the possible reasons that make such substantial improvements?

Answer 1

First of all, I.identity(); does not exist. What you want is either I.setIdentity(), or P = (MatrixXd::Identity(30,30)-M)*P. If you use the first option, Eigen will definitely need to do a full 30x30 subtraction of I and M (it would be very hard for a compiler to see the equivalence to your second expression). Overall, this will result in two temporaries (one for the difference, one for the product).

If you actually used I.Identity() you are calling a static function like a member function, and your compiler should at least warn you about that. This will not actually modify I and you end up with uninitialized values in I, which likely will include some NaN or denormal values, both can be bad for floating point performance. And of course your result would be wrong.

Overall, I'd think the easiest ways to write your equation are

P -= M*P;

or

MatrixXd Pnew = P - M*P;