转载自知乎
感谢 @Sebastian 提供的完整实现https://github.com/sebgao/cTensor
在之前的一篇笔记中推导了卷积运算的im2col方法的实现：https://zhuanlan.zhihu.com/p/63974249
numpy.lib.stride_tricks.as_strided(x, shape=None, strides=None, subok=False, writeable=True) X = np.arange(9, dtype=np.int32).reshape(3,3) print(X) ''' [[0 1 2] [3 4 5] [6 7 8]] ''' A = np.lib.stride_tricks.as_strided(X, shape=(2,2,2,2), strides=(12,4,12,4)) print(A) ''' [[[[0 1] [3 4]] [[1 2] [4 5]]] [[[3 4] [6 7]] [[4 5] [7 8]]]] '''X = np.arange(16, dtype=np.int32).reshape(4,4) print(X) ''' [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15]] ''' A = np.lib.stride_tricks.as_strided(X, shape=(2,2,3,3), strides=(16,4,16,4)) print(A) ''' [[[[ 0 1 2] [ 4 5 6] [ 8 9 10]] [[ 1 2 3] [ 5 6 7] [ 9 10 11]]] [[[ 4 5 6] [ 8 9 10] [12 13 14]] [[ 5 6 7] [ 9 10 11] [13 14 15]]]] '''
X1 = np.arange(1,19, dtype=np.int32).reshape(2,3,3) print(X1) ''' [[[ 1 2 3] [ 4 5 6] [ 7 8 9]] [[10 11 12] [13 14 15] [16 17 18]]] ''' A1 = as_strided(X1, shape=(2,2,2,2,2), strides=(*X1.strides, *X1.strides[-2:])) print(A1) ''' [[[[[ 1 2] [ 4 5]] [[ 2 3] [ 5 6]]] [[[ 4 5] [ 7 8]] [[ 5 6] [ 8 9]]]] [[[[10 11] [13 14]] [[11 12] [14 15]]] [[[13 14] [16 17]] [[14 15] [17 18]]]]] '''
完全正确
桶里，当X是四维数组时（X.shape = (N, C, H, W)），卷积核的宽高为：kh，kw
那么可以确定A的shape为：（N，C，？，？，kh，kw），至于中间的两个数值怎么确定？还是用以前的公式：
def split_by_strides(X, kh, kw, s): N, C, H, W = X.shape oh = (H - kh) // s + 1 ow = (W - kw) // s + 1 strides = (*X.strides[:-2], X.strides[-2]*s, X.strides[-1]*s, *X.strides[-2:]) A = as_strided(X, shape=(N,C,oh,ow,kh,kw), strides=strides)这里没有padding，因为在这之前就X就已经认为是pad过了
还有个小问题，如果X.shape不是上面那个样子怎么弄？这里举个其他的例子：
def split_by_strides(X, kh, kw, s): N, H, W, C = X.shape oh = (H - kh) // s + 1 ow = (W - kw) // s + 1 shape = (N, oh, ow, kh, kw, C) strides = (X.strides[0], X.strides[1]*s, X.strides[2]*s, *X.strides[1:]) A = as_strided(X, shape=shape, strides=strides) return A其实就是将shape改一下，strides参数调整一下就行了
现在蛋糕已经切好了，那么是直接吃吗？不过还有卷积核的问题没有解决，如果把A经过transpose和reshape转换成二维矩阵，再把卷积核如法炮制，那么直接用np.dot就可以实现了，得到结果之后再reshape回来也就是想要的feature map了
可是，这么一顿操作下来，好不容易去掉for循环提高的效率，不是又被拉回来了？有没有一行代码的解决方案？这时候就需要用到numpy的另一个函数tensordot了
不过这样就说来话长了，足够另起一篇了，而且单独一篇篇幅太长也不利于阅读，这里先把代码贴出来，然后回头再介绍tensordot
设卷积核的shape为：（kh，kw，C，kn），X.shape：（N，H，W，C），计算feature_map的代码如下：
kh, kw, C, kn = self.filters.shape X_split = split_by_strides(X, kh, kw, s) # X_split.shape: (N, oh, ow, kh, kw, C) feature_map = np.tensordot(X_split, self.filters, axes=[(3,4,5), (0,1,2)]) + self.bias可以检验一下，featrue_map的shape为：（N，oh，ow，kn）
最后就是上面代码中用到的np.tensordot函数了，这个和np.dot类似，不过np.dot是矩阵乘法，tensordot是张量的点乘，简单的写一下用法，先来看看官方手册上tensordot的介绍：
numpy.tensordot( a, b, axes=2)其中参数a，b都是np数组；axes是给定的轴，axes赋整数的情况不在讨论范围之内，这里讨论的是其为数组的情况，先看最简单的例子：
>>> A = np.random.randint(0,9,(3,4)) >>> B = np.random.randint(0,9,(4,5)) >>> C = np.tensordot(A, B, [1, 0]) >>> C.shape (3, 5) >>> C array([[ 37, 23, 77, 19, 16], [ 83, 32, 116, 40, 51], [ 86, 43, 127, 47, 47]]) >>> np.dot(A, B) array([[ 37, 23, 77, 19, 16], [ 83, 32, 116, 40, 51], [ 86, 43, 127, 47, 47]])这个和np.dot的结果虽然没有区别，不过为了说明其中的不同，还是用另一种角度来进行抽象，先将上面的两个二维数组分别抽象成一个列向量和一个行向量：
说了这么多，主要是为了接下来的内容做铺垫，再来看一个复杂点的例子：
>>> A = np.random.randint(0,9,(3,4,5)) >>> B = np.random.randint(0,9,(4,5,2)) >>> np.tensordot(A, B, [(1,2), (0,1)]) array([[290, 146], [308, 164], [316, 205]])经过上面的解释，想必你肯定知道这些数值是怎么来的了，来验证一下：
>>> np.sum(A[0] * B[:,:,0]) 290 >>> np.sum(A[2] * B[:,:,1]) 205上面分别验证了（0，0）和（2，1）位置的值，都完全正确
所以，tensordot的功能就是对于A, B两个多维数组，按给定的axes得到分别得到子数组subA和subB，然后这两个子数组进行点乘运算
官方手册上有个例子，这里搬过来，应该能更好的理解：
>>> a = np.arange(60.).reshape(3,4,5) >>> b = np.arange(24.).reshape(4,3,2) >>> c = np.tensordot(a,b, axes=([1,0],[0,1])) >>> c.shape (5, 2) >>> c array([[ 4400., 4730.], [ 4532., 4874.], [ 4664., 5018.], [ 4796., 5162.], [ 4928., 5306.]]) >>> # A slower but equivalent way of computing the same... 一个更慢但效果一样的方法 >>> d = np.zeros((5,2)) >>> for i in range(5): ... for j in range(2): ... for k in range(3): ... for n in range(4): ... d[i,j] += a[k,n,i] * b[n,k,j]