So I wanted to test the speed of C++ vs Matlab for solving a linear system of equations. For this purpose I create a random system and measure the time required to solve it using Eigen on Visual Studio:
#include <Eigen/Core>
#include <Eigen/Dense>
#include <chrono>
using namespace Eigen;
using namespace std;
int main()
{
chrono::steady_clock sc; // create an object of `steady_clock` class
int n;
n = 5000;
MatrixXf m = MatrixXf::Random(n, n);
VectorXf b = VectorXf::Random(n);
auto start = sc.now(); // start timer
VectorXf x = m.lu().solve(b);
auto end = sc.now();
// measure time span between start & end
auto time_span = static_cast<chrono::duration<double>>(end - start);
cout << "Operation took: " << time_span.count() << " seconds !!!";
}
Solving this 5000 x 5000 system takes 6.4 seconds on average. Doing the same in Matlab takes 0.9 seconds. The matlab code is as follows:
a = rand(5000); b = rand(5000,1);
tic
x = a\b;
toc
According to this flowchart of the backslash operator:
given that a random matrix is not triangular, permuted triangular, hermitian or upper heisenberg, the backslash operator in Matlab uses a LU solver, which I believe is the same solver that I'm using on the C++ code, that is, lu().solve
Probably there is something that I'm missing, because I thought C++ was faster.
- I am running it with release mode active on the Configuration Manager
- Project Properties - C/C++ - Optimization - /O2 is active
- Tried using Enhanced Instructions (SSE and SSE2). SSE actually made it slower and SSE2 barely made any difference.
- I am using Community version of Visual Studio, if that makes any difference
m.lu().solve(b);
tom.llt().solve(b);
and now it takes 0.1sec. – Outhe-O3 -march=native
to let the C++ compiler make use of your architecture. – Pentosanm.colPivHouseholderQr().solve(b)
, see their documentation. They have various algorithms implemented. Perhaps iterative algorithms are interesting as well, like Conjugate Gradient, BiCGStab or even some fancy multigrid methods? I am not sure what Eigen implements, though. – Pentosan\
operator in Matlab, also known asmldivide
: mathworks.com/help/matlab/ref/mldivide.html If I am not wrong, it seems like Matlab would be doing LU decomposition in your case. – Errandb = rand(5000, 1)
shouldn't this beb = rand(5000, 5000)
to match the C++ code? – Respectfultic
&toc
is not the propper way to measure performance of matlab code. Have a look attimeit()
: de.mathworks.com/help/matlab/ref/timeit.html – Dovetaila
is a5000x5000
matrix andb
is the5000x1
right side vector for the equationa*x=b
. – Dovetail