具有常量的 Python 矢量化

Question

我有一个长度为 n(=300,000) 的序列 X。使用 w (=40) 的窗口长度，我需要实现:

mu(i)= X(i)-X(i-w)

s(i) = sum{k=i-w to i} [X(k)-X(k-1) - mu(i)]^2

我想知道这里是否有防止循环的方法。 mu(i) 在第二个方程中是常数这一事实导致了向量化的复杂化。到目前为止我做了以下事情:

x1=x.shift(1)
xw=x.shift(w)
mu= x-xw 
dx=(x-x1-mu)**2 # wrong because mu wouldn't be constant for each i
s=pd.rolling_sum(dx,w)

上面的代码可以在循环设置中工作(并且正在工作)但花费的时间太长，因此任何有关矢量化或其他速度改进方法的帮助都会有所帮助。我用 mathjax 格式在 crossvalidated 上发布了这个，但这似乎在这里不起作用。

https://stats.stackexchange.com/questions/241050/python-vectorization-with-a-constant

还要澄清一下，我最初没有使用双循环，只是一个循环:

        for i in np.arange(w, len(X)):
            x=X.ix[i-w:i,0] # clip a series of size w
            x1=x.shift(1)   
            mu.ix[i]= x.ix[-1]-x.ix[0]   
            temp= (x-x1-mu.ix[i])**2   # returns a series of size w but now mu is constant
            s.ix[i]= temp.sum()

Answer 1

方法 #1: 一种矢量化方法是使用 broadcasting -

N = X.shape[0]
a = np.arange(N)
k2D = a[:,None] - np.arange(w+1)[::-1]
mu1D = X - X[a-w]
out = ((X[k2D] - X[k2D-1] - mu1D[:,None])**2).sum(-1)

我们可以进一步优化最后一步以获得平方和 np.einsum -

subs = X[k2D] - X[k2D-1] - mu1D[:,None]
out = np.einsum('ij,ij->i',subs,subs)

使用 NumPy strides 可以进一步改进得到 X[k2D] 和 X[k2D-1]。

---

方法 #2: 为了在处理非常大的数组时节省内存，我们可以使用一个循环而不是原始代码中使用的两个循环，就像这样 -

N = X.shape[0]
s = np.zeros((N))
k_idx = np.arange(-w,1)
for i in range(N):
    mu = X[i]-X[i-w]
    s[i] = ((X[k_idx]-X[k_idx-1] - mu)**2).sum()
    k_idx += 1

同样，np.einsum 可以在此处用于计算 s[i]，就像这样 -

subs = X[k_idx]-X[k_idx-1] - mu
s[i] = np.einsum('i,i->',subs,subs)