在 Python（和 Cython）中计算两个矩阵的点积的最快方法是什么
What is the fastest way to compute the dot product of two matrices in Python (and Cython)
我正在尝试用 Cython 创建一个 Python 库，我需要在其中实现点积。我有一个计算点积的非常简单的方法，但是，对于较大的矩阵，它 运行 不够快。
我花了很多时间在谷歌上搜索这个问题，并试图让它尽快解决，但我无法让它更快。
下面的代码显示了我当前如何计算它的Python实现：
a = [[1, 2, 3], [4, 5, 6]]
b = [[1], [2], [3]]
def dot(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
for i in range(len(c)):
for j in range(len(c[i])):
t = 0
for k in range(len(b)):
t += a[i][k] * b[k][j]
c[i][j] = t
return c
print(dot(a, b))
# [[14], [32]]
这确实给出了正确的计算结果 (python [[14], [32]])，但是，对于我要使用它的用途来说，计算时间太长了。非常感谢任何有关如何使它更快的帮助。谢谢
您可以为此使用 numpy。 Numpy 实现了 BLAS 规范（基本线性代数子程序），它们是线性代数库的低级例程（如矩阵乘法）的事实标准。要获得两个矩阵的点积，例如 A 和 B，您可以使用以下代码：
A = [[1, 2, 3], [4, 5, 6]]
B = [[1], [2], [3]]
import numpy as np #Import numpy
numpy_a = np.array(A) #Cast your nested lists to numpy arrays
numpy_b = np.array(B)
print(np.dot(numpy_a, numpy_b)) #Print the result
根据结构的索引成本，您可以通过分解一些操作来提高速度：
def dot(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
bt = transpose(b) # can this be done once cheaply?
for i in range(len(c)):
a1 = a[i]
c1 = c[i]
for j in range(len(c1)):
b1 = bt[j]
t = 0
for k in range(len(b)):
t += a1[k] * b1[k]
c1[j] = t
return c
内部k循环可以写成，惯用的Python为：
for a2, b2 in zip(a1, b1):
t += a2 * b2
我不知道这在 cython 翻译中是否更快。
Fast cython还需要定义各种变量如int、float等，所以它可以直接c翻译，而不是经过一般的，但昂贵的，Python 对象。我不会尝试重复 cython 文档。
您应该注释（即，静态类型）所有可能的变量。如果您愿意，以下是我的解决方案：
# mydot.pyx
import numpy as np
cimport cython
def dot_1(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
for i in range(len(c)):
for j in range(len(c[i])):
t = 0
for k in range(len(b)):
t += a[i][k] * b[k][j]
c[i][j] = t
return c
@cython.boundscheck(False) # turn off bounds-checking
@cython.wraparound(False) # turn off negative index wrapping
def dot_2(double[:, :] A, double[:, :] B):
cdef Py_ssize_t M = A.shape[0]
cdef Py_ssize_t Na = A.shape[1]
cdef Py_ssize_t Nb = B.shape[0]
cdef Py_ssize_t K = B.shape[1]
assert Na == Nb
result = np.empty((M, K), dtype='d')
cdef double[:, :] C = result
cdef double t
for m in range(M):
for k in range(K):
t = 0
for n in range(Na):
t += A[m, n] * B[n, k]
C[m, k] = t
return result
和
# app.py
import pyximport
from numpy import array
from scipy import median
from timeit import repeat
pyximport.install()
from mydot import dot_1, dot_2
a = array([[1, 2, 3], [4, 5, 6]], dtype='d')
b = array([[1], [2], [3]], dtype='d')
dot_1_t = repeat('dot_1(a, b)', repeat=1000, number=1, globals=globals())
dot_2_t = repeat('dot_2(a, b)', repeat=1000, number=1, globals=globals())
print(f'dot_1 took {median(dot_1_t)*1000} ms.')
print(f'dot_2 took {median(dot_2_t)*1000} ms.')
当你运行cython --annotate mydot.pyx时，Cython会生成一个HTML文件注释Cython代码。在那里，黄色突出显示的颜色越深，生成的 C 代码的开销就越大 (Python)。您可以将两个解决方案（尤其是 for 循环）相互比较。
运行 python app.py 应该也会给你更快的结果。当然，如果您提供低于某个阈值的较小尺寸输入，您将看不到两者之间有意义的速度差异，因为您没有进行足够的迭代。然而，在一些阈值之后，速度差异应该是显着的，因为循环中的每次迭代对于您的版本来说都是昂贵的（参见较深的黄线）。
最后要说的是，正如这个问题下的每个人都已经建议的那样，当您提供具有更大维度的矩阵时，numpy 的函数应该具有更高的性能 --- 他们正在使用 blocked (sub)来自底层 BLAS 和 LAPACK 实现的矩阵运算，而不是天真地逐个迭代索引。
P.S：如果你想专注于 dot_2，不仅要专注于 double，还要专注于其他有意义的算术类型，例如 int 和 float s，你应该检查 Cython 的 fused types.
EDIT.因为我的回答后来被选为答案，所以我想举一个更大尺寸输入的例子。如果不使用上面的 app.py，而是使用以下内容：
# app.py
import pyximport
from numpy import array, random as rnd
from scipy import median
from timeit import repeat
pyximport.install()
from mydot import dot_1, dot_2
M = 100
N = 100
K = 1
a = rnd.randn(M, N)
b = rnd.randn(N, K)
dot_1_t = repeat('dot_1(a, b)', repeat=1000, number=1, globals=globals())
dot_2_t = repeat('dot_2(a, b)', repeat=1000, number=1, globals=globals())
print(f'dot_1 took {median(dot_1_t)*1000} ms.')
print(f'dot_2 took {median(dot_2_t)*1000} ms.')
计时应该类似于以下内容：
dot_1 took 5.218300502747297 ms.
dot_2 took 0.013017997844144702 ms.
我正在尝试用 Cython 创建一个 Python 库，我需要在其中实现点积。我有一个计算点积的非常简单的方法，但是，对于较大的矩阵，它 运行 不够快。
我花了很多时间在谷歌上搜索这个问题，并试图让它尽快解决，但我无法让它更快。
下面的代码显示了我当前如何计算它的Python实现：
a = [[1, 2, 3], [4, 5, 6]]
b = [[1], [2], [3]]
def dot(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
for i in range(len(c)):
for j in range(len(c[i])):
t = 0
for k in range(len(b)):
t += a[i][k] * b[k][j]
c[i][j] = t
return c
print(dot(a, b))
# [[14], [32]]
这确实给出了正确的计算结果 (python [[14], [32]])，但是，对于我要使用它的用途来说，计算时间太长了。非常感谢任何有关如何使它更快的帮助。谢谢
您可以为此使用 numpy。 Numpy 实现了 BLAS 规范（基本线性代数子程序），它们是线性代数库的低级例程（如矩阵乘法）的事实标准。要获得两个矩阵的点积，例如 A 和 B，您可以使用以下代码：
A = [[1, 2, 3], [4, 5, 6]]
B = [[1], [2], [3]]
import numpy as np #Import numpy
numpy_a = np.array(A) #Cast your nested lists to numpy arrays
numpy_b = np.array(B)
print(np.dot(numpy_a, numpy_b)) #Print the result
根据结构的索引成本，您可以通过分解一些操作来提高速度：
def dot(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
bt = transpose(b) # can this be done once cheaply?
for i in range(len(c)):
a1 = a[i]
c1 = c[i]
for j in range(len(c1)):
b1 = bt[j]
t = 0
for k in range(len(b)):
t += a1[k] * b1[k]
c1[j] = t
return c
内部k循环可以写成，惯用的Python为：
for a2, b2 in zip(a1, b1):
t += a2 * b2
我不知道这在 cython 翻译中是否更快。
Fast cython还需要定义各种变量如int、float等，所以它可以直接c翻译，而不是经过一般的，但昂贵的，Python 对象。我不会尝试重复 cython 文档。
您应该注释（即，静态类型）所有可能的变量。如果您愿意，以下是我的解决方案：
# mydot.pyx
import numpy as np
cimport cython
def dot_1(a, b):
c = [[0 for j in range(len(b[i]))] for i in range(len(a))]
for i in range(len(c)):
for j in range(len(c[i])):
t = 0
for k in range(len(b)):
t += a[i][k] * b[k][j]
c[i][j] = t
return c
@cython.boundscheck(False) # turn off bounds-checking
@cython.wraparound(False) # turn off negative index wrapping
def dot_2(double[:, :] A, double[:, :] B):
cdef Py_ssize_t M = A.shape[0]
cdef Py_ssize_t Na = A.shape[1]
cdef Py_ssize_t Nb = B.shape[0]
cdef Py_ssize_t K = B.shape[1]
assert Na == Nb
result = np.empty((M, K), dtype='d')
cdef double[:, :] C = result
cdef double t
for m in range(M):
for k in range(K):
t = 0
for n in range(Na):
t += A[m, n] * B[n, k]
C[m, k] = t
return result
和
# app.py
import pyximport
from numpy import array
from scipy import median
from timeit import repeat
pyximport.install()
from mydot import dot_1, dot_2
a = array([[1, 2, 3], [4, 5, 6]], dtype='d')
b = array([[1], [2], [3]], dtype='d')
dot_1_t = repeat('dot_1(a, b)', repeat=1000, number=1, globals=globals())
dot_2_t = repeat('dot_2(a, b)', repeat=1000, number=1, globals=globals())
print(f'dot_1 took {median(dot_1_t)*1000} ms.')
print(f'dot_2 took {median(dot_2_t)*1000} ms.')
当你运行cython --annotate mydot.pyx时，Cython会生成一个HTML文件注释Cython代码。在那里，黄色突出显示的颜色越深，生成的 C 代码的开销就越大 (Python)。您可以将两个解决方案（尤其是 for 循环）相互比较。
运行 python app.py 应该也会给你更快的结果。当然，如果您提供低于某个阈值的较小尺寸输入，您将看不到两者之间有意义的速度差异，因为您没有进行足够的迭代。然而，在一些阈值之后，速度差异应该是显着的，因为循环中的每次迭代对于您的版本来说都是昂贵的（参见较深的黄线）。
最后要说的是，正如这个问题下的每个人都已经建议的那样，当您提供具有更大维度的矩阵时，numpy 的函数应该具有更高的性能 --- 他们正在使用 blocked (sub)来自底层 BLAS 和 LAPACK 实现的矩阵运算，而不是天真地逐个迭代索引。
P.S：如果你想专注于 dot_2，不仅要专注于 double，还要专注于其他有意义的算术类型，例如 int 和 float s，你应该检查 Cython 的 fused types.
EDIT.因为我的回答后来被选为答案，所以我想举一个更大尺寸输入的例子。如果不使用上面的 app.py，而是使用以下内容：
# app.py
import pyximport
from numpy import array, random as rnd
from scipy import median
from timeit import repeat
pyximport.install()
from mydot import dot_1, dot_2
M = 100
N = 100
K = 1
a = rnd.randn(M, N)
b = rnd.randn(N, K)
dot_1_t = repeat('dot_1(a, b)', repeat=1000, number=1, globals=globals())
dot_2_t = repeat('dot_2(a, b)', repeat=1000, number=1, globals=globals())
print(f'dot_1 took {median(dot_1_t)*1000} ms.')
print(f'dot_2 took {median(dot_2_t)*1000} ms.')
计时应该类似于以下内容：
dot_1 took 5.218300502747297 ms.
dot_2 took 0.013017997844144702 ms.