Let me first discuss how to construct the matrix of distances between a set of vectors
\{v_i\}
. The idea is, obviously, to use the fact that the p-distance between two vectors is given by the formula d_p(v_1,v_2) = \|v_1-v_2\|_p = \left(\sum_k (v^k_1-v^k_2)^p\right)^{1/p}
For the Euclidean distance $p=2$ and we have
\|v\|^2 = (v,v)
. So the squared distance is nothing but \|v_1-v_2\|^2 = (v_1-v_2,v_1-v_2)
This is good, then using the linearity of the scalar product, we obtain
\|v_1-v_2\|^2 = \|v_1\|^2 + \|v_2\|^2 -2(v_1,v_2)
This expression can be computed with matrix multiplications. In python you can do it using numpy as follows. First, onstruct the matrix of the positions, i.e. stack all 'size' vectors of lenght 'dimension' on the top of each other import numpyHere I have chosen uniformly distributed vectors, but you can use others of course.
dimension = 2
size = 100
positions = numpy.random.uniform(0,1, (size,dimension))
Now, we construct the matrix
s_{ij} = \|v_i\|^2 +\|v_j\|^2
by repeating, reshaping and transposing the vector of the norms. This is as easy as this# construct the matrix s_ij = |v_i|**2+|v_j|**2'sum_matrix' is what you are looking for. The scalar product is even easier. Indeed the matrix of the products
norms = numpy.sum( positions**2. , axis = 1 )
tmp = numpy.reshape(norms.repeat(size),(size,size))
sum_matrix = tmp + tmp.transpose()
x_{ij} = (v_i,v_j)
is just the multiplication of the vector matrix with its transpose (try on 2x2 example to see that it works). So you can do it easily by# construct the matrix x_ij = (v_i,v_j)
scalars = numpy.dot(positions,positions.transpose())