Python中优化的点积

12

两个n维向量 u=[u1,u2,...un]v=[v1,v2,...,vn] 的点积是 u1*v1 + u2*v2 + ... + un*vn

昨天有一个问题让我寻找使用Python标准库计算点积的最快方法而不借助于第三方模块或C/Fortran/C++的调用。

我测试了四种不同的方法;目前最快的似乎是sum(starmap(mul,izip(v1,v2)))(其中starmapizip来自于itertools模块)。

对于下面所呈现的代码,这些是经过的时间(以秒为单位,对于一百万次运行):

d0: 12.01215
d1: 11.76151
d2: 12.54092
d3: 09.58523

你能想到更快的处理方法吗?

import timeit # module with timing subroutines                                                              
import random # module to generate random numnbers                                                          
from itertools import imap,starmap,izip
from operator import mul

def v(N=50,min=-10,max=10):
    """Generates a random vector (in an array) of dimension N; the                                          
    values are integers in the range [min,max]."""
    out = []
    for k in range(N):
        out.append(random.randint(min,max))
    return out

def check(v1,v2):
    if len(v1)!=len(v2):
        raise ValueError,"the lenght of both arrays must be the same"
    pass

def d0(v1,v2):
    """                                                                                                     
    d0 is Nominal approach:                                                                                 
    multiply/add in a loop                                                                                  
    """
    check(v1,v2)
    out = 0
    for k in range(len(v1)):
        out += v1[k] * v2[k]
    return out

def d1(v1,v2):
    """                                                                                                     
    d1 uses an imap (from itertools)                                                                        
    """
    check(v1,v2)
    return sum(imap(mul,v1,v2))

def d2(v1,v2):
    """                                                                                                     
    d2 uses a conventional map                                                                              
    """
    check(v1,v2)
    return sum(map(mul,v1,v2))

def d3(v1,v2):
    """                                                                                                     
    d3 uses a starmap (itertools) to apply the mul operator on an izipped (v1,v2)                           
    """
    check(v1,v2)
    return sum(starmap(mul,izip(v1,v2)))

# generate the test vectors                                                                                 
v1 = v()
v2 = v()

if __name__ == '__main__':

    # Generate two test vectors of dimension N                                                              

    t0 = timeit.Timer("d0(v1,v2)","from dot_product import d0,v1,v2")
    t1 = timeit.Timer("d1(v1,v2)","from dot_product import d1,v1,v2")
    t2 = timeit.Timer("d2(v1,v2)","from dot_product import d2,v1,v2")
    t3 = timeit.Timer("d3(v1,v2)","from dot_product import d3,v1,v2")

    print "d0 elapsed: ", t0.timeit()
    print "d1 elapsed: ", t1.timeit()
    print "d2 elapsed: ", t2.timeit()
    print "d3 elapsed: ", t3.timeit()
注意文件名必须为dot_product.py才能运行脚本;我在Mac OS X Version 10.5.8上使用Python 2.5.1。
编辑:
我对N=1000运行了这个脚本,并得到了以下结果(以秒为单位,运行一百万次):
d0: 205.35457
d1: 208.13006
d2: 230.07463
d3: 155.29670

我想可以肯定地说,事实上,选项三是最快的选项,选项二是最慢的选项(四个选项中所提供的)。


@Arrieta:你可以通过将“from dot_product”替换为“from __main__”来删除文件名必须为dot_product.py的要求。 - unutbu
@unutbu:当然,我只是认为将文件保存为那个名称以进行快速运行比更改脚本要简单。谢谢。 - Escualo
1
我的结果是: d0 耗时:13.4328830242 d1 耗时:9.52215504646 d2 耗时:10.1050257683 d3 耗时:9.16764998436 请务必检查 d1 和 d3 之间的差异是否具有统计学意义。 - liori
@liori:没错。我正在运行N=1000的问题,并期望看到更大的差异。 - Escualo
1
如果您需要重复执行点积运算,其中一个向量保持不变,那么动态编译方法可能值得探究。所有固定部分为0的项都可以完全删除,如果固定部分为1,则可以消除乘法。 - Nick Johnson
5个回答

9

仅仅为了好玩,我写了一个使用numpy的"d4":

from numpy import dot
def d4(v1, v2): 
    check(v1, v2)
    return dot(v1, v2)

以下是我在使用Python 2.5.1、XP Pro sp3和2GHz Core2 Duo T7200时得出的结果:

d0 elapsed:  12.1977242918
d1 elapsed:  13.885232341
d2 elapsed:  13.7929552499
d3 elapsed:  11.0952246724

d4已过去:56.3278584289 #使用numpy进行优化!

更有趣的是,我还打开了psyco:

d0 elapsed:  0.965477735299
d1 elapsed:  12.5354792299
d2 elapsed:  12.9748163524
d3 elapsed:  9.78255448667

d4 已经过去: 54.4599059378

基于此,我宣布d0为胜利者 :)


更新

@kaiser.se:我应该先提到我已经将所有内容转换为numpy数组:

from numpy import array
v3 = [array(vec) for vec in v1]
v4 = [array(vec) for vec in v2]

# then
t4 = timeit.Timer("d4(v3,v4)","from dot_product import d4,v3,v4")

由于其他测试中已经包括了check(v1, v2),因此我将其包含在内。如果不这样做,会使numpy获得不公平的优势(尽管它看起来可能需要一个)。数组转换减少了大约一秒钟的时间(比我预想的要少得多)。

我所有的测试都是在N=50情况下运行的。

@nikow:我正在使用numpy 1.0.4版本,这无疑是旧版本,他们在过去的一年半里肯定会改进性能。


更新#2

@kaiser.se哇,你完全正确。我一定是想着这些是列表的列表之类的东西(我真的不知道我当时在想什么...因为编程伴侣加1分)。

看起来怎么样:

v3 = array(v1)
v4 = array(v2)

新结果:

d4 elapsed:  3.22535741274

使用Psyco:

d4 elapsed:  2.09182619579
<昨天我有点烦我的numpy结果慢了,因为据说numpy用于大量计算并进行了很多优化。不过,显然我没有足够的烦恼去检查我的结果 :)>


太棒了,Seth!首先,Numpy如此缓慢令人难以置信!我本来期望它会快得多。而且,我对Psyco一无所知(尽管我认为自己是Python迷!)-感谢你指出,我一定会去看看。最后,有趣的是,基本上,Psyco为d0生成了纯优化的C代码,但不知道如何优化d3。我想这个信息是,如果你想使用Psyco,你应该布置算法,使其可以被优化(而不是将其逻辑隐藏在Python结构之后)。再次感谢你的发现! - Escualo
也许你的numpy安装有问题。在我的机器上,numpy比其他选项快得多(我没有尝试psyco)。而且N=50对于numpy来说有点太小了,无法展示它的强大之处。 - nikow
5
你做错了。应该先创建NumPy数组(而不是每次传递列表时都要进行NumPy转换),这样可以大大提高NumPy的速度。此外,去掉检查。 - u0b34a0f6ae
2
你做错了,非常错误。你正在将一个列表传递给numpy。实际上是一个由单元素numpy数组组成的列表。 - u0b34a0f6ae
谢谢更新!这只是又一个例子,说明正确使用numpy很困难。 - u0b34a0f6ae

5
这里是与numpy的比较。我们将快速starmap方法与numpy.dot进行比较。
首先,迭代普通的Python列表:
$ python -mtimeit "import numpy as np; r = range(100)" "np.dot(r,r)"
10 loops, best of 3: 316 usec per loop

$ python -mtimeit "import operator; r = range(100); from itertools import izip, starmap" "sum(starmap(operator.mul, izip(r,r)))"
10000 loops, best of 3: 81.5 usec per loop

那么,NumPy的ndarray:

$ python -mtimeit "import numpy as np; r = np.arange(100)" "np.dot(r,r)"
10 loops, best of 3: 20.2 usec per loop

$ python -mtimeit "import operator; import numpy as np; r = np.arange(100); from itertools import izip, starmap;" "sum(starmap(operator.mul, izip(r,r)))"
10 loops, best of 3: 405 usec per loop

看起来,使用numpy数组的numpy最快,其次是使用列表的python函数构造。

4

我不知道更快的方法,但我建议:

sum(i*j for i, j in zip(v1, v2))

这种写法更易读,而且不需要使用标准库模块。


@SilentGhost:你的方法需要更长时间。对于N=10,它花费了18.0258秒(一百万次运行)。我所追求的是速度;事实上可读性并不重要,因为点积是从一个函数中调用的(udotv=dot(u,v)),而且我可以在定义dot时添加尽可能多的注释。你的方法确实不合适。 - Escualo
1
@SilentGhost,一个快速的观察:将zip更改为itertools.izip可以将时间减少到15.84879。这可能值得知道。 - Escualo
2
这绝对是我会做的事情。如果在Windows上有性能问题,我会加入Psyco。 - hughdbrown
无Psyco:18.6840143091。使用Psyco:25.0433867992。这似乎是Psyco的“最坏情况”优化之一。即使使用izip()(没有Psyco),性能也只有17.4570938485。 - Seth

0
请对这个"d2a"函数进行基准测试,并将其与您的"d3"函数进行比较。
from itertools import imap, starmap
from operator import mul

def d2a(v1,v2):
    """
    d2a uses itertools.imap
    """
    check(v1,v2)
    return sum(imap(mul, v1, v2))

map(在 Python 2.x 中,我假设您使用的是这个版本)在计算之前不必要地创建了一个虚拟列表。


-2
在Mathematica中(10^12次加法和乘法):
a = RandomReal[1,{10^4,10^4}];
b = RandomReal[1,{10^4,10^4}];
AbsoluteTiming[c=a.b;]//First

9.65295 秒

(Windows 10, Dell XPS17 9700, Mathematica 12.3)


你的回答不够清晰,有一些术语需要支持性信息,请阅读:* 如何写一个好的回答? - user11717481

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接