Vectorise Python code
<p>I have coded a kriging algorithm but I find it quite slow. Especially, do you have an idea on how I could vectorise the piece of code in the cons function below:</p>

<pre><code>import time
import numpy as np

B = np.zeros((200, 6))
P = np.zeros((len(B), len(B)))

def cons():
  time1=time.time()
  for i in range(len(B)):
    for j in range(len(B)):
      P[i,j] = corr(B[i], B[j])
  time2=time.time()
  return time2-time1

def corr(x,x_i):
  return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))))    

time_av = 0.
for i in range(30):
  time_av+=cons()
print "Average=", time_av/100.
</code></pre>

<hr>

<p>Edit: Bonus questions</p>

<ol>
<li>What happens to the broadcasting solution if I want <code>corr(B[i], C[j])</code> with C the same dimension than B</li>
<li><p>What happens to the scipy solution if my p-norm orders are an array:</p>

<pre><code>p=np.array([1.,2.,1.,2.,1.,2.])
def corr(x, x_i):
  return np.exp(-np.sum(np.abs(np.array(x) - np.array(x_i))**p))  
</code></pre>

<p>For 2., I tried <code>P = np.exp(-cdist(B, C,'minkowski', p))</code> but scipy is expecting a scalar.</p></li>
</ol>

<p>Your problem seems very simple to vectorize. For each pair of rows of <code>B</code> you want to compute</p>

<pre><code>P[i,j] = np.exp(-np.sum(np.abs(B[i,:] - B[j,:])))
</code></pre>

<p>You can make use of array broadcasting and introduce a third dimension, summing along the last one:</p>

<pre><code>P2 = np.exp(-np.sum(np.abs(B[:,None,:] - B),axis=-1))
</code></pre>

<p>The idea is to reshape the first occurence of <code>B</code> to shape <code>(N,1,M)</code> while the second <code>B</code> is left with shape <code>(N,M)</code>. With array broadcasting, the latter is equivalent to <code>(1,N,M)</code>, so</p>

<pre><code>B[:,None,:] - B
</code></pre>

<p>is of shape <code>(N,N,M)</code>. Summing along the last index will then result in the <code>(N,N)</code>-shape correlation array you're looking for.</p>

<hr>

<p>Note that if you were using <code>scipy</code>, you would be able to do this using <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html" rel="nofollow"><code>scipy.spatial.distance.cdist</code></a> (or, equivalently, a combination of <code>scipy.spatial.distance.pdist</code> and <code>scipy.spatial.distance.squareform</code>), without unnecessarily computing the lower triangular half of this symmetrix matrix. Using <a href="http://stackoverflow.com/users/3293881/divakar">@Divakar</a>'s suggestion in comments for the simplest solution this way:</p>

<pre><code>from scipy.spatial.distance import cdist
P3 = 1/np.exp(cdist(B, B, 'minkowski',1))
</code></pre>

<p><code>cdist</code> will compute the Minkowski distance in 1-norm, which is exactly the sum of the absolute values of coordinate differences.</p>

python numpy vectorization