Быстрая итерация по массиву для квадратов остатков

Question

Мне нравится, когда данные методом наименьших квадратов (массив с плавающей точкой) совпадают со многими известными формами сигналов. Мой код работает, но он слишком медленный для многих запусков, которые я планирую сделать:

import numpy
import time

samples = 50000
width_signal = 100
data = numpy.random.normal(0, 1, samples)
signal = numpy.random.normal(0, 1, width_signal)  # Placeholder

t0 = time.clock()
for i in range(samples - width_signal):
    data_chunk = data[i:i + width_signal]
    residuals = data_chunk - signal
    squared_residuals = residuals**2
    summed_residuals = numpy.sum(squared_residuals)
t1 = time.clock()
print('Time elapsed (sec)', t1-t0)

РЕДАКТИРОВАТЬ: Исправлена ​​ошибка: сначала квадрат остатков, а затем их суммировать.

Это займет около 0,2 секунды для запуска на моей машине. Поскольку у меня много наборов данных и форм сигналов, это слишком медленно. Моя конкретная проблема не учитывает типичные методы MCMC, потому что формы сигналов слишком разные. Это должна быть грубая сила.

Типичные объемы - 50000 поплавков для данных и 100 для сигнала. Они могут варьироваться в несколько раз.

Мои тесты показывают, что:

• Суммирование данных numpy.sum(residuals) кушает 90% времени. Я попробовал Python sum(residuals), и он быстрее для маленьких массивов (~ <50 элементов) и медленнее для больших массивов. Должен ли я вставить условие <code>if?
• Я попытался numpy.roll() вместо прямой выборки данных, а .roll() медленнее.

Вопросы:

• Есть ли логическое улучшение для ускорения?
• Есть ли более быстрый способ суммирования массивов? Я не знаю C, но если бы это было намного быстрее, я мог бы попробовать это.
• Может ли помочь графический процессор? У меня много пробежек. Если да, где я могу найти фрагмент кода для этого?

Answer 1

Помимо версий, предоставленных Divakar, вы можете просто использовать компилятор, такой как Numba или Cython.

Exmaple

import numba as nb
@nb.njit(fastmath=True,parallel=True)
def sq_residuals(data,signal):
  summed_residuals=np.empty(data.shape[0]+1-signal.shape[0],dtype=data.dtype)
  for i in nb.prange(data.shape[0] - signal.shape[0]+1):
      sum=0.
      for j in range(signal.shape[0]):
        sum+=(data[i+j]-signal[j])**2
      summed_residuals[i]=sum
  return summed_residuals

Задержка

Numba 0.4dev, Python 3.6, Quadcore i7
MSD_uniffilt_conv(Divakar): 2.4ms

after the first call which invokes some compilation overhead:
sq_residuals              : 1.7ms

Answer 2

Основываясь на различных методах, предложенных в Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy, мы рассмотрим наш случай здесь.

Подход № 1

Мыможет использовать np.lib.stride_tricks.as_strided на основе scikit-image's view_as_windows, чтобы получить раздвижные окна и, таким образом, получить наше первое решение здесь, вот так -

from skimage.util import view_as_windows

d = view_as_windows(data,(width_signal))-signal # diffs
out = np.einsum('ij,ij->i',d,d)

Подробнеена использовании as_strided на основе view_as_windows.

подход № 2

Опять же, основываясь на приеме умножения матриц в этом ответном посте, мы могли быулучшить производительность, вот так -

def MSD_strided(data, signal):
    w = view_as_windows(data,(width_signal))
    return (w**2).sum(1) + (signal**2).sum(0) - 2*w.dot(signal)

Подход № 3

Мы будем совершенствовать подход № 2 путем обеспечения равномерной фильтрации и свертки -

from scipy.ndimage.filters import uniform_filter 

def MSD_uniffilt_conv(data, signal):
    hW = width_signal//2
    l = len(data)-len(signal)+1
    parte1 = uniform_filter(data**2,width_signal)[hW:hW+l]*width_signal
    parte3 = np.convolve(data, signal[::-1],'valid')    
    return parte1 + (signal**2).sum(0) - 2*parte3

Бенчмаркинг

Времена на вывешенном образце -

In [117]: %%timeit
     ...: for i in range(samples - width_signal + 1):
     ...:     data_chunk = data[i:i + width_signal]
     ...:     residuals = data_chunk - signal
     ...:     squared_residuals = residuals**2
     ...:     summed_residuals = numpy.sum(squared_residuals)
1 loop, best of 3: 239 ms per loop

In [118]: %%timeit
     ...: d = view_as_windows(data,(width_signal))-signal
     ...: np.einsum('ij,ij->i',d,d)
100 loops, best of 3: 11.1 ms per loop

In [209]: %timeit MSD_strided(data, signal)
10 loops, best of 3: 18.4 ms per loop

In [210]: %timeit MSD_uniffilt_conv(data, signal)
1000 loops, best of 3: 1.71 ms per loop

~140x Ускорение с третьим!