My question is what kind of black magic is Mathematica doing to obtain the correct answer so quickly compared to other programming languages?
Details:
I've written a Mathematica notebook to find the smallest (in magnitude) eigenvalue of a general real sparse matrix. In this experiment, I chose the matrix M to be a symmetric matrix with ones on the main diagonal, and twos just above and below the main diagonal. Here is the code:
n = 10000;
M = N[SparseArray[
Join[Table[{i, i} -> 1.0, {i, 1, n}],
Table[{i, i - 1} -> 2.0, {i, 2, n}],
Table[{i, i + 1} -> 2.0, {i, 1, n - 1}]]]];
val = Eigenvalues[M, -1][[1]]
For n=10000, the smallest eigenvalue is found almost instantly (80ms) to be val=-0.000137886. As a comparison, I tried solving the same problem in an iPython notebook using numpy and scipy.sparse.linalg. After generating the same matrix in Python, I solve it using eigs, which uses the "shift invert mode" via Arnoldi iteration:
A=coo_matrix((data, (row, col)), shape=(nn, nn)).toarray()
eigs(A, 1, sigma=0)[0][0]
which produces the correct answer:
(-0.00013788630102868301+0j)
Astonishingly, eigs takes a whopping 6 seconds!
I did a similar test in C++ using Eigen (another sparse matrix library) and the solution took a similarly unacceptable amount of time.
What I've tried
- Quitting the Mathematica Kernel to ensure that results aren't being cached.
- Introducing intentional asymmetry to eliminate the possibility that mathematica is using a fast algorithm designed for symmetric matrices (it does not affect speed).
- Asking Python for largest eigenvalue instead does not improve speed (some sources imply that the largest eigenvalue is easier than the smallest).
I would appreciate any thoughts on this matter.
Bandfor defining diagonals inSparseArray.SparseArray[{Band[{1, 1}] -> 1., Band[{2, 1}] -> 2., Band[{1, 2}] -> 2.}, {n, n}]$\endgroup$ – Edmund Jun 20 '17 at 20:03toarray()converts to a dense matrix which might account for the relative slowness. $\endgroup$ – Daniel Lichtblau Jul 26 '17 at 16:05