Python: Speeding up vector and matrix operations
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