numpy.linalg.cholesky performance

[pv]

Both scipy.linalg.cholesky and numpy.linalg.cholesky call DPOTRF. However, the timings are

In [1]: A = np.random.rand(500,500); A = A.dot(A.T)
In [2]: %timeit scipy.linalg.cholesky(A)
100 loops, best of 3: 2.2 ms per loop
In [3]: %timeit numpy.linalg.cholesky(A)
100 loops, best of 3: 4.19 ms per loop

Perhaps things can be improved in Numpy?

[bsipocz]

timings seem to have improved over the years, I don't see the 2x factor any more for the example above:

In [4]: A = np.random.rand(500,500); A = A.dot(A.T)                                                                     

In [5]: %timeit scipy.linalg.cholesky(A)                                                                                
2.26 ms ± 180 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [6]: %timeit numpy.linalg.cholesky(A)                                                                                
2.08 ms ± 43.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [7]: A = np.random.rand(1500,1500); A = A.dot(A.T)                                                                   

In [8]: %timeit numpy.linalg.cholesky(A)                                                                                
34.6 ms ± 607 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [9]: %timeit scipy.linalg.cholesky(A)                                                                                
27.3 ms ± 904 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

[melissawm]

The issue is still there for me, but the difference between NumPy and SciPy seems to get smaller as the matrix size increases. On (Manjaro) Linux, with NumPy 1.19.1/OpenBLAS from conda-forge:

In [1]: A = np.random.rand(500,500); A = A.dot(A.T)

In [2]: %timeit np.linalg.cholesky(A)
6.53 ms ± 78.7 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [3]: %timeit sp.linalg.cholesky(A)
3.52 ms ± 36.2 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [4]: A = np.random.rand(2000,2000); A = A.dot(A.T)

In [5]: %timeit np.linalg.cholesky(A)
126 ms ± 372 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [6]: %timeit sp.linalg.cholesky(A)
84 ms ± 1.27 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [7]: A = np.random.rand(10000,10000); A = A.dot(A.T)

In [8]: %timeit np.linalg.cholesky(A)
8.01 s ± 111 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [9]: %timeit sp.linalg.cholesky(A)
7.44 s ± 213 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[bsipocz]

Ahh, interesting. I'm on OSX, using the wheels from PyPI.

[rgommers]

The difference is smaller for me too. This is simply a consequence of how we do wrapping of BLAS/LAPACK libraries I think:

>>> %timeit scipy.linalg.cholesky(A)
3.52 ms ± 52.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
>>> %timeit numpy.linalg.cholesky(A)
4.12 ms ± 49 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

The difference is smaller than that between using different BLAS/LAPACK libraries. So let's close this, I don't think there is anything still actionable here.