Python: Speeding up vector and matrix operations

Question

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?

Answer 1

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 &= \text{trace}((a^Ta-b^Tb)^2)\\ &=\text{trace}((a^Ta-b^Tb)(a^Ta-b^Tb))\\ &=\text{trace}(a^Taa^Ta-a^Tab^Tb-b^Tba^Ta+b^Tbb^Tb))\\ &=(aa^T)^2-2(ab^T)^2+(bb^T)^2\\ &=\|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

Answer 2

Here's an explicit working-out of the maths needed for the answer by LutzL.

Let's start by writing \$I\$ and \$J\$ for the two vectors I[i] and J[j] respectively and \$n\$ for the length of these vectors (800 in the example in the post). Then the distance \$d\$ that we want to compute is $$\sum_k \sum_l \left(I_k I_l − J_k J_l\right)^2$$ where \$k\$ and \$l\$ range over the indices from \$0\$ to \$n-1\$. Expand the square: $$\sum_k \sum_l \left(I_k^2 I_l^2 + J_k^2 J_l^2 − 2 I_k I_l J_k J_l\right).$$ Extract the sub-expressions that don't depend on \$l\$: $$\sum_k \left(I_k^2 \sum_l I_l^2 + J_k^2 \sum_l J_l^2 − 2 I_k J_k \sum_l I_l J_l\right).$$ Multiply out: $$ \left(\sum_k I_k^2\right) \left(\sum_l I_l^2\right) + \left(\sum_k J_k^2\right) \left(\sum_l J_l^2\right) − 2 \left(\sum_k I_k J_k\right) \left(\sum_l I_l J_l\right).$$ Now there are no nested summations, so rename \$l\$ to \$k\$ and then combine identical sums: $$ \left(\sum_k I_k^2\right)^2 + \left(\sum_k J_k^2\right)^2 - 2 \left(\sum_k I_k J_k\right)^2.$$ Finally, rewrite the sums as dot products: $$\left(I·I\right)^2 + \left(J·J\right)^2 - 2\left(I·J\right)^2$$