Быстрая накопленная сумма и мощность оператора - PullRequest
Новая
       19

Быстрая накопленная сумма и мощность оператора

0 голосов
/ 18 апреля 2019

У меня есть алгоритм прогноза, который обрабатывает тренд временных рядов до заданного горизонта, используя следующий код: 

import numpy as np
horizon = 91
phi = 0.2
trend = -0.004
trend_up_to_horizon = np.cumsum(phi ** np.arange(horizon) + 1) * self.trend

В этом примере первые два значения trend_up_horizon: 

array([-0.008 , -0.0128])

Есть ли в вычислительном отношении более быстрый способ добиться этого? На данный момент это занимает много времени, так как я полагаю, использование метода np.cumsum и оператора ** обходятся дорого.

Спасибо за любую помощь

Ответы [ 2 ]

2 голосов
/ 19 апреля 2019

Вы можете сделать эту операцию немного быстрее. Как вы уже предполагали, основной проблемой здесь является (ненужный) оператор питания.

В дополнение к этому у Numpy нет специальной реализации для power (float64, int64), где показатель степени представляет собой небольшое положительное целое число. Вместо этого Numpy всегда вычисляет мощность (float64, float64), что является гораздо более сложной задачей.

В Numba и Numexpr есть специальная реализация для простого случая (float64, int64), поэтому давайте попробуем это на первом шаге.

Первый подход 

import numpy as np
import numba as nb

horizon = 91
phi = 0.2
trend = -0.004

@nb.njit()
def pow_cumsum(horizon,phi,trend):
    out=np.empty(horizon)
    csum=0.
    for i in range(horizon):
        csum+=phi**i+1
        out[i]=csum*trend
    return out

Как уже говорилось, перед непосредственным вычислением мощности нет необходимости, алгоритм может быть переписан, чтобы избежать этого полностью.

Второй подход 

@nb.njit()
def pow_cumsum_2(horizon,phi,trend):
    out=np.empty(horizon)

    out[0]=2.*trend
    TMP=2.
    val=phi
    for i in range(horizon-1):
        TMP=(val+1+TMP)
        out[i+1]=TMP*trend
        val*=phi
    return out

Задержка 

%timeit np.cumsum(phi ** np.arange(horizon) + 1) * trend
7.44 µs ± 89.1 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit pow_cumsum(horizon,phi,trend)
882 ns ± 4.91 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
%timeit pow_cumsum_2(horizon,phi,trend)
559 ns ± 3.8 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)
1 голос
/ 18 апреля 2019

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

работает %timeit на вашем базовом np.cumsum(phi ** np.arange(horizon) + 1) * trend говорит, что на моем ноутбуке требуется 17,5 мкс, что немного

версия Cython, которая делает эквивалент: 

import numpy as np
cimport numpy as np
cimport cython

@cython.boundscheck(False)
def do_cumsum(size_t horizon, double phi, double trend):
    cdef np.ndarray[double, ndim=1] out = np.empty(horizon, dtype=np.float)
    cdef double csum = 0
    cdef int i

    for i in range(horizon):
        csum += phi ** i + 1
        out[i] = csum * trend

    return out

, что сокращает время do_cumsum(horizon, phi, trend) до 6,9 мкс, в то время как при переключении на 32-битные операции с одинарной точностью это уменьшается до 4,5 мкс

при этом микросекунды невелики, и вам, вероятно, лучше сосредоточить свои усилия в другом месте

...