SuperLU is too slow

[qnzhou]

I have noticed that spsolve hangs during solve. Here are my system setting:

OS: mac 10.9.4
Python version: 2.7.6
Scipy version: 0.14.0

I have created a simple test script to illustrate the bug:
https://dl.dropboxusercontent.com/u/29899857/scipy_test.zip

With scipy 0.13, the script finishes within seconds:

  > ./scipy_test.py 
  Before solve
  After solve
  545.636841049

However, with scipy 0.14, it hangs after printing "Before solve". I have tried other matrices as well, besides the simple matrices (diagonal matrix), scipy hangs every time.

I have also observed this problem on Linux systems.

[pv]

This is due to spsolve switching from UMFPACK to SuperLU due to license reasons. You get the same issue on earlier Scipy versions with spsolve(A, b, use_umfpack=False), or if UMFPACK was not installed on the machine.

It probably does not actually hang, but the SuperLU performance for this problem is too slow to be practical. Probable culprit is column ordering algorithm, as spsolve(A, b, permc_spec="MMD_AT_PLUS_A") gives somewhat better performance than the default value.

[qnzhou]

Got it. I didn't realize umfpack support is gone since spsolve still takes the flag use-umfpack as input.

[nonhermitian]

I am not sure why SuperLU is too slow in your case. On my machine:

Numpy Version:       1.8.2
Scipy Version:       0.14.0
Cython Version:      0.20.2
Matplotlib Version:  1.3.1
Python Version:      3.4.1
Platform Info:       Darwin (x86_64)

The timings for the various re-orderings, averaged over five runs, are as follows:

COLAMD-Ordering
Solution time: 1.6238735675811768

MMD_AT_PLUS_A-Ordering
Solution time: 1.6608563899993896

MMD_ATA-Ordering
Solution time: 1.6819881916046142

NATURAL-Ordering
Solution time: 1.6404422283172608

RCM-Ordering
Solution time: 2.2440507888793944

The last case uses the RCM ordering from #3751. So it seems that the default 'COLAMD' ordering is the fastest. RCM would be the fastest, but the timing includes the amount of time it takes to find the correct permutation. These results are not so surprising as COLAMD is a fill-in minimizing reordering and the solution time scales with the amount of fill-in (NNZ) in the LU factors.

[nonhermitian]

I take my last post back, I forgot that I had scikits.umfpack installed and it is called by default if present. Indeed switching to SuperLU one finds that the solution takes much much longer to find.

Update:

The only solution I could get to finish is the MMD_AT_PLUS_A method that was ~220 times slower than the same ordering using umfpack.

[pv]

Ok, it seems then that the issue in Superlu is not only in the column ordering.
In this case, if we do not want to include non-bsd compatible things in scipy itself, scikits.umfpack should likely be polished up.

[rgommers]

@rc are you still planning to look at scikits.umfpack?

[rc]

Yes, but not before beginning of October. I have several deadlines to meet by the end of September.

[argriffing]

I'm not sure what were the plans mentioned here, but I just thought I'd bump/ping this old discussion...

[pv]

Probably the plan is to give sufficient amount of polish to https://github.com/rc/scikit-umfpack
and then make sure scipy.sparse works together with it sensibly.

Alternatively, there might be other license-compatible and more efficient sparse lu codes out there, but this doesn't seem super likely...

[rc]

FYI: The umfpack scikit is living: https://scikits.appspot.com/scikit-umfpack

[nonhermitian]

Recently, the Anaconda Python distro and Intel distro make use of the Intel MKL by default. It is possible to use ctypes to access the MKL sparse solver. In our work, we see timings many times faster than SuperLU, even when doing single-thread comparisons.

[pv]

Yes, and as noted above, UMFPACK is also faster and moreover free software and has existing bindings.

[ev-br]

Apparently, there are multiple free and non-free alternatives, all external. There's not much to do within scipy, so closing.