Python: Speeding up vector and matrix operations The core of my program are the following lines

import numpy as np
import time

I, J = np.random.randn(20000,800), np.random.randn(10000,800)

d = np.zeros((len(I),1))
for j in range(len(J)):
    t0 = time.time()
    # Compute distance
    matrix_j = np.array([J[j]]).T*np.array([J[j]])
    for i in range(len(I)): 
        matrix_i = np.array([I[i]]).T*np.array([I[i]])
        d[i] = np.sum((matrix_i-matrix_j)**2)
    print(time.time()-t0)
    # Do stuff with vector d


I and J are some matrices with data. I basically compare vectors of each matrix where I first extract a vector and than multiply it with its transposed version to get a matrix matrix_i and matrix_j. Afterwards I subtract both matrices from each other and sum all elements. But it takes ages even for such a small data set.

What can I do to speed things a little bit up? Any ideas?
What you compute with the presented code is 

matrix_j = array([ J[j,k]**2 for k in range(N) ]) 


(where N=800 is the vector length) so that 

d[i] = sum( (J[j,k]**2 - I[i,k]**2)**2 for k in range(N) )


This does not make much sense geometrically. 



What you could have intended, by what might be implied by the formulas and variable names, is to have J[j] and I[i] as actual row vectors which can be achieved by transforming the tables I,J into matrices,

I,J = np.matrix(I), np.matrix(J)


so that then d[i] is the square of the Frobenius norm of J[j].T*J[j]-I[i].T*I[i], which is also the sum over the squares over all matrix elements. As this matrix is symmetric, this quantity can be computed as 

d[i] = np.sum(np.diag( (J[j].T*J[j]-I[i].T*I[i])**2 ))


which has the form for row vectors \$a,b\$ of equal dimension
\begin{align}\|a^Ta-b^Tb\|_F^2 &amp;= \text{trace}((a^Ta-b^Tb)^2)\\ &amp;=\text{trace}((a^Ta-b^Tb)(a^Ta-b^Tb))\\ &amp;=\text{trace}(a^Taa^Ta-a^Tab^Tb-b^Tba^Ta+b^Tbb^Tb))\\ &amp;=(aa^T)^2-2(ab^T)^2+(bb^T)^2\\ &amp;=\|a\|^4-2\langle a,b\rangle ^2+\|b\|^4 \end{align}

Re-expressed in the rows of the original matrices, that should be 2D-arrays, not converted to matrix objects, this is

d[i] = np.dot(I[i],I[i])**2 - 2*np.dot(I[i],J[j])**2 + np.dot(J[j],J[j])**2


For a component-wise derivation see the answer of Gareth Rees.



These norm squares and scalar products should be faster to compute than the original formula. Also, pre-computing the norm squares np.dot(I[i],I[i]) allows re-use so that all len(I)+len(J) norm squares are only computed once and the main computational cost is from the len(I)*len(J) scalar products.

scipy.linalg.norm might be faster in computing the norm, as squaring can be faster than multiplying the same number to itself, but I'm not sure how that translates over the several layers of interpretation and data encapsulation.



The scalar products are elements of a matrix product of I and J.T, so that a compact computation proceeds by (all are np.array objects)

dotIJ = np.dot(I,J.T); 
normI2 = [ np.dot(row,row) for row in I ]; 
normJ2 = [ np.dot(row,row) for row in J ]; 
d = np.zeros((len(I),1)); 
for j in range(len(J)):
    d[:,0] = [ normI2[i]**2 + normJ2[j]**2 - 2*dotIJ[i,j]**2 for i in range(len(I)) ]
    # process values of d