Python中用于计算余弦距离的优化方法

7
我写了一个计算两个数组之间余弦距离的方法:
def cosine_distance(a, b):
    if len(a) != len(b):
        return False
    numerator = 0
    denoma = 0
    denomb = 0
    for i in range(len(a)):
        numerator += a[i]*b[i]
        denoma += abs(a[i])**2
        denomb += abs(b[i])**2
    result = 1 - numerator / (sqrt(denoma)*sqrt(denomb))
    return result

在大型数组上运行它可能非常缓慢。是否有一种优化版本的方法可以更快地运行?

更新:我已经尝试了到目前为止所有的建议,包括scipy。这是需要超越的版本,结合了Mike和Steve的建议:

def cosine_distance(a, b):
    if len(a) != len(b):
        raise ValueError, "a and b must be same length" #Steve
    numerator = 0
    denoma = 0
    denomb = 0
    for i in range(len(a)):       #Mike's optimizations:
        ai = a[i]             #only calculate once
        bi = b[i]
        numerator += ai*bi    #faster than exponent (barely)
        denoma += ai*ai       #strip abs() since it's squaring
        denomb += bi*bi
    result = 1 - numerator / (sqrt(denoma)*sqrt(denomb))
    return result

a和b是复数数组吗? - John La Rooy
我不明白为什么你要先取绝对值再平方。 - John La Rooy
1
我刚刚进行了一个快速测试,当列表中有大约1000个元素时,使用numpy更快。 - John La Rooy
2
NumPy 在处理小数组时较慢的原因是由于转换为 NumPy 数组的开销。 - John La Rooy
只要你试图刮掉代码中的东西,如果你仍在使用Python 2.x,则可以尝试使用xrange()而不是range()。如果你正在使用Python 3,则只有range()并且它返回一个迭代器。 - steveha
显示剩余9条评论
8个回答

7
如果你可以使用SciPy,你可以从`spatial.distance`中使用`cosine`函数:
请参考http://docs.scipy.org/doc/scipy/reference/spatial.distance.html
如果你不能使用SciPy,你可以尝试通过重写Python代码来获得一些小的加速(编辑:但它并没有像我想象的那样有效果,详见下文)。
from itertools import izip
from math import sqrt

def cosine_distance(a, b):
    if len(a) != len(b):
        raise ValueError, "a and b must be same length"
    numerator = sum(tup[0] * tup[1] for tup in izip(a,b))
    denoma = sum(avalue ** 2 for avalue in a)
    denomb = sum(bvalue ** 2 for bvalue in b)
    result = 1 - numerator / (sqrt(denoma)*sqrt(denomb))
    return result

当a和b的长度不匹配时,最好引发异常。

通过在调用sum()时使用生成器表达式,您可以使用Python内部的C代码完成大部分工作来计算您的值。这应该比使用for循环更快。

我没有计时,所以无法猜测它可能会快多少。但是,SciPy代码几乎肯定是用C或C++编写的,因此它应该是您可以得到的最快速度。

如果您在Python中进行生物信息学,则确实应该使用SciPy。

编辑:Darius Bacon测试了我的代码并发现它更慢。所以我测试了我的代码,...是的,它更慢。对于所有人的教训:当您尝试加快速度时,请勿猜测,而要进行测量。

我对于为什么我的尝试将更多工作放在Python的C内部上变慢感到困惑。我尝试了长度为1000的列表,它还是更慢。

我不能再花更多时间来试图巧妙地黑客Python了。如果您需要更快的速度,我建议您尝试SciPy。

编辑:我刚刚手动测试了一下,没有使用timeit。我发现对于短的a和b,旧代码更快;对于长的a和b,新代码更快;在两种情况下,差异不大。(我现在在想我是否可以相信我的Windows计算机上的timeit;我想在Linux上再次尝试此测试。)我不会更改工作代码以尝试使其更快。再次敦促您尝试SciPy。 :-)


分子线是不正确的:它执行的是嵌套循环而不是并行循环。 - Darius Bacon
1
另外,当我修复了那行代码以获得正确的答案时,它仍然比原始代码慢。无论如何,关于SciPy我们达成一致!(分子= sum(avalue * bvalue for avalue, bvalue in zip(a, b))) - Darius Bacon
使用SciPy很好。不幸的是,您的非SciPy重写返回了错误的值。将分子行替换为gnibbler的结果可以得到正确的答案,但实际上比我的原始代码慢得多。 - Dan
1
有趣的是,scipy实际上要慢得多。为了测试,我将一些小的数组运行100,000次迭代。原始代码运行约1.3秒,而scipy则在约7.5秒运行完成。不知道在更大的数组上会不会有所改观? - Dan
scipy.spatial.distance非常慢。这并不奇怪:它只是一个普通的Python函数 :(。 - nh2
显示剩余2条评论

6

起初我认为,如果你不像numpy或scipy那样跳出Python语言界限,或者不改变计算方式,是不可能大幅提高运算速度的。但是,以下是我尝试解决这个问题的方案:

from itertools import imap
from math import sqrt
from operator import mul

def cosine_distance(a, b):
    assert len(a) == len(b)
    return 1 - (sum(imap(mul, a, b))
                / sqrt(sum(imap(mul, a, a))
                       * sum(imap(mul, b, b))))

使用500k元素数组,Python 2.6的速度大约是原来的两倍。(在Jarret Hardie的建议下将map更改为imap后)

以下是原帖修订版代码的微调版本:

from itertools import izip

def cosine_distance(a, b):
    assert len(a) == len(b)
    ab_sum, a_sum, b_sum = 0, 0, 0
    for ai, bi in izip(a, b):
        ab_sum += ai * bi
        a_sum += ai * ai
        b_sum += bi * bi
    return 1 - ab_sum / sqrt(a_sum * b_sum)

虽然它看起来不太好看,但确实更快……

编辑:还可以尝试使用Psyco!它可以将最终版本的速度提高4倍。我怎么会忘记呢?


很不错的补充 - 很高兴听到使用imap对mul over ** 2有优势。 - Jarret Hardie
我不认为它那么丑 :p - John La Rooy
我有点失望地看到命令式代码击败了更直接表达问题的纯函数式代码。 - Darius Bacon

2

如果你要对a[i]b[i]进行平方,则不需要取绝对值。

a[i]b[i]存储在临时变量中,以避免多次索引。也许编译器可以优化此操作,但也可能不行。

检查**2运算符。它是否简化为乘法,还是使用了一般的幂函数(对数-乘以2-反对数)。

不要两次执行平方根(尽管其成本很小)。请执行sqrt(denoma * denomb)


好的调用...这些都稍微节省了一点时间。 - Dan
@Dan:欢迎。接下来我会看看是否可以展开循环,以防迭代器成为瓶颈(它们往往会这样)。然后我会进行一些堆栈采样,以查看函数是否被调用得比必要的次数多(或者是否存在其他未被注意到的时间瘤)。 - Mike Dunlavey

1

对于大约1000个以上的数组,这种方法更快。

from numpy import array
def cosine_distance(a, b):
    a=array(a)
    b=array(b)
    numerator=(a*b).sum()
    denoma=(a*a).sum()
    denomb=(b*b).sum()
    result = 1 - numerator / sqrt(denoma*denomb)
    return result

1
类似于Darius Bacon的回答,我一直在使用operator和itertools来产生更快的答案。根据timeit的测试结果,在一个500项的数组中,以下方法似乎比原来的方法快1/3:
from math import sqrt
from itertools import imap
from operator import mul

def op_cosine(a, b):
    dot_prod = sum(imap(mul, a, b))
    a_veclen = sqrt(sum(i ** 2 for i in a))
    b_veclen = sqrt(sum(i ** 2 for i in b))

    return 1 - dot_prod / (a_veclen * b_veclen)

1

在处理长输入数组时,使用SciPy内部的C代码效果最佳。而在处理短输入数组时,使用简单直接的Python代码效果最佳;Darius Bacon基于izip()的代码表现最佳。因此,最终的解决方案是根据输入数组的长度,在运行时决定使用哪种方法:

from scipy.spatial.distance import cosine as scipy_cos_dist

from itertools import izip
from math import sqrt

def cosine_distance(a, b):
    len_a = len(a)
    assert len_a == len(b)
    if len_a > 200:  # 200 is a magic value found by benchmark
        return scipy_cos_dist(a, b)
    # function below is basically just Darius Bacon's code
    ab_sum = a_sum = b_sum = 0
    for ai, bi in izip(a, b):
        ab_sum += ai * bi
        a_sum += ai * ai
        b_sum += bi * bi
    return 1 - ab_sum / sqrt(a_sum * b_sum)

我制作了一个测试套件,用不同长度的输入对函数进行测试,并发现在长度约为200左右时,SciPy函数开始获胜。输入数组越大,它的优势就越大。对于非常短的长度数组,比如3,简单的代码会赢得胜利。该函数添加了一点点开销来决定最佳方法,然后执行它。
如果您感兴趣,这里是测试套件:
from darius2 import cosine_distance as fn_darius2
fn_darius2.__name__ = "fn_darius2"

from ult import cosine_distance as fn_ult
fn_ult.__name__ = "fn_ult"

from scipy.spatial.distance import cosine as fn_scipy
fn_scipy.__name__ = "fn_scipy"

import random
import time

lst_fn = [fn_darius2, fn_scipy, fn_ult]

def run_test(fn, lst0, lst1, test_len):
    start = time.time()
    for _ in xrange(test_len):
        fn(lst0, lst1)
    end = time.time()
    return end - start

for data_len in range(50, 500, 10):
    a = [random.random() for _ in xrange(data_len)]
    b = [random.random() for _ in xrange(data_len)]
    print "len(a) ==", len(a)
    test_len = 10**3
    for fn in lst_fn:
        n = fn.__name__
        r = fn(a, b)
        t = run_test(fn, a, b, test_len)
        print "%s:\t%f seconds, result %f" % (n, t, r)

0
def cd(a,b):
    if(len(a)!=len(b)):
        raise ValueError, "a and b must be the same length"
    rn = range(len(a))
    adb = sum([a[k]*b[k] for k in rn])
    nma = sqrt(sum([a[k]*a[k] for k in rn]))
    nmb = sqrt(sum([b[k]*b[k] for k in rn]))

    result = 1 - adb / (nma*nmb)
    return result

你正在调用sum()函数时使用列表推导式。这将创建一个列表,然后sum()函数将使用该列表一次,然后该列表将被垃圾回收。Python有一个很棒的功能叫做“生成器表达式”,您可以使用与列表推导式相同的语法,但它将创建一个迭代器。如果您只是从调用sum()函数的内部删除[],那么您现在将使用生成器表达式。在此处阅读更多信息:http://docs.python.org/howto/functional.html#generator-expressions-and-list-comprehensions - steveha
@steveha:这取决于输入长度和函数。我不知道这里的情况,但对于短输入(len ~100),str.join(..)与列表理解比genexps更快。 - u0b34a0f6ae
@kaizer.se:str.join是一个特殊情况,因为当它有一个列表时,它首先求出长度的总和,然后创建一个总大小的字符串并用这些部分填充它;否则,它按照显而易见的方式构建字符串(对于可迭代的每个部分:result+= part)。 - tzot

0

您更新的解决方案仍然有两个平方根。您可以通过将sqrt行替换为以下内容将其减少到一个:

result = 1 - numerator / (sqrt(denoma*denomb))

乘法通常比sqrt快得多。虽然在函数中只调用一次,但它听起来像是在计算大量余弦距离,因此改进会累积。

您的代码看起来应该适合矢量优化。因此,如果跨平台支持不是问题,并且您想进一步加速它,则可以使用C编写余弦距离代码,并确保编译器积极地对生成的代码进行矢量化(即使Pentium II也能够进行一些浮点矢量化)


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