使用SciPy计算k组合数和不使用的方法

Question

使用SciPy计算k组合数和不使用的方法

6

我对于 SciPy 中的 comb 函数比起 Python 的朴素实现速度更慢感到困惑。下面是两个程序求解欧拉计划 53 题所需的时间。请注意，本文保留了 HTML 标签。

使用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

输出：

1 loops, best of 3: 483 ms per loop

没有使用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

输出：

10 loops, best of 3: 86.9 ms per loop

第二个版本快了大约5倍（在两台Mac上测试，python-2.7.9-1，IPython 2.3.1-py27_0）。这是SciPy的comb函数的预期行为吗（源代码）？我做错了什么吗？

编辑（来自Anaconda 3.7.3-py27_0分发的SciPy）：

import scipy; print scipy.version.version
0.12.0

编辑（在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

- Aristide

我在ipython控制台（命令行）中尝试了这个，两个版本的结果都很相似（每次循环约90毫秒）。 - ErikR

谢谢。我已经添加了我的SciPy版本，即0.12.0。这也是你的版本吗？ - Aristide

我正在使用0.14.0版本 - 但我想创建一个简单的脚本并从命令行计时它：time python simple-script.py。 - ErikR

3个回答

4

回答自己的问题。看起来在SciPy中存在两个不同的函数来做同样的事情。我不太确定为什么会这样。但是另一个函数binom比comb快8.5倍，比我的快1.5倍，这让人放心：

%%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

输出：

10 loops, best of 3: 56.3 ms per loop

SciPy 0.14.0 的朋友们，这对你们也适用吗？

- Aristide

嗨，我觉得这个问题很有趣。我的结果是使用SciPy 0.14.0、IPython 2.3.0和Python 2.7.8在Fedora 21上运行，每次循环15.3毫秒，取3次最佳结果。希望这可以帮到你。 - skytux

@skytux 它在您的计算机上与 comb 和朴素实现相比如何？ - Aristide

1

使用 comb 函数，我得到了 10个循环，3次中的最佳结果：每个循环55.3毫秒。 - skytux

2

我相信这比其他纯Python方法提到的方法更快：

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)

与NumPy解决方案相比，它的速度较慢，但请注意，NumPy无法像Python一样使用“无限整数”。这意味着，虽然速度慢，但Python将设法返回正确的结果。对于组合数学来说，NumPy的默认行为并不总是理想的（至少在涉及非常大的整数时）。

例如：

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

这个误差很大并不令人惊讶，这只是截断的影响。为什么要截断？NumPy限制使用64位来表示整数。然而，对于这样一个大的整数需求还远远不止这些：

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

- mathandy

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Hooked · Accepted Answer

看起来您可能在错误地运行计时器，并测量将 scipy 加载到内存中所需的时间。当我运行以下代码：

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]

这表明scipy比Python版本稍微快一些。