python中巨大矩阵的点积的行和

Posted 2023-03-13

【中文标题】python中巨大矩阵的点积的行和

【英文标题】：Row Sum of a dot product for huge matrix in python

【发布时间】：2017-07-30 04:19:15

【问题描述】：

我有 2 个矩阵 100kx200 和 200x100k 如果它们是小矩阵，我会使用 numpy 点积

sum(a.dot(b), axis = 0)

但是矩阵太大了，而且我不能使用循环，有没有聪明的方法可以做到这一点？

【参考方案1】：

一个可能的优化是

>>> numpy.sum(a @ b, axis=0)
array([  1.83633615,  18.71643672,  15.26981078, -46.33670382,  13.30276476])
>>> numpy.sum(a, axis=0) @ b
array([  1.83633615,  18.71643672,  15.26981078, -46.33670382,  13.30276476])

计算 a @ b 需要 10k×200×10k 次运算，而先对行求和会将乘法运算减少到 1×200×10k 次运算，从而提高 10k×。

这主要是因为认识

   numpy.sum(x, axis=0) == [1, 1, ..., 1] @ x
=> numpy.sum(a @ b, axis=0) == [1, 1, ..., 1] @ (a @ b)
                            == ([1, 1, ..., 1] @ a) @ b
                            == numpy.sum(a, axis=0) @ b

其他轴也类似。

>>> numpy.sum(a @ b, axis=1)
array([  2.8794171 ,   9.12128399,  14.52009991,  -8.70177811, -15.0303783 ])
>>> a @ numpy.sum(b, axis=1)
array([  2.8794171 ,   9.12128399,  14.52009991,  -8.70177811, -15.0303783 ])

（注意：x @ y 等同于x.dot(y) for 2D matrixes and 1D vectors on Python 3.5+ with numpy 1.10.0+）

---

$ INITIALIZATION='import numpy;numpy.random.seed(0);a=numpy.random.randn(1000,200);b=numpy.random.rand(200,1000)'

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.einsum("ij,jk->k", a, b)'
10 loops, best of 3: 87.2 msec per loop

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.sum(a@b, axis=0)'
100 loops, best of 3: 12.8 msec per loop

$ python3 -m timeit -s "$INITIALIZATION" 'numpy.sum(a, axis=0)@b'
1000 loops, best of 3: 300 usec per loop

---

插图：

In [235]: a = np.random.rand(3,3)
array([[ 0.465,  0.758,  0.641],
       [ 0.897,  0.673,  0.742],
       [ 0.763,  0.274,  0.485]])

In [237]: b = np.random.rand(3,2)
array([[ 0.303,  0.378],
       [ 0.039,  0.095],
       [ 0.192,  0.668]])

现在，如果我们只做a @ b，我们将需要18 个乘法和6 个加法操作。另一方面，如果我们执行np.sum(a, axis=0) @ b，我们只需要6 个乘法和2 个加法操作。提高了 3 倍，因为我们在 a 中有 3 行。至于 OP 的情况，这应该比简单的 a @ b 计算提高 10k 倍，因为他在 a 中有 10k 行。

【讨论】：

非常感谢，运行大约需要 0.026 秒

只是对您的注释的小评论：@ 运算符并不完全是 python 3 的功能，因为它仅在 3.5 release 中引入。

@mgc 谢谢。已更新。

另外，还有一点需要注意：您的第二种方法只有在我们（最终）想做sum 和axis=0 时才有可能

【参考方案2】：

发生了两个sum-reductions - 一个来自与np.dot 的marix-multilication，然后是显式的sum。

我们可以使用np.einsum 一次性完成这两项工作，就像这样 -

np.einsum('ij,jk->k',a,b)

示例运行 -

In [27]: a = np.random.rand(3,4)

In [28]: b = np.random.rand(4,3)

In [29]: np.sum(a.dot(b), axis = 0)
Out[29]: array([ 2.70084316,  3.07448582,  3.28690401])

In [30]: np.einsum('ij,jk->k',a,b)
Out[30]: array([ 2.70084316,  3.07448582,  3.28690401])

运行时测试-

In [45]: a = np.random.rand(1000,200)

In [46]: b = np.random.rand(200,1000)

In [47]: %timeit np.sum(a.dot(b), axis = 0)
100 loops, best of 3: 5.5 ms per loop

In [48]: %timeit np.einsum('ij,jk->k',a,b)
10 loops, best of 3: 71.8 ms per loop

遗憾的是，看起来我们在使用 np.einsum 时并没有做得更好。

---

要更改为np.sum(a.dot(b), axis = 1)，只需在此处交换输出字符串表示法-np.einsum('ij,jk->i',a,b)，就像这样-

In [42]: np.sum(a.dot(b), axis = 1)
Out[42]: array([ 3.97805141,  3.2249661 ,  1.85921549])

In [43]: np.einsum('ij,jk->i',a,b)
Out[43]: array([ 3.97805141,  3.2249661 ,  1.85921549])

【讨论】：

谢谢，如果我想为axis = 1制作它我会改变什么？

我试过了，它已经运行了 5 分钟多，现在我猜它对于大矩阵来说很慢

@AyaAbdelsalam 是的！这就是我发现的。

ij,jk->k 证明了成功的解决方案。我们可以先对i 求和，然后对j 点求和。我想知道np.einsum('j,jk->k', np.einsum('ij->j', a), b) 会怎么做。

【参考方案3】：

使用我添加到 Divakar 答案中的想法进行一些快速测试：

In [162]: a = np.random.rand(1000,200)
In [163]: b = np.random.rand(200,1000)

In [174]: timeit c1=np.sum(a.dot(b), axis=0)
10 loops, best of 3: 27.7 ms per loop

In [175]: timeit c2=np.sum(a,axis=0).dot(b)
1000 loops, best of 3: 432 µs per loop

In [176]: timeit c3=np.einsum('ij,jk->k',a,b)
10 loops, best of 3: 170 ms per loop

In [177]: timeit c4=np.einsum('j,jk->k', np.einsum('ij->j', a), b)
1000 loops, best of 3: 353 µs per loop

In [178]: timeit np.einsum('ij->j', a) @b
1000 loops, best of 3: 304 µs per loop

einsum 实际上比np.sum 快！

In [180]: timeit np.einsum('ij->j', a)
1000 loops, best of 3: 173 µs per loop
In [181]: timeit np.sum(a,0)
1000 loops, best of 3: 312 µs per loop

对于较大的数组，einsum 的优势会降低

In [183]: a = np.random.rand(100000,200)
In [184]: b = np.random.rand(200,100000)
In [185]: timeit np.einsum('ij->j', a) @b
10 loops, best of 3: 51.5 ms per loop
In [186]: timeit c2=np.sum(a,axis=0).dot(b)
10 loops, best of 3: 59.5 ms per loop

【讨论】：

对于你的第一组，我得到了np.einsum('ij->j', a)：96.8 µs和np.sum(a,0)：74.8 µs。我在NumPy 1.12.0。我想我记得看到einsum 比sum 快，但是如果我没记错的话，在以前版本的 NumPy 上。

这里有相同的 numpy 版本。