# How can I improve the speed of eigenvalue decompositions for large matrices?

I often need to compute the eigenvalues of large matrices, and I invariably resort to MATLAB for these, simply because it is much faster. I'd like to change that, so that I can work entirely inside my notebook.

Here's a plot comparing the timings between the two for eigenvalue decompositions for matrices of varying sizes (left). The y-axis shows the time in seconds. As you can see, there's about a factor 3 difference between the two (right).

Here's a sample code in Mathematica to generate timings:

timings = With[{x = RandomReal[NormalDistribution[], {#, #}]},
Eigenvalues[x]; // Timing // First] & /@ Range[500,5000,500]


and its equivalent in MATLAB:

s = 500:500:5000;
t = zeros(numel(s),1);
for i = 1:numel(s)
x=randn(s(i));
t1=tic;eig(x);t(i)=toc(t1);
end


I do not think that Mathematica's algorithms are inefficient, as the fastest algorithms for eigenvalue decompositions (in the general case, not exploiting symmetry and such) are $\mathcal{O}(N^{2.376})$ and the timings both MATLAB's and Mathematica's implementations have the same correct slope on a log-log plot.

I suspected unpacking in the background during the call to Eigenvalues and turning On["Packing"] confirms this. However, I don't think this alone could be the cause for a 3 fold speed reduction. I'm not expecting the timings to be exact either, as I understand that arrays and matrices are baked into the core of one and not the other, which can lead to performance differences.

However, I'm interested in knowing if

1. there are reasons other than the simplified one I gave above for the difference in timings and
2. there ways in which I can improve the speeds or at least, reduce the difference by some amount. Or is this something that one has to accept as a fact of life?
• It seems like Matlab and Mathematica use the same core algorithms, but the constants involved are different. Jan 18, 2012 at 2:52
• I get different numbers. M.8.0 vs Matlab R2009b - the results are roughly similar, with Mathematica even a tiny bit faster. I tested for matrices of size 1000, 2000 and 4000, but separately, not in a loop. As for unpacking, it indeed happens, but on the level of the final list of eigenvalues only, since eigenvalues are generally complex. This does not contribute much to running time. So, either the newer versions of Matlab got dramatically faster, or I don't know what could cause such a difference. I am on Win7, 64 bits, AMD Phenom II X6 2.8 GHz, 8Gb RAM. Jan 18, 2012 at 4:01
• Which version of Mathematica and Matlab have been used for this test?
– user21
Jan 18, 2012 at 8:10
• @ruebenko Mathematica v8 and MATLAB 2011b. However, I must say that they were not "tested" at the same time. I've been collecting stuff and data for questions for this site over the course of time and the state of my machine (other processes) was different for each. So it is possible that the timings are not as drastic, but I do know that MATLAB certainly runs faster on mine. I will update when I get a chance to start with a clean slate. However, I accepted Mark's answer, because that's sort of what I was getting at. Perhaps I'll remove the "how to improve" part to be in tune with his answer
– rm -rf
Jan 18, 2012 at 10:00

Mathematica is every bit as fast as Matlab for these types of computations. The source of the discrepancy arises from the fact that Timing keeps track of total time used by all processors when Mathematica distributes the computation across them. We can examine a fair comparison using AbsoluteTiming, which is more comparable to Matlab's tic and toc.

Consider the following computed on my Macbook Pro:

t1 = First[Timing[Eigenvalues[RandomReal[{0, 1},
{1000, 1000}]]]];
t2 = First[AbsoluteTiming[Eigenvalues[RandomReal[{0, 1},
{1000, 1000}]]]];
{t1, t2}


{5.16576, 1.329784}

Again, the only difference is the use of Timing versus AbsoluteTiming. You can watch the wall clock to convince yourself that the faster time is accurate. Let's try this with with the OP's code:

timingsGood = With[{x = RandomReal[NormalDistribution[], {#, #}]},
Eigenvalues[x]; // AbsoluteTiming // First] & /@
Range[500, 5000, 500];
timingsBad = With[{x = RandomReal[NormalDistribution[], {#, #}]},
Eigenvalues[x]; // Timing // First] & /@
Range[500, 5000, 500];

Note that the (incorrect) Timing result is always consistently about three times longer than the (correct) AbsoluteTiming result, which accounts just about exactly for the OP's observations.
Somewhat tangentially, if you open the help center and paste tutorial/LinearAlgebraInMathematicaOverview in search (link), you'll find a rather extensive discussion of linear algebra in mma. While I have not used mma seriously for dense matrices, I know that there are options for the sparse solvers described in this tutorial which are not described elsewhere (eg the options to the Arnoldi method--at least, I do not remember finding them elsewhere).