在三维矩阵的情况下，用numpy操作替换for循环

Question

我有这个数学公式要实现！[https://i.stack.imgur.com/S28BA.png]其中例如w_fk表示形状(F，K(的矩阵。我已经将其实现为

gamma_dashed_lft = np.zeros((L, F, T))
for l in range(L):
for f in range(F):
for t in range(T):
temp = 0
for k in range(K):
temp = temp + (q_lk[l, k] * w_fk[f, k] * h_kt[k, t])
gamma_dashed_lft[l, f, t] = temp
return gamma_dashed_lft

在给定公式的情况下，用矩阵乘法代替循环的方法是什么？

Answer 1

您应该提供一个具体的例子。幸运的是，读取尺寸并创建并不太难

In [302]: L,F,T,K=2,3,4,5
In [303]: q_lk=np.arange(L*K).reshape(L,K)
In [304]: w_fk=np.arange(F*K).reshape(F,K)
In [305]: h_kt=np.arange(K*T).reshape(K,T)

当应用于您的代码时会产生：

In [306]: gamma_dashed_lft = np.zeros((L, F, T))
...: for l in range(L):
...:     for f in range(F):
...:         for t in range(T):
...:             temp = 0
...:             for k in range(K):
...:                 temp = temp + (q_lk[l, k] * w_fk[f, k] * h_kt[k, t])
...:             gamma_dashed_lft[l, f, t] = temp
...: 
In [308]: gamma_dashed_lft
Out[308]: 
array([[[ 400.,  430.,  460.,  490.],
[1000., 1080., 1160., 1240.],
[1600., 1730., 1860., 1990.]],
[[1000., 1080., 1160., 1240.],
[2600., 2855., 3110., 3365.],
[4200., 4630., 5060., 5490.]]])

充分利用broadcasting的等价表达式是：

In [309]: arr =(q_lk[:,None,:,None]*w_fk[None,:,:,None]*h_kt[None,None,:,:]).sum(axis=2)
In [310]: arr.shape
Out[310]: (2, 3, 4)
In [311]: np.allclose(arr,gamma_dashed_lft)
Out[311]: True

在设置广播时，我的目标是一个形状为(L，F，K，T(的阵列，其总和在K上减少。

既然你让我创建测试用例，我就让你计算广播细节。这对你来说将是一次很好的锻炼。

einsum

In [446]: D=np.einsum('lk,fk,kt->lft', q_lk, w_fk, h_kt)
In [447]: D.shape
Out[447]: (2, 3, 4)
In [448]: arr =(q_lk[:,None,:,None]*w_fk[None,:,:,None]*h_kt[None,None,:,:]).sum
...: (axis=2)
In [449]: np.allclose(arr,D)
Out[449]: True
In [450]: timeit arr =(q_lk[:,None,:,None]*w_fk[None,:,:,None]*h_kt[None,None,:,
...: :]).sum(axis=2)
22.4 µs ± 2.02 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [451]: timeit D=np.einsum('lk,fk,kt->lft', q_lk, w_fk, h_kt)
12.2 µs ± 40.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

matmul

In [458]: timeit E=((q_lk[:,None,None,:]*w_fk[None,:,None,:])@h_kt[None,None,:,:,: ]).squeeze()
10.6 µs ± 44.5 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

为此，我用(1,1，K，T(对@进行(L，F，1，K(数组，得到(L，F，1，T(。LF是matmul的"批次"维度，而K是产品维度的总和。

Answer 2

产品q[l,k]*w[f,k]是针对t的每个值进行评估的，所以我将其处理为：

for l,k:
qw[k] = q[l,k] * w[f,k]

然后在内部循环中使用该临时阵列CCD_ 8。现在有了一堆内积，实际上是一个矩阵向量积，并且在运算中保存了T因子。

然而，这并不能为您提供矩阵矩阵乘法所提供的缓存优化。为此，为每个l:创建一个矩阵qw_l

for l:
qw_l[f,k] = q[l,k]*w[f,k]
matrix-matrix-product: qw_l times h

(注意，我并没有拼写出所有的循环！(

这会额外花费一个F x K大小的数组，但这可能不是问题。如果你想知道，创建它的成本是F x K，而矩阵矩阵乘积是F x T x K，所以你不必担心。