13

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.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Daniel Walsh
  • 131
  • 3
  • 1
    Speed is generally one of the strengths of low level, strongly typed, compiled languages. – ktm Jun 20 '17 at 19:28
  • 3
    From Some Notes on Internal Implementation: "Eigenvalues and Eigenvectors use ARPACK Arnoldi methods." – Michael E2 Jun 20 '17 at 19:33
  • Related: (42736), (43851), (103559), (132176) – Michael E2 Jun 20 '17 at 19:40
  • 1
    General comment: See Band for defining diagonals in SparseArray. SparseArray[{Band[{1, 1}] -> 1., Band[{2, 1}] -> 2., Band[{1, 2}] -> 2.}, {n, n}] – Edmund Jun 20 '17 at 20:03
  • As I understand, the Arnoldi method is an iterative approach. The number of iterations required is likely to be the same in both languages. Each iteration involves a matrix multiply. I would therefore hypothesise that Mathematica is able to perform that matrix multiply faster than Python. Possibly, the multiply is done as a specialised sparse matrix multiply in Mathematica, but handed-off to ARPACK by Python (perhaps as a non-sparse matrix). You might investigate this by comparing their run times as a function of size. The sparse, banded approach should go as n, the full as n^2 – mikado Jun 20 '17 at 20:38
  • My understanding is that the problem is solved by the Arnoldi method by passing the matrix to ARPACK by Mathematica or Python, so neither language is actually handling the iteration procedure directly, rather they are passing it off to the pre-existing ARPACK routine. Why would Mathematica or Python be doing matrix multiplications instead of ARPACK? – Daniel Walsh Jun 20 '17 at 21:30
  • 1
    @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. – mikado Jun 20 '17 at 22:15
  • Also note that as soon as matrix multiplications are involved, Mathematica uses low level parallelization (multithreaded algorithm, Lapack library) to compute the products. Unless the competitor program uses similar algorithms for multiplications, it can lead to a serious advantage for MMA on most modern multicore machines. – Francois Vigneron Jun 27 '17 at 23:16
  • 3
    It appears that the scipy toarray() converts to a dense matrix which might account for the relative slowness. – Daniel Lichtblau Jul 26 '17 at 16:05
  • This also does not always hold. Sometime scipy will be faster owing to different choice of the Mathematica "BasisSize" option. It's all very problem dependent. Mathematica can also sometimes hang forever because it does something inefficient in terms of memory handling. I cannot count the times I've had to restart the entire program because of that. – b3m2a1 Aug 01 '18 at 23:56

0 Answers0