python - 使用稀疏矩阵优化大型 numpy 数组乘法

Question

我正在做以下形式的统计计算：

在哪里

• d (data) 是大小为 [i=30, 100x100] 的矩阵
• m (model) 是大小为 [i=30, 100x100] 的矩阵
• C ^-1（协方差）是 100 ² x 100 ²对称矩阵

d和C ^-1是常数，除了C ^-1是对称的之外没有任何结构。m改变，但总是稀疏的。计算的输出只是一个浮点数。

我需要在蒙特卡洛模拟中多次执行此计算，因此速度至关重要。使用带有m的稀疏数组乘法技术可以大大加快简单矩阵点积的速度。但是，下面的慢函数的每次迭代仍然需要大约 0.1 秒才能运行。绝大多数时间 (>98%) 都花在了矩阵乘法上，而不是模型生成函数 ( generate_model) 上。如果可能的话，我想把它加快一个数量级。

代码和输出粘贴在下面。

不起作用的事情包括：

• 升级到英特尔 MKL 线性代数例程（加速百分之几，小得惊人）
• 使用numpy.linalg.multi_dot
• 利用C ^-1是对称的这一事实（即使在原则上这也不起作用，请参阅这个 mathoverflow 问题）

这类工作包括：

• 预先计算C ^-1 d，加速约 40%

如何加快此代码的速度？ 依赖于 , 等软件包的解决方案cython以及numba“标准” scipy/numpy 解决方案是最受欢迎的。提前致谢！

from __future__ import division
import numpy as np
import scipy.sparse
import sys 
import timeit

def generate_model(n, size, hw = 8): 
    #model for the data--squares at random locations
    output = np.zeros((n, size, size))
    for i in range(n):
        randx = np.random.randint(hw, size-hw)
        randy = np.random.randint(hw, size-hw)
        output[i,(randx-hw):(randx+hw), (randy-hw):(randy+hw)]=np.random.random((hw*2, hw*2))
    return output

def slow_function(datacube, invcovmatrix, size):
    model = generate_model(30, size) 
    output = 0 
    for i in range(model.shape[0]):
        data = datacube[i,:,:].flatten()
        mu = model[i,:,:].flatten()
        sparsemu = scipy.sparse.csr_matrix(mu)
        output += -0.5* (
                  np.float(-2.0*sparsemu.dot(invcovmatrix).dot(data)) +
                  np.float(sparsemu.dot(sparsemu.dot(invcovmatrix).T))
                  )   
    return output

def wrapper(func, *args, **kwargs):
    def wrapped():
        return func(*args, **kwargs)
    return wrapped

if __name__ == "__main__":
    size = 100 
    invcovmat    = np.random.random((size**2, size**2))
    #make symmetric for consistency
    invcovmat    = (invcovmat+invcovmat.T)/2
    datacube     = np.random.random((30, size, size))
    #TIMING
    wrapped = wrapper(slow_function, datacube, invcovmat, size)
    times = []
    for i in range(20):
        print i
        times.append(timeit.timeit(wrapped, number = 1)) 
    times.sort()
    print '\n', np.mean(times[0:3]), ' s/iteration; best of 3'

输出：

0.10408163070678711  s/iteration; best of 3

标签： pythonnumpymatrixsparse-matrix