calculating the number of k-combinations with and without SciPy

Question

I'm puzzled by the fact that the function comb of SciPy appears to be slower than a naive Python implementation. This is the measured time for two equivalent programs solving the Problem 53 of Project Euler.

With SciPy:

%%timeit
from scipy.misc import comb

result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if comb(n, k) > 1000000:
            result += 1
result

Output:

1 loops, best of 3: 483 ms per loop

Without SciPy:

%%timeit
from math import factorial

def comb(n, k):
    return factorial(n) / factorial(k) / factorial(n - k)

result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if comb(n, k) > 1000000:
            result += 1
result

Output:

10 loops, best of 3: 86.9 ms per loop

The second version is about 5 times faster (tested on two Macs, python-2.7.9-1, IPython 2.3.1-py27_0). Is this an expected behaviour of the comb function of SciPy (source code)? What am I doing wrong?

Edit (SciPy from the Anaconda 3.7.3-py27_0 distribution):

import scipy; print scipy.version.version
0.12.0

Edit (same difference outside IPython):

$ time python with_scipy.py
real    0m0.700s
user    0m0.610s
sys     0m0.069s

$ time python without_scipy.py
real    0m0.134s
user    0m0.099s
sys     0m0.015s

I tried this in the ipython console (command line) and got similar results for both versions (~ 90 ms per loop). — ErikR, Dec 17 '14 at 18:27
Thanks. I have added my version of SciPy, namely 0.12.0. Is this yours too? — Aristide, Dec 17 '14 at 18:32
I'm using 0.14.0 - but I would create a simple script and time it from the command line: `time python simple-script.py` — ErikR, Dec 17 '14 at 18:33

score 8 · Answer 1 · answered Dec 17 '14 at 18:48

8

It looks like you may be running the timings incorrectly and measuring the time it takes to load scipy into memory. When I run:

import timeit
from scipy.misc import comb
from math import factorial

def comb2(n, k):
    return factorial(n) / factorial(k) / factorial(n - k)

def test():
    result = 0
    for n in range(1, 101):
        for k in range(1, n + 1):
            if comb(n, k) > 1000000:
                result += 1

def test2():
    result = 0
    for n in range(1, 101):
        for k in range(1, n + 1):
            if comb2(n, k) > 1000000:
                result += 1

T = timeit.Timer(test)
print T.repeat(3,50)

T2 = timeit.Timer(test2)
print T2.repeat(3,50)

I get:

[2.2370951175689697, 2.2209839820861816, 2.2142510414123535]
[2.136591911315918, 2.138144016265869, 2.1437559127807617]

which indicates that scipy is slightly faster than the python version.

answered Dec 17 '14 at 18:48

Hooked

84,485
43
192
261

Excellent remark. Following your advice, the times are now a bit shorter, respectively 480 ms and 83 ms. The difference remains though. What is your version of SciPy? – Aristide Dec 17 '14 at 18:55
@Aristide `0.14.0`. Did you run the code I provided? What did it output for you? – Hooked Dec 17 '14 at 18:58
I have just suppressed my importations in the IPython cells, which should make the trick. Running your code actually gives: [23.93229913711548, 25.72079610824585, 24.305974006652832] and [4.27461314201355, 4.307721138000488, 4.3079118728637695]. – Aristide Dec 17 '14 at 19:03
This is quite strange, I'm not sure what to make of it, though it could be something specific to a Mac build. I tried again using `ipython` instead of python and I can repeat my results only. – Hooked Dec 17 '14 at 19:15
You say that scipy is "slightly faster", but actually `test` (around 2.22) is the scipy version and `test2` (around 2.13) is the non-scipy version, so zou should have said scipy is slightly slower. On my machine the difference is even bigger, scipy around 2.8, non-scipy around 0.8 – isarandi Sep 03 '16 at 10:39

score 4 · Accepted Answer · edited May 23 '17 at 12:00

4

Answering my own question. It seems that there exist two different functions for the same thing in SciPy. I'm not quite sure why. But the other one, binom, is 8.5 times faster than comb, and 1.5 times faster than mine, which is reassuring:

%%timeit
from scipy.special import binom
result = 0
for n in range(1, 101):
    for k in range(1, n + 1):
        if binom(n, k) > 1000000:
            result += 1
result

Output:

10 loops, best of 3: 56.3 ms per loop

SciPy 0.14.0 guys, does this work for you too?

edited May 23 '17 at 12:00

Community

1
1

answered Dec 17 '14 at 19:19

Aristide

3,606
2
30
50

Hi, I found this question interesting. My result is `100 loops, best of 3: 15.3 ms per loop` using SciPy 0.14.0, IPython 2.3.0 and Python 2.7.8 in Fedora 21. Hope this helps. – skytux Dec 17 '14 at 22:24
@skytux How does it compare to `comb` and naive implementation on your machine? – Aristide Dec 18 '14 at 07:42
1

using `comb` function I get `10 loops, best of 3: 55.3 ms per loop`. – skytux Dec 18 '14 at 20:46

score 2 · Answer 3 · answered Mar 15 '17 at 08:26

I believe this is faster than the other pure-python methods mentioned:

from math import factorial
from operator import mul
def n_choose_k(n, k):
    if k < n-k:
        return reduce(mul, xrange(n-k+1, n+1)) // factorial(k)
    else:
        return reduce(mul, xrange(k+1, n+1)) // factorial(n-k)

It's slow compared to numpy solutions, but notice that NumPy doesn't work in "unlimited integers" like Python. This means, while slow, Python will manage to return correct results. NumPy's default behavior is not always ideal for combinatorics (at least not when very large integers are involved).

E.g.

In [1]: from scipy.misc import comb

In [2]: from scipy.special import binom

In [3]: %paste
    from math import factorial
    from operator import mul
    def n_choose_k(n, k):
        if k < n-k:
            return reduce(mul, xrange(n-k+1, n+1)) // factorial(k)
        else:
            return reduce(mul, xrange(k+1, n+1)) // factorial(n-k)

## -- End pasted text --

In [4]: n_choose_k(10, 3), binom(10, 3), comb(10, 3)
Out[4]: (120, 120.0, 120.0)

In [5]: n_choose_k(1000, 250) == factorial(1000)//factorial(250)//factorial(750)
Out[5]: True

In [6]: abs(comb(1000, 250) - n_choose_k(1000, 250))
Out[6]: 3.885085558125553e+230

In [7]: abs(binom(1000, 250) - n_choose_k(1000, 250))
Out[7]: 3.885085558125553e+230

It's not surprising the error is so large, it's just the effect of truncation. Why truncate? NumPy limits itself to using 64 bits to represent an integer. The requirement for such a large integer however is quite a bit more:

In [8]: from math import log

In [9]: log((n_choose_k(1000, 250)), 2)  
Out[9]: 806.1764820287578

In [10]: type(binom(1000, 250))
Out[10]: numpy.float64

calculating the number of k-combinations with and without SciPy

With SciPy:

Without SciPy:

3 Answers3

Linked