{"id":48,"date":"2012-12-19T20:38:00","date_gmt":"2012-12-19T20:38:00","guid":{"rendered":"https:\/\/sandrock.co.za\/carl\/2012\/12\/19\/numpy-performance-games\/"},"modified":"2012-12-19T20:38:00","modified_gmt":"2012-12-19T20:38:00","slug":"numpy-performance-games","status":"publish","type":"post","link":"https:\/\/sandrock.co.za\/carl\/2012\/12\/numpy-performance-games\/","title":{"rendered":"Numpy performance games"},"content":{"rendered":"

I was reading a blog post about Numba vs Cython<\/a>, which reported a huge speedup in using Numba<\/a> with a simply Python code for calculating the pairwise distance between a number of points. The code they posted was as follows:<\/p>\n

<\/p>\n

import numpy as np\n\ndef pairwise_python(X, D):\n    M = X.shape[0]\n    N = X.shape[1]\n    for i in xrange(M):\n        for j in xrange(M):\n            d = 0.0\n            for k in xrange(N):\n                tmp = X[i, k] - X[j, k]\n                d += tmp * tmp\n            D[i, j] = np.sqrt(d)\n<\/pre>\n
\nThey tested it with the following code:<\/p>\n
N = 1000\nX = np.random.random((N, 3))\nD = np.empty((N, N))\n<\/pre>\n
%timeit pairwise_python(X, D)\n\n 1 loops, best of 3: 8.3 s per loop\n<\/pre>\n
I thought I could probably do a bit better by making more use of the built-in numpy functions, so I rewrote it like thise:
\n<\/p>\n
def pairwise_numpy(X, D):\n    M, N = X.shape\n    for i in xrange(M):\n        for j in xrange(M):\n            D[i, j] = np.linalg.norm(X[i, :] - X[j, :])<\/pre>\n
%timeit pairwise_numpy(X, D)\n\n1 loops, best of 3: 13.8 s per loop\n<\/pre>\n
Imagine my surprise. My “optimised” version takes almost twice as long! Shocked, I thought I shold try to go all the way
\n<\/p>\n
def pairwise_numpy_fast(X, D):\n    for i, row in enumerate(X):\n        D[i, :] = np.sqrt(np.sum((X - row)**2, 1))\n<\/pre>\n
%timeit pairwise_numpy_fast(X, D)\n\n10 loops, best of 3: 46.9 ms per loop\n<\/pre>\n
That’s a much more satisfying result. Goes to show how much one can gain without even going too crazy with tricks.<\/p>\n
Then there is my current approach, using Fortran for the tight loops. Notice how much the Fortran source resembles the Python ones. This is not your daddy’s Fortran, it’s Fortran 90.
\n<\/p>\n
subroutine pairwise(X, M, N, D)\n  double precision, intent (IN) :: X(M, N)\n  integer, intent(IN) :: M, N\n  double precision, intent (OUT) :: D(M, M)\n  integer :: i, j\n\n  forall (i = 1:M, j = 1:M)\n     D(i, j) = sqrt(sum((X(i, :) - X(j, :))**2))\n  end forall\nend subroutine pairwise\n<\/pre>\n
After compiling this file to a Python module with f2py (<\/p>\n
f2py -c -m pairwise_fortran pairwise.f90<\/pre>\n
), we can check it with
\n<\/p>\n
from pairwise_fortran import pairwise as pairwise_fortran\n%timeit pairwise_fortran(X, M, N)\n\n100 loops, best of 3: 6.58 ms per loop\n<\/pre>\n
Now that’s more like it!<\/p>\n
So, from 8.3 s down to 6.6 ms. That’s three orders of magnitude from the original naive Python code, but the “clever” numpy version is less than one order slower. My feeling is that the benefits gained from a language like Python more than make up for any speed deficits. Of course, there’s the distinct advantage of having a really quick way to incorporate the efficient compilers available for Fortran using a tool like f2py.<\/p>\n","protected":false},"excerpt":{"rendered":"
I was reading a blog post about Numba vs Cython, which reported a huge speedup in using Numba with a simply Python code for calculating the pairwise distance between a number of points. The code they posted was as follows: import numpy as np def pairwise_python(X, D): M = X.shape[0] N = X.shape[1] for i in […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_feature_clip_id":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_post_was_ever_published":false},"categories":[1],"tags":[],"class_list":["post-48","post","type-post","status-publish","format-standard","hentry","category-uncategorized"],"jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/posts\/48","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/comments?post=48"}],"version-history":[{"count":0,"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/posts\/48\/revisions"}],"wp:attachment":[{"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/media?parent=48"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/categories?post=48"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/sandrock.co.za\/carl\/wp-json\/wp\/v2\/tags?post=48"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}