最后一次更新日期： 2019/3/21

此篇文章提供一个将scipy库的稀疏矩阵运用于卷积运算上，以期在不增加过多内存开销的情况下提高性能的入门思路。

先导入以下模块：
import numpy as np
from scipy import sparse
import time

1.基于滑窗的卷积运算

卷积核

卷积核的三个轴分别对应：高，宽，颜色通道

In [5]: kernel=np.eye(3).repeat(3).reshape((3,3,3))

In [6]: kernel[:,:,0]
Out[6]: 
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

图片数据

以cifar10数据集作为示例，加载数据集的方法此处不作讲解。
数据集的四个轴分别对应：样本，高，宽，颜色通道

In [11]: cifar=CifarManager()
    ...: train_images,train_labels,test_images,test_labels=cifar.read_as_array(chinese_label=True)
    ...: The default path of image data is set to:
    ...: D:\training_data\used\cifar-10-batches-py
    ...: reading data---
    ...: completed

In [12]: test_images.shape
Out[12]: (10000, 32, 32, 3)

In [13]: Image.fromarray(test_images[0].astype('uint8')).resize((200,200))
Out[13]: 

图片数据集的数据结构如下图所示：

单图卷积运算

#array: 单张图片数据
#kernel: 卷积核
#step: 步长
def conv1(array,kernel,step=1):
    #图片的高和宽
    h,w=array.shape[:2]
    #卷积核的高和宽
    kh,kw=kernel.shape[:2]
    #纵向和横向的移动步数
    hs,ws=(h-kh)//step+1,(w-kw)//step+1
    #初始化输出数组
    out=np.empty((hs,ws))
    #纵向滑动卷积核
    for i in range(hs):
        #卷积核的纵向覆盖范围
        start0,end0=i*step,i*step+kh
        #横向滑动卷积核
        for j in range(ws):
            #卷积核的横向覆盖范围
            start1,end1=j*step,j*step+kw
            #该位置的卷积运算
            out[i,j]=np.vdot(array[start0:end0,start1:end1],kernel)
    return out

In [37]: out=conv1(test_images[0],kernel)

In [38]: out.shape
Out[38]: (30, 30)

In [39]: out_=(out-out.min())/(out.max()-out.min())*255
    ...: Image.fromarray(out_.astype('uint8')).resize((200,200))
Out[39]: 

np.vdot用于将数组展开为向量后进行点积运算。

卷积运算最为简单的实现方式是滑窗，用卷积核在每个位置上覆盖图像数据，对应元素相乘再求和，存入结果数组的相对位置。

批量卷积运算

#array: 图片数据集
#kernel: 卷积核
#step: 步长
def conv2(array,kernel,step=1):
    #图片的张数、高、宽
    n,h,w=array.shape[:3]
    #卷积核的高、宽、颜色通道
    kh,kw,kc=kernel.shape[:3]
    #纵向和横向的移动步数
    hs,ws=(h-kh)//step+1,(w-kw)//step+1
    #初始化输出数组
    out=np.empty((n,hs,ws))
    #纵向滑动卷积核
    for i in range(hs):
        #卷积核的纵向覆盖范围
        start0,end0=i*step,i*step+kh
        #横向滑动卷积核
        for j in range(ws):
            #卷积核的横向覆盖范围
            start1,end1=j*step,j*step+kw
            #该位置的卷积运算
            array_=array[:,start0:end0,start1:end1].reshape((-1,kh*kw*kc))
            kernel_=kernel.ravel()
            out[:,i,j]=np.dot(array_,kernel_)
    return out

In [50]: start=time.clock()
    ...: out2=conv2(test_images,kernel)
    ...: print('\ntime used: %f s'%(time.clock()-start))

time used: 2.255798 s

In [51]: out2.shape
Out[51]: (10000, 30, 30)

In [52]: (out==out2[0]).all()
Out[52]: True

np.dot可以用于完成 向量点积，矩阵与向量的乘积，矩阵乘法 等运算。

在进行某一位置的卷积运算时，先在数据集上提取卷积核覆盖范围的数据，如上方的例子，可以得到一个形状为(10000,3,3,3)的张量，第一个轴对应每张图片，后三个轴对应卷积核的三个轴。将子数据集的后三个轴展开，形状变为(10000,27)，卷积核展开为向量，形状变为(27,)，然后进行一次矩阵和列向量的乘积运算就得到卷积结果了。

（为方便描述，公式中皆以1为起始索引）
以位置为例，n,h,w,i,j分别对应样本数,图片的高,图片的宽,纵向滑动的步数,横向滑动的步数，计算的主要过程如下：

2.基于权重矩阵的卷积运算

def conv3(array,kernel,step=1):
    #图片的张数、高、宽
    n,h,w=array.shape[:3]
    #卷积核的高、宽、颜色通道
    kh,kw,kc=kernel.shape[:3]
    #纵向和横向的移动步数
    hs,ws=(h-kh)//step+1,(w-kw)//step+1
    #滑窗卷积运算等效的权重矩阵
    weights=np.zeros((h*w*3,hs*ws))
    #纵向位置变化
    for i in range(hs):
        #卷积核的纵向覆盖范围
        start0,end0=i*step,i*step+kh
        #横向位置变化
        for j in range(ws):
            #卷积核的横向覆盖范围
            start1,end1=j*step,j*step+kw
            #当前位置的权重向量
            weights_=np.zeros((h,w,3))
            weights_[start0:end0,start1:end1]=kernel
            #合并至矩阵
            weights[:,i*ws+j]=weights_.ravel()
    #图片数据集变形
    array_=array.reshape((n,h*w*3))
    #矩阵乘法计算卷积结果
    result=np.dot(array_,weights)
    return result.reshape((-1,hs,ws)),weights

In [18]: start=time.clock()
    ...: out3,weights3=conv3(test_images,kernel)
    ...: print('\ntime used: %f s'%(time.clock()-start))

time used: 0.560485 s

In [19]: (out2==out3).all()
Out[19]: True

In [20]: weights.shape
Out[20]: (3072, 900)

In [21]: weights3.size
Out[21]: 2764800

可以看到，在消去循环体，将原本基于滑窗的卷积运算转换为权重矩阵运算后，性能得到了明显的提升(约3倍)。
但随之而来的问题就是，构造权重矩阵空间开销过大，权重矩阵形状为(输入总大小，输出总大小)，即使在低分辨率的cifar数据集上，该矩阵也达到了2764800个元素的大小，随着图片分辨率的提高，开销会增大到无法接受，乃至有性能倒退的可能。
考虑到卷积运算的权重矩阵有大量的零元素，这些元素对计算没有贡献，应当采用稀疏矩阵的形式进行处理。

对基于滑窗的计算方式进行转换的思路，以起始位置为例：

3.稀疏矩阵的定义

列压缩存储

In [54]: mt=np.array([[1,2,0,0],[0,4,0,0],[0,0,0,3],[0,0,0,0]])

In [55]: mt
Out[55]: 
array([[1, 2, 0, 0],
       [0, 4, 0, 0],
       [0, 0, 0, 3],
       [0, 0, 0, 0]])

In [59]: cscm=sparse.csc_matrix(mt)

In [60]: cscm.data
Out[60]: array([1, 2, 4, 3], dtype=int32)

In [61]: cscm.indices
Out[61]: array([0, 0, 1, 2], dtype=int32)

In [62]: cscm.indptr
Out[62]: array([0, 1, 3, 3, 4], dtype=int32)

In [63]: cscm
Out[63]: 
<4x4 sparse matrix of type '<class 'numpy.int32'>'
    with 4 stored elements in Compressed Sparse Column format>

稀疏矩阵用于在矩阵中零元素占多数的情况下减少存储和运算开销，零元素越多，提升越大。
稀疏矩阵有多种压缩数据的方式，本文中因为会将稀疏矩阵用于右乘，该运算下是按列提取向量，故采用列压缩存储。
可通过传入一个稠密矩阵创建对应的稀疏矩阵，但在原矩阵构造出来会很庞大时不建议这么做，可通过直接构造稀疏矩阵的关键参数来创建。

csc_matrix有三个关键参数：
(1). data参数是存储非零元素的值的一维数组，元素按先列后行排列；
(2). indices参数是存储非零元素的列坐标的一维数组，比如一个元素处在某一列的第一个位置，对应的列坐标就是0；
(3). indptr参数是存储每列的首个非零元素在data中位置的一维数组，再额外加上data的长度；当某一列没有非零元素时，对应的indptr元素和上一列相同，如果是第一列则为0。

csc_matrix的toarray和todense方法分别可以将稀疏矩阵还原为稠密矩阵的ndarray形式和matrix形式。

4.基于稀疏矩阵的卷积运算

由于先构造稠密矩阵再转换为稀疏矩阵的方式与使用稀疏矩阵减少存储开销的目的相违背，此处会通过构造关键参数直接创建稀疏矩阵。

data参数的构造较为简单，权重矩阵中的非零元素就是卷积核的元素，每个位置卷积运算的权重矩阵会被展开为列向量，csc_matrix中data的元素又是按先列后行排列的，所以只要将卷积核展开成的向量按位置总数重复拼接就行了；

indices参数的构造需要了解numpy中元素的排列顺序：
首先看起始位置，权重矩阵展开为列向量后，原本卷积核元素的位置如下所示：

indptr参数的构造也很简单，因为每列只有卷积核元素为非零元素(注意，即使卷积核元素可能为0，也要当作非零元素)，所以每列首个非零元素在data中的位置可以构成以卷积核大小为步长的等差数列。

def conv4(array,kernel,step=1):
    #图片的张数、高、宽
    n,h,w,c=array.shape
    #卷积核的高、宽、颜色通道
    kh,kw,kc=kernel.shape
    #纵向和横向的移动步数
    hs,ws=(h-kh)//step+1,(w-kw)//step+1
    #1.非零值
    data=np.tile(kernel.ravel(),hs*ws)
    #2.行索引
    #(即卷积核元素在对应的权重矩阵转换成的向量中的位置)
    #(1)起始位置索引
    #由三个轴方向上的索引偏移量组合得到
    idx0=np.arange(0,kh*w*c,w*c).reshape((-1,1,1))
    idx1=np.arange(0,kw*c,c).reshape((1,-1,1))
    idx2=np.arange(0,c,1).reshape((1,1,-1))
    loc_base_=idx0+idx1+idx2
    loc_base=np.tile(loc_base_.ravel(),hs*ws)
    #(2)横向和纵向移动的索引偏移
    loc_offset0=np.arange(0,hs*step*w*c,step*w*c).reshape((-1,1))
    loc_offset1=np.arange(0,ws*step*c,step*c).reshape((1,-1))
    loc_offset_=loc_offset0+loc_offset1
    loc_offset=loc_offset_.repeat(kh*kw*kc)
    indices=loc_base+loc_offset
    #3.列偏移
    #(即每列第一个非零值在values中的位置)
    indptr=np.arange(hs*ws+1)*(kh*kw*kc)
    #构造稀疏矩阵
    weights=sparse.csc_matrix((data,indices,indptr))
    #图片数据集变形
    array_=array.reshape((n,h*w*3))
    #稀疏矩阵乘法计算卷积结果
    result=array_*weights
    return result.reshape((-1,hs,ws)),weights

In [77]: start=time.clock()
    ...: out4,weights4=conv4(test_images,kernel)
    ...: print('\ntime used: %f s'%(time.clock()-start))

time used: 0.728124 s

In [77]: (out2==out4).all()
Out[77]: True

In [86]: weights4.data.size+weights4.indices.size+weights4.indptr.size
Out[86]: 49501

可以看到，虽然性能上稍有折损，但关键的内存开销降下来了，权重矩阵从原本的输入大小*输出大小个元素，减少到卷积核大小*输出大小*2+输出大小+1个元素。

在模型训练完毕后，可将所有卷积核转换为稀疏矩阵存储，避免每次前向传播都要重新构造。