"vectorize"在Numpy中进行大量矩阵运算最快的方法是什么?
What is the fastest way to "vectorize" a large number of matrix operations in Numpy?
假设我想在 Numpy 中做很多矩阵乘法;什么是最快的方法?
具体地说,这就是问题所在:我有两个很长的矩阵列表,我想将它们按元素相乘。也就是说,我有
[a_1, a_2, a_3, ..., a_N]
和
[b_1, b_2, b_3, ..., b_N],
其中每个 a_i
、b_i
是一个 nxn
矩阵(n
很小,比如 n=2
),而 N
是大(比如 N = 100000
),我想找到矩阵乘积 a_1 * b_1, a_2 * b_2, ...
使用 Python 和 Numpy/Scipy 最快的方法是什么?
一些选项是:
- 使用
for
循环——这很慢,因为 Python 循环很慢。
- 将小矩阵放入两个
NxN
块对角矩阵 A
和 B
--这将导致必须乘以一个比需要的大得多的矩阵。
- 使用
vectorize
-- 这是最容易编码的,但并不比 for
循环快。
只需像往常一样使用 numpy.matmul
或 @
,它是一个 ufunc 并且可以进行广播,只是不在数组元素上而是在矩阵子数组上。在您的情况下,您只需要将 (n,m) 矩阵堆叠在 (N,n,m) numpy 数组中,将 (m,p) 矩阵堆叠在 (N,m,p) numpy 数组中。
m = np.array([[1, 3], [2,4]])
m
Out[12]:
array([[1, 3],
[2, 4]])
m @ m
Out[13]:
array([[ 7, 15],
[10, 22]])
stackedm = np.stack([m,m,m])
stackedm
Out[15]:
array([[[1, 3],
[2, 4]],
[[1, 3],
[2, 4]],
[[1, 3],
[2, 4]]])
stackedm @ stackedm
Out[16]:
array([[[ 7, 15],
[10, 22]],
[[ 7, 15],
[10, 22]],
[[ 7, 15],
[10, 22]]])
您已经可以乘以 3D 数组,只需将您的数组列表放入 numpy 数组中,例如,
A = np.array([a_1, a_2, ..., a_N])
B = np.array([b_1, b_2, ..., b_N])
然后乘以A @ B
(@
是矩阵乘法运算符)。这是一个使用两个 3x3 数组“列表”的示例:
In [1]: import numpy as np
In [2]: x = np.random.randint(0, 9, (2, 3, 3))
In [3]: y = np.random.randint(0, 9, (2, 3, 3))
In [4]: x
Out[4]:
array([[[0, 4, 8],
[2, 5, 5],
[3, 0, 5]],
[[7, 6, 1],
[7, 0, 7],
[5, 2, 8]]])
In [5]: y
Out[5]:
array([[[7, 2, 6],
[6, 1, 4],
[6, 8, 5]],
[[8, 5, 4],
[8, 2, 7],
[3, 7, 0]]])
In [7]: x @ y
Out[7]:
array([[[ 72, 68, 56],
[ 74, 49, 57],
[ 51, 46, 43]],
[[107, 54, 70],
[ 77, 84, 28],
[ 80, 85, 34]]])
为了证明这一切所做的是每个矩阵在相应索引处的乘积:
In [8]: x[0]
Out[8]:
array([[0, 4, 8],
[2, 5, 5],
[3, 0, 5]])
In [9]: y[0]
Out[9]:
array([[7, 2, 6],
[6, 1, 4],
[6, 8, 5]])
In [10]: x[0] @ y[0]
Out[10]:
array([[72, 68, 56],
[74, 49, 57],
[51, 46, 43]])
In [11]: (x @ y)[0]
Out[11]:
array([[72, 68, 56],
[74, 49, 57],
[51, 46, 43]])
假设我想在 Numpy 中做很多矩阵乘法;什么是最快的方法?
具体地说,这就是问题所在:我有两个很长的矩阵列表,我想将它们按元素相乘。也就是说,我有
[a_1, a_2, a_3, ..., a_N]
和
[b_1, b_2, b_3, ..., b_N],
其中每个 a_i
、b_i
是一个 nxn
矩阵(n
很小,比如 n=2
),而 N
是大(比如 N = 100000
),我想找到矩阵乘积 a_1 * b_1, a_2 * b_2, ...
使用 Python 和 Numpy/Scipy 最快的方法是什么?
一些选项是:
- 使用
for
循环——这很慢,因为 Python 循环很慢。 - 将小矩阵放入两个
NxN
块对角矩阵A
和B
--这将导致必须乘以一个比需要的大得多的矩阵。 - 使用
vectorize
-- 这是最容易编码的,但并不比for
循环快。
只需像往常一样使用 numpy.matmul
或 @
,它是一个 ufunc 并且可以进行广播,只是不在数组元素上而是在矩阵子数组上。在您的情况下,您只需要将 (n,m) 矩阵堆叠在 (N,n,m) numpy 数组中,将 (m,p) 矩阵堆叠在 (N,m,p) numpy 数组中。
m = np.array([[1, 3], [2,4]])
m
Out[12]:
array([[1, 3],
[2, 4]])
m @ m
Out[13]:
array([[ 7, 15],
[10, 22]])
stackedm = np.stack([m,m,m])
stackedm
Out[15]:
array([[[1, 3],
[2, 4]],
[[1, 3],
[2, 4]],
[[1, 3],
[2, 4]]])
stackedm @ stackedm
Out[16]:
array([[[ 7, 15],
[10, 22]],
[[ 7, 15],
[10, 22]],
[[ 7, 15],
[10, 22]]])
您已经可以乘以 3D 数组,只需将您的数组列表放入 numpy 数组中,例如,
A = np.array([a_1, a_2, ..., a_N])
B = np.array([b_1, b_2, ..., b_N])
然后乘以A @ B
(@
是矩阵乘法运算符)。这是一个使用两个 3x3 数组“列表”的示例:
In [1]: import numpy as np
In [2]: x = np.random.randint(0, 9, (2, 3, 3))
In [3]: y = np.random.randint(0, 9, (2, 3, 3))
In [4]: x
Out[4]:
array([[[0, 4, 8],
[2, 5, 5],
[3, 0, 5]],
[[7, 6, 1],
[7, 0, 7],
[5, 2, 8]]])
In [5]: y
Out[5]:
array([[[7, 2, 6],
[6, 1, 4],
[6, 8, 5]],
[[8, 5, 4],
[8, 2, 7],
[3, 7, 0]]])
In [7]: x @ y
Out[7]:
array([[[ 72, 68, 56],
[ 74, 49, 57],
[ 51, 46, 43]],
[[107, 54, 70],
[ 77, 84, 28],
[ 80, 85, 34]]])
为了证明这一切所做的是每个矩阵在相应索引处的乘积:
In [8]: x[0]
Out[8]:
array([[0, 4, 8],
[2, 5, 5],
[3, 0, 5]])
In [9]: y[0]
Out[9]:
array([[7, 2, 6],
[6, 1, 4],
[6, 8, 5]])
In [10]: x[0] @ y[0]
Out[10]:
array([[72, 68, 56],
[74, 49, 57],
[51, 46, 43]])
In [11]: (x @ y)[0]
Out[11]:
array([[72, 68, 56],
[74, 49, 57],
[51, 46, 43]])