最后一次更新日期： 2019/3/21
此篇文章提供一個(gè)將scipy庫(kù)的稀疏矩陣運(yùn)用于卷積運(yùn)算上，以期在不增加過多內(nèi)存開銷的情況下提高性能的入門思路。
先導(dǎo)入以下模塊：
import numpy as np
from scipy import sparse
import time
1.基于滑窗的卷積運(yùn)算
卷積核
卷積核的三個(gè)軸分別對(duì)應(yīng)：高，寬，顏色通道
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.]])
圖片數(shù)據(jù)
以cifar10數(shù)據(jù)集作為示例，加載數(shù)據(jù)集的方法此處不作講解。
數(shù)據(jù)集的四個(gè)軸分別對(duì)應(yīng)：樣本熟呛，高，寬尉姨，顏色通道
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]:
圖片數(shù)據(jù)集的數(shù)據(jù)結(jié)構(gòu)如下圖所示：
axis0對(duì)應(yīng)樣本庵朝，axis1對(duì)應(yīng)圖片的高，axis2對(duì)應(yīng)圖片的寬，axis3對(duì)應(yīng)顏色通道九府。從中抽取一張圖片椎瘟，則
axis0,axis1,axis2分別對(duì)應(yīng)高,寬,通道。
單圖卷積運(yùn)算
#array: 單張圖片數(shù)據(jù)
#kernel: 卷積核
#step: 步長(zhǎng)
def conv1(array,kernel,step=1):
#圖片的高和寬
h,w=array.shape[:2]
#卷積核的高和寬
kh,kw=kernel.shape[:2]
#縱向和橫向的移動(dòng)步數(shù)
hs,ws=(h-kh)//step+1,(w-kw)//step+1
#初始化輸出數(shù)組
out=np.empty((hs,ws))
#縱向滑動(dòng)卷積核
for i in range(hs):
#卷積核的縱向覆蓋范圍
start0,end0=i*step,i*step+kh
#橫向滑動(dòng)卷積核
for j in range(ws):
#卷積核的橫向覆蓋范圍
start1,end1=j*step,j*step+kw
#該位置的卷積運(yùn)算
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用于將數(shù)組展開為向量后進(jìn)行點(diǎn)積運(yùn)算侄旬。
卷積運(yùn)算最為簡(jiǎn)單的實(shí)現(xiàn)方式是滑窗肺蔚，用卷積核在每個(gè)位置上覆蓋圖像數(shù)據(jù)，對(duì)應(yīng)元素相乘再求和儡羔，存入結(jié)果數(shù)組的相對(duì)位置宣羊。
單張圖片，單個(gè)卷積核汰蜘，在某一個(gè)位置的卷積運(yùn)算如下圖所示：
批量卷積運(yùn)算
#array: 圖片數(shù)據(jù)集
#kernel: 卷積核
#step: 步長(zhǎng)
def conv2(array,kernel,step=1):
#圖片的張數(shù)仇冯、高、寬
n,h,w=array.shape[:3]
#卷積核的高族操、寬苛坚、顏色通道
kh,kw,kc=kernel.shape[:3]
#縱向和橫向的移動(dòng)步數(shù)
hs,ws=(h-kh)//step+1,(w-kw)//step+1
#初始化輸出數(shù)組
out=np.empty((n,hs,ws))
#縱向滑動(dòng)卷積核
for i in range(hs):
#卷積核的縱向覆蓋范圍
start0,end0=i*step,i*step+kh
#橫向滑動(dòng)卷積核
for j in range(ws):
#卷積核的橫向覆蓋范圍
start1,end1=j*step,j*step+kw
#該位置的卷積運(yùn)算
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可以用于完成 向量點(diǎn)積，矩陣與向量的乘積色难，矩陣乘法 等運(yùn)算泼舱。
在進(jìn)行某一位置的卷積運(yùn)算時(shí)，先在數(shù)據(jù)集上提取卷積核覆蓋范圍的數(shù)據(jù)枷莉，如上方的例子娇昙，可以得到一個(gè)形狀為(10000,3,3,3)的張量，第一個(gè)軸對(duì)應(yīng)每張圖片笤妙，后三個(gè)軸對(duì)應(yīng)卷積核的三個(gè)軸冒掌。將子數(shù)據(jù)集的后三個(gè)軸展開，形狀變?yōu)?code>(10000,27)危喉，卷積核展開為向量宋渔，形狀變?yōu)?code>(27,)州疾，然后進(jìn)行一次矩陣和列向量的乘積運(yùn)算就得到卷積結(jié)果了辜限。
（為方便描述，公式中皆以1為起始索引）
以位置為例严蓖，
n,h,w,i,j分別對(duì)應(yīng)樣本數(shù),圖片的高,圖片的寬,縱向滑動(dòng)的步數(shù),橫向滑動(dòng)的步數(shù)薄嫡，計(jì)算的主要過程如下：
2.基于權(quán)重矩陣的卷積運(yùn)算
def conv3(array,kernel,step=1):
#圖片的張數(shù)、高颗胡、寬
n,h,w=array.shape[:3]
#卷積核的高毫深、寬、顏色通道
kh,kw,kc=kernel.shape[:3]
#縱向和橫向的移動(dòng)步數(shù)
hs,ws=(h-kh)//step+1,(w-kw)//step+1
#滑窗卷積運(yùn)算等效的權(quán)重矩陣
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
#當(dāng)前位置的權(quán)重向量
weights_=np.zeros((h,w,3))
weights_[start0:end0,start1:end1]=kernel
#合并至矩陣
weights[:,i*ws+j]=weights_.ravel()
#圖片數(shù)據(jù)集變形
array_=array.reshape((n,h*w*3))
#矩陣乘法計(jì)算卷積結(jié)果
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
可以看到毒姨，在消去循環(huán)體哑蔫，將原本基于滑窗的卷積運(yùn)算轉(zhuǎn)換為權(quán)重矩陣運(yùn)算后，性能得到了明顯的提升(約3倍)。
但隨之而來的問題就是闸迷，構(gòu)造權(quán)重矩陣空間開銷過大嵌纲，權(quán)重矩陣形狀為(輸入總大小，輸出總大小)腥沽，即使在低分辨率的cifar數(shù)據(jù)集上逮走，該矩陣也達(dá)到了2764800個(gè)元素的大小，隨著圖片分辨率的提高今阳，開銷會(huì)增大到無法接受师溅，乃至有性能倒退的可能。
考慮到卷積運(yùn)算的權(quán)重矩陣有大量的零元素盾舌，這些元素對(duì)計(jì)算沒有貢獻(xiàn)墓臭，應(yīng)當(dāng)采用稀疏矩陣的形式進(jìn)行處理。
對(duì)基于滑窗的計(jì)算方式進(jìn)行轉(zhuǎn)換的思路矿筝，以起始位置為例：
3.稀疏矩陣的定義
列壓縮存儲(chǔ)
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>
稀疏矩陣用于在矩陣中零元素占多數(shù)的情況下減少存儲(chǔ)和運(yùn)算開銷起便，零元素越多，提升越大窖维。
稀疏矩陣有多種壓縮數(shù)據(jù)的方式榆综，本文中因?yàn)闀?huì)將稀疏矩陣用于右乘，該運(yùn)算下是按列提取向量铸史，故采用列壓縮存儲(chǔ)鼻疮。
可通過傳入一個(gè)稠密矩陣創(chuàng)建對(duì)應(yīng)的稀疏矩陣，但在原矩陣構(gòu)造出來會(huì)很龐大時(shí)不建議這么做琳轿，可通過直接構(gòu)造稀疏矩陣的關(guān)鍵參數(shù)來創(chuàng)建判沟。
csc_matrix有三個(gè)關(guān)鍵參數(shù)：
(1). data參數(shù)是存儲(chǔ)非零元素的值的一維數(shù)組，元素按先列后行排列崭篡；
(2). indices參數(shù)是存儲(chǔ)非零元素的列坐標(biāo)的一維數(shù)組挪哄，比如一個(gè)元素處在某一列的第一個(gè)位置，對(duì)應(yīng)的列坐標(biāo)就是0琉闪；
(3). indptr參數(shù)是存儲(chǔ)每列的首個(gè)非零元素在data中位置的一維數(shù)組迹炼，再額外加上data的長(zhǎng)度；當(dāng)某一列沒有非零元素時(shí)颠毙，對(duì)應(yīng)的indptr元素和上一列相同斯入，如果是第一列則為0。
csc_matrix的toarray和todense方法分別可以將稀疏矩陣還原為稠密矩陣的ndarray形式和matrix形式蛀蜜。
4.基于稀疏矩陣的卷積運(yùn)算
由于先構(gòu)造稠密矩陣再轉(zhuǎn)換為稀疏矩陣的方式與使用稀疏矩陣減少存儲(chǔ)開銷的目的相違背刻两，此處會(huì)通過構(gòu)造關(guān)鍵參數(shù)直接創(chuàng)建稀疏矩陣。
data參數(shù)的構(gòu)造較為簡(jiǎn)單滴某，權(quán)重矩陣中的非零元素就是卷積核的元素磅摹，每個(gè)位置卷積運(yùn)算的權(quán)重矩陣會(huì)被展開為列向量滋迈，csc_matrix中data的元素又是按先列后行排列的，所以只要將卷積核展開成的向量按位置總數(shù)重復(fù)拼接就行了户誓；
indices參數(shù)的構(gòu)造需要了解numpy中元素的排列順序：
首先看起始位置杀怠，權(quán)重矩陣展開為列向量后，原本卷積核元素的位置如下所示：
圖上的數(shù)字標(biāo)識(shí)的即是展開為向量后元素的位置厅克，numpy數(shù)組上對(duì)元素的訪問是從最后一個(gè)軸開始的赔退，在RGB圖片形狀
(h,w,c)下，2軸的索引每加1证舟，展開后的索引加1硕旗，1軸的索引每加1，展開后的索引加c女责，0軸索引每加1漆枚，展開后的索引加w*c；在橫向移動(dòng)的過程中抵知，每移動(dòng)一步墙基，展開后索引會(huì)全體增加
step*c；在縱向移動(dòng)的過程中刷喜，每移動(dòng)一步残制，展開后索引會(huì)全體增加
step*w*c。因此
indices可以分解為三部分：起始位置索引+橫向移動(dòng)帶來的索引偏移+縱向移動(dòng)帶來的索引偏移掖疮。
indptr參數(shù)的構(gòu)造也很簡(jiǎn)單初茶，因?yàn)槊苛兄挥芯矸e核元素為非零元素(注意，即使卷積核元素可能為0浊闪，也要當(dāng)作非零元素)恼布，所以每列首個(gè)非零元素在data中的位置可以構(gòu)成以卷積核大小為步長(zhǎng)的等差數(shù)列。
def conv4(array,kernel,step=1):
#圖片的張數(shù)搁宾、高折汞、寬
n,h,w,c=array.shape
#卷積核的高、寬盖腿、顏色通道
kh,kw,kc=kernel.shape
#縱向和橫向的移動(dòng)步數(shù)
hs,ws=(h-kh)//step+1,(w-kw)//step+1
#1.非零值
data=np.tile(kernel.ravel(),hs*ws)
#2.行索引
#(即卷積核元素在對(duì)應(yīng)的權(quán)重矩陣轉(zhuǎn)換成的向量中的位置)
#(1)起始位置索引
#由三個(gè)軸方向上的索引偏移量組合得到
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)橫向和縱向移動(dòng)的索引偏移
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.列偏移
#(即每列第一個(gè)非零值在values中的位置)
indptr=np.arange(hs*ws+1)*(kh*kw*kc)
#構(gòu)造稀疏矩陣
weights=sparse.csc_matrix((data,indices,indptr))
#圖片數(shù)據(jù)集變形
array_=array.reshape((n,h*w*3))
#稀疏矩陣乘法計(jì)算卷積結(jié)果
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
可以看到爽待，雖然性能上稍有折損，但關(guān)鍵的內(nèi)存開銷降下來了奸忽，權(quán)重矩陣從原本的輸入大小*輸出大小個(gè)元素堕伪，減少到卷積核大小*輸出大小*2+輸出大小+1個(gè)元素揖庄。
在模型訓(xùn)練完畢后栗菜，可將所有卷積核轉(zhuǎn)換為稀疏矩陣存儲(chǔ)，避免每次前向傳播都要重新構(gòu)造蹄梢。
