从卷积核生成线性运算矩阵的 Numpy 方法答案

【问题标题】：Numpy way to generate linear operation matrix from a convolution kernel从卷积核生成线性运算矩阵的 Numpy 方法
【发布时间】：2020-01-18 11:16:09
【问题描述】：

形状为(k1, k2, n_channel, n_filter) 的二维卷积核K 应用于形状为(m1, m2, n_channel) 的二维向量A 并生成另一个形状为(m1 - k1 + 1, m2 - k2 + 1, n_filter) 的二维向量B（带有有效填充）。

对于每个K，也确实存在一个W_K，其形状为(m1 - k1 + 1, m2 - k2 + 1, n_filter, m1, m2, n_channel)，因此W_K和A的张量点等于B。即B = np.tensordot(W_K, A, 3)。

我正在尝试找到一个纯 NumPy 解决方案来从 K 生成此 W_K 而不使用任何 python 循环。

我可以看到W_K[i,j,f] == np.pad(K[...,f], ((i,m1-i-k1), (j,m2-j-k2)), 'constant', constant_values=0) 或只是W_K[i, j, f, i:i+k1, j:j+k2, ...] == K[..., f]。

我正在寻找的几乎类似于 Toeplitz 矩阵。但我需要多维度的。

循环代码示例：

import numpy as np

# 5x5 image with 3-channels
A = np.random.random((5,5,3))
# 2x2 Conv2D kernel with 2 filters for A  
K = np.random.random((2,2,3,2))

# It should be of (4,4,2,5,5,3), but I create this way for convenience. I move the axis at the end.
W_K = np.empty((4,4,5,5,3,2))
for i, j in np.ndindex(4, 4):
  W_K[i, j] = np.pad(K, ((i, 5-i-2),(j, 5-j-2), (0, 0), (0, 0)), 'constant', constant_values=0)

# above lines can also be rewritten as
W_K = np.zeros((4,4,5,5,3,2))
for i, j in np.ndindex(4, 4):
  W_K[i, j, i:i+2, j:j+2, ...] = K[...]

W_K = np.moveaxis(W_K, -1, 2)

# now I can do
B = np.tensordot(W_K, A, 3)

【问题讨论】：

标签： python numpy tensorflow vectorization convolution

【解决方案1】：

你想要的东西需要一点fancy indexing 体操，但编写代码并不是很麻烦。这个想法是创建 4 维索引数组，应用您的第二个循环示例的 W_K[i, j, i:i+2, j:j+2, ...] 部分。

这是您示例的略微修改版本，只是为了确保某些相关维度不同（因为这使错误更容易找到：它们将是正确的错误而不是损坏的值）：

import numpy as np

# parameter setup
k1, k2, nch, nf = 2, 4, 3, 2 
m1, m2 = 5, 6 
w1, w2 = m1 - k1 + 1, m2 - k2 + 1 
K = np.random.random((k1, k2, nch, nf)) 
A = np.random.random((m1, m2, nch)) 

# your loopy version for comparison
W_K = np.zeros((w1, w2, nf, m1, m2, nch)) 
for i, j in np.ndindex(w1, w2): 
    W_K[i, j, :, i:i+k1, j:j+k2, ...] = K.transpose(-1, 0, 1, 2) 

W_K2 = np.zeros((w1, w2, m1, m2, nch, nf))  # to be transposed back
i,j = np.mgrid[:w1, :w2][..., None, None]  # shape (w1, w2, 1, 1) 
k,l = np.mgrid[:k1, :k2]  # shape (k1, k2) ~ (1, 1, k1, k2)  

W_K2[i, j, i+k, j+l, ...] = K 
W_K2 = np.moveaxis(W_K2, -1, 2) 

print(np.array_equal(W_K, W_K2))  # True

我们首先创建一个跨越W_K 的前两个维度的索引网格i,j，然后创建跨越其（moveaxis 之前的）第二个维度和第三个维度的两个相似网格。通过将两个尾随单例维度注入前者，我们最终得到 4d 索引数组，它们一起跨越 W_K 的前四个维度。

剩下的就是使用原始的K 分配给这个切片，然后移回维度。由于当表达式中的切片（非高级）索引并非全部相邻时，高级索引会如何改变行为，因此使用 moveaxis 方法更容易做到这一点。我首先尝试使用其最终尺寸创建W_K2，但随后我们将拥有具有细微不同行为（特别是不同形状）的W_K[i, j, :, i+k, j+l, ...]。

【讨论】：