Update/expand the PerformanceTips wiki page

[ev-br]

Performance Tips page in the scipy wiki(
http://wiki.scipy.org/PerformanceTips#head-a5a8d210cc636efc34184be66171fd55b9aeb35e) is very useful but might benefit from a bit a of face-lift. The wiki does not seem to be easily editable anymore, so maybe the content should be moved to scipy central and/or linalg tutorials.

[ilayn]

@ev-br Since wiki is now dumped should we close this?

[mdhaber]

That link no longer works, but it's available at
https://scipy.github.io/old-wiki/pages/PerformanceTips

I suppose the bit about linear algebra on large arrays might be useful in the linalg tutorials. Do you still think so @ev-br?

[ev-br]

Agreed, this content (or a part which is still relevant, I did not check for quite a while), could be moved to the linalg tutorials.

[lucascolley]

Here is the linalg section:

Linear Algebra on Large Arrays
Sometimes the performance or memory behaviour of linear algebra on large arrays can suffer because of copying behind the scenes. Consider for example, computing the dot product of two large matrices with a small result:

>>> N = 1e6
>>> n = 40
>>> A = np.ones((N,n))
>>> C = np.dot(A.T, A)

Although C is only 40 by 40, inspecting the memory usage during the operation of dot will indicate that a copy is being made. The reason is that the dot product uses underlying BLAS operations which depend on the matrices being stored in contiguous C order. The transpose of operation does not effect a copy, but places the transposed array in Fortran order:

>>> A.flags
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
   ...
>>> A.T.flags
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True
   ...

If we give dot two matrices that are both C_CONTIGUOUS, then the performance is better:

>>> from time import time
>>> AT_F = np.ones((n,N), order='F')
>>> AT_C = np.ones((n,N), order='C')
>>> t = time();C = np.dot(A.T, A);time() - t
3.9203271865844727
>>> t = time();C = np.dot(AT_F, A);time() - t
3.9461679458618164
>>> t = time();C = np.dot(AT_C, A);time() - t
2.4167969226837158

This is not ideal, however, because it required us to start with two arrays. Is there not some way of computing dot(A.T, A) without making an extra copy? There is support for this in the BLAS through the _SYRK method (see http://www.netlib.org/blas/blasqr.pdf). SciPy exposes this through scipy.linalg.blas (which can be used with C_CONTIGUOUS arrays and F_CONTIGUOUS arrays respectively). Alas, this function is not yet implemented, however, _GEMM is it is now (scipy 0.13.0). Here we use it:

>>> import scipy.linalg.blas
>>> t = time();C = scipy.linalg.blas.dgemm(alpha=1.0, a=A.T, b=A.T, trans_b=True);time() - t
2.412520170211792

This gives the same performance as dot but with the advantage that we did not need to make an extra copy. Note that we passed A.T here so that the array was in Fortran order because blas is the Fortran BLAS. One must be careful about this: passing the wrong type of array will not realize the performance gains. Here, for example, we pass the C_CONTIGUOUS arrays for which copies must me made.

>>> t = time();C = scipy.linalg.blas.dgemm(alpha=1.0, a=A, b=A, trans_a=True);time() - t
5.0581560134887695

There is not much documentation about this, but some useful discussions can be found in the mailing list archives. Another related question is about inplace operations (again, avoiding copying). Some functions provide an option for this, and again, some of the BLAS functions can be used for this. Here are some related discussions:

Enhancing dot()

numpy large arrays?

dgemm overwrite_c bug?

inplace dot products

making numpy.dot faster

Inplace dot and blas access