# Pairwise distance and residual calculation

I have a code for calculating the pairwise distances and the residuals of my data (X, Y, Z). Data is quite large (average of 7000 rows) and so my interest is code efficiency. My initial code is

import tkinter as tk
from tkinter import filedialog
import pandas as pd
import, numpy as np
from scipy.spatial.distance import pdist, squareform

root = tk.Tk()
root.withdraw()

data = np.array(data, dtype=np.float)
npoints, cols = data.shape

pwdistance = np.zeros((npoints, npoints))
pwresidual = np.zeros((npoints, npoints))
for i in range(npoints):
for j in range(npoints):
pwdistance[i][j] = np.sqrt((data[:,0][i]-data[:,0][j])**2 + (data[:,1][i]-data[:,1][j])**2)
pwresidual[i][j] = (data[:,2][i]-data[:,2][j])**2


With the pwdistance, I changed it to the following below which works extremely good.

pwdistance = squareform(pdist(data[:,:2]))


Is there a pythonic way of calculating my pwresidual, so I do not need to use a loop and to enable my code to run faster?

• The comma in import, numpy as np is a typo, right? – Stop ongoing harm to Monica Jan 18 '17 at 20:14
• I had the impression that CR was more strict than SO about the completeness of the code. This example is not self contained, and cannot be tested without constructing data, and correcting at least one typo. – hpaulj Jan 19 '17 at 20:27

Let's start by setting

col2 = data[:,2]


to avoid some repetition. (If I knew what was in the data array then maybe I could pick a better name, but all I have to go on is what was in the post.)

Now here are two possible approaches to computing pwresidual:

1. pwresidual[i,j] is the result of an operation on col2[i] and col2[j]. You might recognize this as being similar to an outer product in mathematics — except that the operation is subtraction rather than multiplication.

But in NumPy every universal function has an "outer" counterpart, and np.subtract is a universal function, so all we need is:

pwresidual = np.subtract.outer(col2, col2) ** 2

2. Alternatively, we can use NumPy's broadcasting mechanism, combined with np.newaxis to ensure every pair of items gets operated on:

pwresidual = (col2[:,np.newaxis] - col2[np.newaxis,:]) ** 2


Usually in numpy we write

pwresidual[i][j] = (data[:,2][i]-data[:,2][j])**2


as

pwresidual[i,j] = (data[i,2] - data[j,2])**2


Here's it's a style issue more than a functional one, but sometimes the distinction is important.

But as Gareth shows you don't need to iterate - just broadcast

pwresidual = (data[:,None,2]-data[None,:,2])**2


It subtracts a row array from a column array:

In [95]: data=np.arange(12).reshape(4,3)
In [96]: data
Out[96]:
array([[ 0,  1,  2],
[ 3,  4,  5],
[ 6,  7,  8],
[ 9, 10, 11]])
In [97]: data[:,None,2]
Out[97]:
array([[ 2],
[ 5],
[ 8],
[11]])
In [98]: data[None,:,2]
Out[98]: array([[ 2,  5,  8, 11]])
In [99]: data[:,None,2]-data[None,:,2]
Out[99]:
array([[ 0, -3, -6, -9],
[ 3,  0, -3, -6],
[ 6,  3,  0, -3],
[ 9,  6,  3,  0]])


as for the other expression:

pwdistance[i][j] = np.sqrt((data[:,0][i]-data[:,0][j])**2 + (data[:,1][i]-data[:,1][j])**2)

pwdistance[i, j] = np.sqrt((data[i,0]-data[j,0])**2 + (data[i,1]-data[j,1])**2)

pwdistance = np.sqrt( (data[:,None,0]-data[None,:,0])**2 + (data[:,None,1]-data[None,:,1])**2)


Together we can do

temp = (data[:,None,:] - data[None,:,:])**2
pwdistance = np.sqrt( temp[:,:,0] + temp[:,:,1])
pwresidual = temp[:,:,2]


WIth the above test data:

In [101]: pwdistance
Out[101]:
array([[  0.        ,   4.24264069,   8.48528137,  12.72792206],
[  4.24264069,   0.        ,   4.24264069,   8.48528137],
[  8.48528137,   4.24264069,   0.        ,   4.24264069],
[ 12.72792206,   8.48528137,   4.24264069,   0.        ]])
In [102]: pwresidual
Out[102]:
array([[ 0,  9, 36, 81],
[ 9,  0,  9, 36],
[36,  9,  0,  9],
[81, 36,  9,  0]])

• I think that using np.newaxis (instead of None) in index expressions makes it clearer what the intention is. – Gareth Rees Jan 23 '17 at 10:47