11
$\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.

$\endgroup$
  • 1
    $\begingroup$ Speed is generally one of the strengths of low level, strongly typed, compiled languages. $\endgroup$ – user6014 Jun 20 '17 at 19:28
  • 2
    $\begingroup$ From Some Notes on Internal Implementation: "Eigenvalues and Eigenvectors use ARPACK Arnoldi methods." $\endgroup$ – Michael E2 Jun 20 '17 at 19:33
  • $\begingroup$ Related: (42736), (43851), (103559), (132176) $\endgroup$ – Michael E2 Jun 20 '17 at 19:40
  • 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
  • 2
    $\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.