12
$\begingroup$

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

  1. Quitting the Mathematica Kernel to ensure that results aren't being cached.
  2. Introducing intentional asymmetry to eliminate the possibility that mathematica is using a fast algorithm designed for symmetric matrices (it does not affect speed).
  3. 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.

Improve this question
$\endgroup$
10
  • 1
    $\begingroup$ Speed is generally one of the strengths of low level, strongly typed, compiled languages. $\endgroup$ – ktm Jun 20 '17 at 19:28
  • 3
    $\begingroup$ From Some Notes on Internal Implementation: "Eigenvalues and Eigenvectors use ARPACK Arnoldi methods." $\endgroup$ – Michael E2 Jun 20 '17 at 19:33
  • 1
    $\begingroup$ General comment: See Band for defining diagonals in SparseArray. SparseArray[{Band[{1, 1}] -> 1., Band[{2, 1}] -> 2., Band[{1, 2}] -> 2.}, {n, n}] $\endgroup$ – Edmund Jun 20 '17 at 20:03
  • 1
    $\begingroup$ @DanielWalsh my understanding is that an Arnoldi method does not need to be given the matrix M explicitly. It only needs a function x->M.x. I believe this is referred to as "Reverse Communication" in the ARPACK user guide. Presumably, Mathematica is making more intelligent use of the library than Python. $\endgroup$ – mikado Jun 20 '17 at 22:15
  • 3
    $\begingroup$ It appears that the scipy toarray() converts to a dense matrix which might account for the relative slowness. $\endgroup$ – Daniel Lichtblau Jul 26 '17 at 16:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.