Hey guys,

I have been working on a project with Matlab, and the main time-consuming step is matrix inversion. I thought I should be able to speed this up by migrating to c++, and I have just rewritten the code into c++ with Eigen library. To my big surprise, with Eigen it takes twice or more time to compute the inverse. As a test, one can try the following c++ file and Matlab file:
---------------------------
//test.cpp
#include <Eigen/Core>
#include <Eigen/LU>
#include <iostream>
#include <ctime>
using namespace std;
USING_PART_OF_NAMESPACE_EIGEN
#define N 500
int main(int, char *[])
{
Matrix<double,N,N> m;
for (int i = 0; i < N; i++){
for (int j = 0; j < N; j++)
m(i,j) = ( (double)rand() / (double)RAND_MAX ) - 0.5;
}
Matrix<double,N,N> mInv;
time_t t1 = clock();
mInv = m.inverse();
time_t t2 = clock();
cout << "Time elapsed: " << (double)(t2 - t1) / (double) CLOCKS_PER_SEC << endl;
}
-------------------------------------
%test.m
clear all;
N = 500;
for i = 1 : N
for j = 1 : N
A(i,j) = rand(1) - 0.5;
end
end
tic
B = inv(A);
toc
------------------------------------------
By the way my laptop has core 2 Duo CPU T9300 @ 2.50GHz. I have set SSE2 and tried -O2 -O3 with g++ or release mode with VC++ 2008. I have tried the devel branch and tried both inverse() and computeInverse(). In addition, I am not solving linear equations, so I really need matrix inversion. Am I missing something here? Thanks a lot!

Rick

First of all, matrices of such a large size should be dynamic-size, not fixed-size! Matlab definitely uses dynamic-size. So:

Code: Select all

MatrixXd m(N,N), mInv(N,N);


Other things to try: check gcc version >= 4.2, pass -DNDEBUG, pass just -O2 not -O3, make sure you've passed -msse2 (just checking).

Finally:

In addition, I am not solving linear equations, so I really need matrix inversion. Am I missing something here?

You definitely don't need matrix inversion to solve systems: for that, you should use solve() directly. Like this: to solve AX=B, do:

Code: Select all

X = A.solve(B);


I get similar results. I changed the test program to compute the inverse 100 times, to make sure I get meaningful timings. I also changed to dynamic matrices in Eigen.

The results are 3.3 sec with Matlab (R2009b) and 8.7 sec with Eigen (development branch, GCC 4.3.3). It may be possible that a newer version of GCC yields faster programs. Eigen's stable branch is vastly slower (28 sec), because it uses complete pivoting (I think).

A partial explanation is that fairly modern versions of Matlab use Intel Math Kernel Library (MKL) with support for parallellization. My computer is also Core 2 Duo, so that explains a factor of 2 at most. It seems there is still room for optimization.

If matlab is using a multicore implementation, that should be noticable in the task manager.

Incidentally, what do you need the inverse for? Apparently, for most uses, there are better alternatives.

Thanks everyone! I tried dynamic size but it does not help much. And thanks jitseniesen for testing it out. Just found out that Matlab inv() actually calls LAPACK routines, which probably also contributes to its speed. (http://www.mathworks.com/access/helpdes ... f/inv.html)

EamonNerbonne wrote:If matlab is using a multicore implementation, that should be noticable in the task manager.

Depends. Most task managers (don't know about windows) only show processes, not threads. Matlab probably uses only one process, but that uses N threads.

Yes, that's likely the explanation. We know that eventually we have to parallelize algorithms in Eigen, it's just not among our 3.0 priorities. That said, the algorithm is already blocked, so more than 50% of the work is already done.

well, the decomposition is fast, the solver is fast too, but the inverse function is not fully optimized. Currently it is implemented as lu.solve(identity) while a faster approach would be to implement a fast inverse function for triangular matrices, and then do something like L.solve(U.inverse()). This is not hard to do but definitely not a priority as in most cases computing the inverse matrix is not needed.

Thanks, but I would really appreciate it if some one could optimize this. In my situation, I have to compute a physics formula which is a series of matrix inversion operation...

ggael wrote:well, the decomposition is fast, the solver is fast too, but the inverse function is not fully optimized. Currently it is implemented as lu.solve(identity) while a faster approach would be to implement a fast inverse function for triangular matrices, and then do something like L.solve(U.inverse()). This is not hard to do but definitely not a priority as in most cases computing the inverse matrix is not needed.

bjacob wrote:

EamonNerbonne wrote:If matlab is using a multicore implementation, that should be noticable in the task manager.

Depends. Most task managers (don't know about windows) only show processes, not threads. Matlab probably uses only one process, but that uses N threads.

I was referring to the fact that you'll see a spike in usage for several seconds (since that's how long this particular inverse takes to compute, apparently) across either 1 core or several cores - assuming they're otherwise idle.