用Scipy计算两矩阵行积的矢量化方法

bn31dyow  于 2023-03-18  发布在  其他
关注(0)|答案(5)|浏览(151)

我想尽可能快地计算两个同维矩阵的行点积,我的方法是这样的:

import numpy as np
a = np.array([[1,2,3], [3,4,5]])
b = np.array([[1,2,3], [1,2,3]])
result = np.array([])
for row1, row2 in a, b:
    result = np.append(result, np.dot(row1, row2))
print result

当然输出是:

[ 26.  14.]
ie3xauqp

ie3xauqp1#

简单明了的方法是:

import numpy as np
a=np.array([[1,2,3],[3,4,5]])
b=np.array([[1,2,3],[1,2,3]])
np.sum(a*b, axis=1)

这避免了Python循环并且在以下情况下更快:

def npsumdot(x, y):
    return np.sum(x*y, axis=1)

def loopdot(x, y):
    result = np.empty((x.shape[0]))
    for i in range(x.shape[0]):
        result[i] = np.dot(x[i], y[i])
    return result

timeit npsumdot(np.random.rand(500000,50),np.random.rand(500000,50))
# 1 loops, best of 3: 861 ms per loop
timeit loopdot(np.random.rand(500000,50),np.random.rand(500000,50))
# 1 loops, best of 3: 1.58 s per loop
vuv7lop3

vuv7lop32#

查看numpy.einsum了解另一种方法:

In [52]: a
Out[52]: 
array([[1, 2, 3],
       [3, 4, 5]])

In [53]: b
Out[53]: 
array([[1, 2, 3],
       [1, 2, 3]])

In [54]: einsum('ij,ij->i', a, b)
Out[54]: array([14, 26])

einsum似乎比inner1d快一些:

In [94]: %timeit inner1d(a,b)
1000000 loops, best of 3: 1.8 us per loop

In [95]: %timeit einsum('ij,ij->i', a, b)
1000000 loops, best of 3: 1.6 us per loop

In [96]: a = random.randn(10, 100)

In [97]: b = random.randn(10, 100)

In [98]: %timeit inner1d(a,b)
100000 loops, best of 3: 2.89 us per loop

In [99]: %timeit einsum('ij,ij->i', a, b)
100000 loops, best of 3: 2.03 us per loop

注:NumPy是不断发展和完善的;上面所示的函数的相对性能可能在过去的几年中发生了变化。2如果性能对你来说很重要,请用你将要使用的NumPy版本运行你自己的测试。

5m1hhzi4

5m1hhzi43#

我尝试了一下,发现inner1d是最快的,但是这个函数是内部的,所以更健壮的方法是使用

numpy.einsum("ij,ij->i", a, b)

更好的方法是调整你的记忆,使求和发生在第一维,例如,

a = numpy.random.rand(3, n)
b = numpy.random.rand(3, n)

numpy.einsum("ij,ij->j", a, b)

对于10 ** 3 <= n <= 10 ** 6,这是最快的方法,比非转置的方法快两倍,最大值出现在2级缓存达到最大值时,大约为2 * 10 ** 4
还应注意,转置的sum移动比未转置的移动快得多。

该地块是用perfplot(我的一个小项目)创建的

import numpy
from numpy.core.umath_tests import inner1d
import perfplot

def setup(n):
    a = numpy.random.rand(n, 3)
    b = numpy.random.rand(n, 3)
    aT = numpy.ascontiguousarray(a.T)
    bT = numpy.ascontiguousarray(b.T)
    return (a, b), (aT, bT)

b = perfplot.bench(
    setup=setup,
    n_range=[2 ** k for k in range(1, 25)],
    kernels=[
        lambda data: numpy.sum(data[0][0] * data[0][1], axis=1),
        lambda data: numpy.einsum("ij, ij->i", data[0][0], data[0][1]),
        lambda data: numpy.sum(data[1][0] * data[1][1], axis=0),
        lambda data: numpy.einsum("ij, ij->j", data[1][0], data[1][1]),
        lambda data: inner1d(data[0][0], data[0][1]),
    ],
    labels=["sum", "einsum", "sum.T", "einsum.T", "inner1d"],
    xlabel="len(a), len(b)",
)

b.save("out1.png")
b.save("out2.png", relative_to=3)
lrl1mhuk

lrl1mhuk4#

你最好避免append,但是我想不出一个方法来避免python循环。也许是一个自定义的Ufunc?我不认为numpy.vectorize会在这里帮助你。

import numpy as np
a=np.array([[1,2,3],[3,4,5]])
b=np.array([[1,2,3],[1,2,3]])
result=np.empty((2,))
for i in range(2):
    result[i] = np.dot(a[i],b[i]))
print result

编辑

基于this answer,如果实际问题中的向量是1D,inner1d似乎可以工作。

from numpy.core.umath_tests import inner1d
inner1d(a,b)  # array([14, 26])
xxe27gdn

xxe27gdn5#

我偶然发现了这个答案,并使用运行在Python 3.5中的Numpy 1.14.3重新验证了结果。尽管我发现对于非常大的矩阵(见下面的例子),除了一个方法之外,所有方法都非常接近,以至于性能差异没有意义,但在我的系统上,上述答案大部分都成立。
对于较小的矩阵,我发现einsum是最快的,而且速度快得多,在某些情况下甚至是两倍。
我的大型矩阵示例:

import numpy as np
from numpy.core.umath_tests import inner1d
a = np.random.randn(100, 1000000)  # 800 MB each
b = np.random.randn(100, 1000000)  # pretty big.

def loop_dot(a, b):
    result = np.empty((a.shape[1],))
    for i, (row1, row2) in enumerate(zip(a, b)):
        result[i] = np.dot(row1, row2)

%timeit inner1d(a, b)
# 128 ms ± 523 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit np.einsum('ij,ij->i', a, b)
# 121 ms ± 402 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit np.sum(a*b, axis=1)
# 411 ms ± 1.99 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit loop_dot(a, b)  # note the function call took negligible time
# 123 ms ± 342 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

所以einsum在非常大的矩阵上仍然是最快的,但是快了一点,虽然这看起来是一个统计上显著的(微小的)量!

相关问题