Python vectorization with a constant
<p>I have a series X of length n(=300,000). Using a window length of w (=40), I need to implement:</p>

<p>mu(i)= X(i)-X(i-w)</p>

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

<p>I was wondering if there's a way to prevent loops here. The fact that mu(i) is constant in second equation is causing complications in vectorization. I did the following so far:</p>

<pre><code>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)
</code></pre>

<p>The above code would work (and was working) in a loop setting but takes too long, so any help regarding vectorization or other speed improvement methods would be helpful. I posted this on crossvalidated with mathjax formatting but that doesn't seem to work here.</p>

<p><a href="http://stats.stackexchange.com/questions/241050/python-vectorization-with-a-constant">http://stats.stackexchange.com/questions/241050/python-vectorization-with-a-constant</a></p>

<p>Also just to clarify, I wasn't using a double loop, just a single one originally:</p>

<pre><code>        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()
</code></pre>

<p><strong>Approach #1 :</strong> One vectorized approach would be using <a href="https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html" rel="nofollow"><code>broadcasting</code></a> -</p>

<pre><code>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)
</code></pre>

<p>We can further optimize the last step to get squared summations with <a href="https://docs.scipy.org/doc/numpy/reference/generated/numpy.einsum.html" rel="nofollow"><code>np.einsum</code></a> -</p>

<pre><code>subs = X[k2D] - X[k2D-1] - mu1D[:,None]
out = np.einsum('ij,ij-&gt;i',subs,subs)
</code></pre>

<p>Further improvement is possible with the use of <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.strides.html" rel="nofollow"><code>NumPy strides</code></a> to get <code>X[k2D]</code> and <code>X[k2D-1]</code>.</p>

<hr>

<p><strong>Approach #2 :</strong> To save on memory when working very large arrays, we can use one loop instead of two loops used in the original code, like so -</p>

<pre><code>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
</code></pre>

<p>Again, <code>np.einsum</code> could be used here to compute <code>s[i]</code>, like so -</p>

<pre><code>subs = X[k_idx]-X[k_idx-1] - mu
s[i] = np.einsum('i,i-&gt;',subs,subs)
</code></pre>

python performance pandas numpy vectorization