title: Compute the pairwise distance in scipy with missing values
body: <p>I'm a bit stumped by how <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.pdist.html" rel="noreferrer"><code>scipy.spatial.distance.pdist</code></a> handles missing (<code>nan</code>) values.</p>
  
  <p>So just in case I messed up the dimensions of my matrix, let's get that out of the way. From the docs:</p>
  
  <blockquote>
    <p>The points are arranged as m n-dimensional row vectors in the matrix X.</p>
  </blockquote>
  
  <p>So let's generate three points in 10 dimensional space with missing values:</p>
  
  <pre><code>numpy.random.seed(123456789)
  data = numpy.random.rand(3, 10) * 5
  data[data &lt; 1.0] = numpy.nan
  </code></pre>
  
  <p>If I compute the Euclidean distance of these three observations:</p>
  
  <pre><code>pdist(data, "euclidean")
  </code></pre>
  
  <p>I get:</p>
  
  <pre>array([ nan,  nan,  nan])</pre>
  
  <p>However, if I filter all the columns with missing values I do get proper distance values:</p>
  
  <pre><code>valid = [i for (i, col) in enumerate(data.T) if ~numpy.isnan(col).any()]
  pdist(data[:, valid], "euclidean")
  </code></pre>
  
  <p>I get:</p>
  
  <pre>array([ 3.35518662,  2.35481185,  3.10323893])</pre>
  
  <p>This way, I throw away more data than I'd like since I don't need to filter the whole matrix but only the pairs of vectors being compared at a time. Can I make <code>pdist</code> or a similar function perform pairwise masking, somehow?</p>
  
  <hr>
  
  <p><strong>Edit:</strong></p>
  
  <p>Since my full matrix is rather large, I did some timing tests on the small data set provided here.</p>
  
  <p>1.) The scipy function.</p>
  
  <pre><code>%timeit pdist(data, "euclidean")
  </code></pre>
  
  <pre>10000 loops, best of 3: 24.4 µs per loop</pre>
  
  <p>2.) Unfortunately, the solution provided so far is roughly 10 times slower.</p>
  
  <pre><code>%timeit numpy.array([pdist(data[s][:, ~numpy.isnan(data[s]).any(axis=0)], "euclidean") for s in map(list, itertools.combinations(range(data.shape[0]), 2))]).ravel()
  </code></pre>
  
  <pre>1000 loops, best of 3: 231 µs per loop</pre>
  
  <p>3.) Then I did a test of "pure" Python and was pleasantly surprised:</p>
  
  <pre><code>from scipy.linalg import norm
  
  %%timeit
  m = data.shape[0]
  dm = numpy.zeros(m * (m - 1) // 2, dtype=float)
  mask = numpy.isfinite(data)
  k = 0
  for i in range(m - 1):
      for j in range(i + 1, m):
          curr = numpy.logical_and(mask[i], mask[j])
          u = data[i][curr]
          v = data[j][curr]
          dm[k] = norm(u - v)
          k += 1
  </code></pre>
  
  <pre>10000 loops, best of 3: 98.9 µs per loop</pre>
  
  <p>So I think the way to go forward is to Cythonize the above code in a function.</p>