在使用scipy和numpy做数据计算时，感觉运行速度较慢，但是程序已经到了使用多数计算使用内积运算地步了，真的不知道该如何优化。如果能够优化下内积运算该有多好啊，奔着这个目标，希望能够写一篇文章盘点各种内积优化方法，也算是贡献自己的微薄之力。
开篇我写两点自己经验，抛砖引玉，希望大家多多提意见。由于自己对对于Scipy和Numpy熟悉度不够，所以有不正确的地方，还请大家多多斧正。
在说我的优化之前，先啰嗦下：scipy.sparse的矩阵包中，牵扯到矩阵运算，矩阵的格式优选csr_matrix和csc_matrix。不然速度肯定慢的你怀疑人生。
特别说明1本文的实验在ipython或者jupyter环境进行，时间消耗测试使用的是“%timeit”命令，Scipy版本为“0.19.1”。
特别说明2在程序中，很多非计算操作，如：list转稀疏矩阵、矩阵转置、矩阵拼接和矩阵更新等，由于它们具有内存操作，所以时间代价相当昂贵，并且可以提前处理，所以在测量时间消耗时，无需将他们的时间消耗也计算在内。在性能优化中，有两条原则相当重要：减少内存操作和减少CPU命令数。更多详情查看《Python高性能编程》第6章。
特别说明3如果你是计算专业的在读生，那么学好《计算机架构导论》、《操作系统》、《数据结构》、《离散数学》。前两本书让你在硬件和操作系统层次明白编程语言的特性，配上一些相关书籍，你会很快明白为什么会快，为什么会慢，为什么有些语言风格会快，有些则慢。后两本则告诉你如何优化你的算法，好比：现在从山北到山南，你可以从北山脚爬到山顶再到南山脚，也可以围着山跑，从北山脚跑到南山脚。当然，这些书的用处，绝不仅于此，它也是科班生与培训班生的区别。计算机编程不是学好几门编程语言和数据结构那么简单。
一、大小矩阵内积运算
当两个规模相当的矩阵做内积时，选择CSC或CSR并没有太大差别，时间效果相当。但是当为一大一小矩阵时，就有一些技巧，可以节约时间。假设B为大矩阵，S为小矩阵。
- 当CSR格式时，S×B速度较快，与B×S相比节约了一半时间。
- 当CSC格式时，B×S速度较快，与S×B相比节约一半时间。
上述两种方法，时间相近，不分伯仲之间。
以下是我的计算例子。
import scipy.sparse as sp
def is_csr_instance(mtx):
if isinstance(mtx, sp.csr_matrix):
return True
else:
return False
def is_csc_instance(mtx):
if isinstance(mtx, sp.csc_matrix):
return True
else:
return False
a_mtx = sp.csc_matrix([[1., 1., 3.]*120])
mtx = sp.csc_matrix([[1., 0., 0.]*120]*30000)
is_csc_instance(a_mtx), is_csc_instance(mtx)
mtx.shape, a_mtx.shape
mtx_T = mtx.T
mtx_T = mtx_T.tocsc()
print is_csc_instance(mtx_T), is_csr_instance(mtx_T)
print u"\n\ncsc little×big"
print type(a_mtx), type(mtx_T)
print a_mtx.shape, mtx_T.shape
%timeit c = a_mtx.dot(mtx_T)
print u"\n\ncsr little×big"
a_mtx_r = a_mtx.tocsr()
mtx_T_r = mtx_T.tocsr()
print type(a_mtx_r), type(mtx_T_r)
print a_mtx_r.shape, mtx_T_r.shape
%timeit c = a_mtx_r.dot(mtx_T_r)
a_mtx_T = a_mtx.T
a_mtx_T = a_mtx_T.tocsc()
mtx_T.shape, a_mtx_T.shape
print "\n\ncsc big×little"
print type(mtx), type(a_mtx_T)
print mtx.shape, a_mtx_T.shape
%timeit c = mtx.dot(a_mtx_T)
print "\n\ncsr big×little"
mtx = mtx.tocsr()
a_mtx_T = a_mtx_T.tocsr()
print type(mtx), type(a_mtx_T)
print mtx.shape, a_mtx_T.shape
%timeit c = mtx.dot(a_mtx_T)
输出如下：
csc little×big
<class 'scipy.sparse.csc.csc_matrix'> <class 'scipy.sparse.csc.csc_matrix'>
(1, 360) (360, 30000)
100 loops, best of 3: 17.4 ms per loop
csr little×big
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'>
(1, 360) (360, 30000)
100 loops, best of 3: 8.13 ms per loop
csc big×little
<class 'scipy.sparse.csc.csc_matrix'> <class 'scipy.sparse.csc.csc_matrix'>
(30000, 360) (360, 1)
100 loops, best of 3: 8.31 ms per loop
csr big×little
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'>
(30000, 360) (360, 1)
100 loops, best of 3: 17.6 ms per loop
二 多矩阵内积优化
不好意思，这条优化有时有效有时无效，所以暂时不要完全相信，欢迎各位对此条多提意见。
当有多个矩阵进行内积计算时，可以通过矩阵拼接将多次内积计算合并为一次节约时间。时间优化效果与矩阵的中需要计算的非零数据次数成反比，需要计算的次数越多，节约的时间越少。假设稀疏矩阵中，非零元素随机出现，那么需要计算的非零数据次数非常少，所以有近似结论：矩阵越稀疏，需要计算的非零数据越少，节约的时间越多。矩阵稠密度是非零元素个数与矩阵总元素数的比值。
本实验有两个组，对照组为一个1×N与一个M×N的矩阵做四次内积，实验组为一个1×4N的矩阵与一个M×4N的矩阵做一次内积。实验分3次：例1，例2和例3：
- 例1中，两个矩阵稠密度为100%，对照组时间消耗略高。
- 例2中，两个矩阵稠密度为33.34%，对照组时间较高。
- 例3中，两个矩阵稠密度分别为16.7%和8.3%，对照组时间消耗明显很高。
实验公共代码
import scipy.sparse as sp
def quadra_dot(a_mtx, b_mtx):
a = a_mtx * b_mtx
b = a_mtx * b_mtx
c = a_mtx * b_mtx
d = a_mtx * b_mtx
def uni_dot(a_mtx, b_mtx):
a = a_mtx * b_mtx
def density(mtx):
non_zeros_numbers = len(mtx.data) * 1.0
m, n = mtx.shape
print non_zeros_numbers / (m*n)
例1
a_mtx = sp.csr_matrix([[2.23, 1.56, 3.47]*120]*300)
mtx = sp.csr_matrix([[1.07, 2.19, 3.12]*120]*30000)
print(u"对照组：")
b_mtx = mtx.T
b_mtx = b_mtx.tocsr()
print type(a_mtx), type(b_mtx), a_mtx.shape, b_mtx.shape
# 测试时间消耗
%timeit quadra_dot(a_mtx, b_mtx)
print(u"实验组：")
c_mtx = sp.vstack((b_mtx, b_mtx))
c_mtx = sp.vstack((c_mtx, b_mtx))
c_mtx = sp.vstack((c_mtx, b_mtx))
a_mtx = sp.hstack((a_mtx, a_mtx))
a_mtx = sp.hstack((a_mtx, a_mtx))
c_mtx = c_mtx.tocsr()
a_mtx = a_mtx.tocsr()
print type(a_mtx), type(c_mtx), a_mtx.shape, c_mtx.shape
%timeit uni_dot(a_mtx, c_mtx)
例1输出：
对照组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 360) (360, 30000)
1 loop, best of 3: 29.8 s per loop
实验组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 1440) (1440, 30000)
1 loop, best of 3: 28 s per loop
例2
a_mtx = sp.csr_matrix([[2.23, 1.56, 3.47]*120]*300)
mtx = sp.csr_matrix([[1.07, 2.19, 3.12]*120]*30000)
density(a_mtx)
density(mtx)
# 代码与例1的对应部分相同，不在重复
...
例2输出：
density 0.3333
density 0.3333
对照组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 360) (360, 30000)
1 loop, best of 3: 9.06 s per loop
实验组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 1440) (1440, 30000)
1 loop, best of 3: 8.85 s per loop
例3
a_mtx = sp.csr_matrix([[0., 0., 0., 0., 13.23, 0., 0., 0., 1.32, 0., 0., 0., 0., 0., 0., 0., 13.23, 0., 0., 0., 1.32, 0., 0., 0.]*5]*300)
mtx = sp.csr_matrix([[1.07, 0., 0., 0., 1.30, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.0, 0., 0., 0.]*5]*30000)
density(a_mtx)
density(mtx)
# 代码与例1的对应部分相同，不在重复
...
例2输出：
density 0.166666
density 0.083333
对照组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 120) (120, 30000)
1 loop, best of 3: 559 ms per loop
实验组：
<class 'scipy.sparse.csr.csr_matrix'> <class 'scipy.sparse.csr.csr_matrix'> (300, 480) (480, 30000)
1 loop, best of 3: 374 ms per loop
三 稀疏矩阵归一化和转置，不会影响矩阵计算性能
相同格式的稀疏矩阵做点乘速度很快，不同格式速度仅仅慢一丢丢。比如归一化和转置之后, 不转格式不会影响速度.
某些情况下在点乘计算前，需要进行归一化操作，比如计算cosine相似度，需要对两个稀疏矩阵分别做行归一化和列归一化，或者转置。在进行归一化或者转置后，矩阵的格式可能会发生改变.
这里使用的是sklearn.preprocessing.normalize函数进行归一化的。对于稀疏矩阵，行归一化的返回值是CSR矩阵，列归一化的返回值是CSC矩阵（实验结果见下面代码）；之所以这么这么做，是为了提高计算速度，同时也降低计算难度，sklearn的做法是：如果是sparse矩阵，当是行归一化时，就将原始矩阵转为CSR格式，这样就可以对矩阵的data(data是sparse.csr_matrix的一个属性)中的每行的元素，进行快速归一化。当列归一化时，转为CSC矩阵，然后对data中的列元素进行快速归一化。如果你不明白为什么如此操作的好处，请参看稀疏矩阵压缩原理。
转置操作输入CSR矩阵返回CSC矩，阵输入CSC矩阵返回CSR矩阵。至于转置为何也会改变矩阵格式，答案也是速度快，编码简单，为什么呢？自己动手计算一下吧。
