Оптимизация умножения больших массивов с разреженными матрицами - PullRequest
Новая
       52

Оптимизация умножения больших массивов с разреженными матрицами

0 голосов
/ 03 мая 2018

Я делаю статистический расчет по следующей форме:

enter image description here

, где

  • d (данные) - матрица размера [i = 30, 100x100]
  • м (модель) - матрица размера [i = 30, 100x100]
  • C -1 (ковариация) представляет собой 100 2 на 100 2 симметричной матрицы

d и C -1 являются константами и не имеют никакой структуры, кроме C -1 быть симметричным. m изменяется, но всегда редкий . Результатом вычисления является просто число с плавающей точкой.

Мне нужно выполнить этот расчет много раз в симуляции Монте-Карло, поэтому скорость важна. Использование методов умножения разреженных массивов с m значительно ускоряет работу по сравнению с наивным матричным точечным произведением. Однако каждая итерация медленной функции, представленной ниже, по-прежнему занимает около 0,1 секунды. Подавляющее большинство времени (> 98%) тратится на умножения матриц, а не на функцию генерации моделей (generate_model). Я хотел бы ускорить это на порядок, если это возможно.

Код и выходные данные вставлены ниже.

Вещи, которые не работают, включают:

  • Обновление до подпрограмм линейной алгебры Intel MKL (ускорение на несколько процентов, на удивление небольшое)
  • Использование numpy.linalg.multi_dot
  • Воспользовавшись тем, что C -1 является симметричным (это не работает даже в принципе, см. этот вопрос о математическом потоке )

Вещи, которые работают , включают:

  • Предварительный расчет 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
...