 <p>Suggested by remi, I implemented the same matrix multiplication using Eige. Here it is:</p> <pre><code>const int descrSize = 128; MatrixXi a(nframes, descrSize); MatrixXi b(vocabulary_size, descrSize); MatrixXi ab(nframes, vocabulary_size); unsigned char* dataPtr = DATAdescr; for (int i=0; i&lt;nframes; ++i) { for (int j=0; j&lt;descrSize; ++j) { a(i,j)=(int)*dataPtr++; } } unsigned char* vocPtr = vocabulary; for (int i=0; i&lt;vocabulary_size; ++i) { for (int j=0; j&lt;descrSize; ++j) { b(i,j)=(int)*vocPtr ++; } } ab = a*b.transpose(); a.cwiseProduct(a); b.cwiseProduct(b); MatrixXi aa = a.rowwise().sum(); MatrixXi bb = b.rowwise().sum(); MatrixXi d = (aa.replicate(1,vocabulary_size) + bb.transpose().replicate(nframes,1) - 2*ab).cwiseAbs2(); </code></pre> <p>The key line is the line that says </p> <pre><code>ab = a*b.transpose(); </code></pre> <p>vocabulary an DATAdescr are arrays of unsigned char. DATAdescr is 2782x128 and vocabulary is 4000x128. I saw at implementation that I can use Map, but I failed at first to use it. The initial loops for assigment are 0.001 cost, so this is not a bottleneck. The whole process is about 1.23 s</p> <p>The same implementation in matlab (0.05s.) is:</p> <pre><code>aa=sum(a.*a,2); bb=sum(b.*b,2); ab=a*b'; d = sqrt(abs(repmat(aa,[1 size(bb,1)]) + repmat(bb',[size(aa,1) 1]) - 2*ab)); </code></pre> <p>Thanks remi in advance for you help.</p>