Is it possible to make these matrix operations without loops in python 2.7 (numpy)?

Question

I want to implement the loops below using matrices in Python:

import numpy as np

n = 5 # samples
k = 2 # inputs
m = 3 # gaussians

# X is nxk
X = np.array([[0.0, 10.0], [20.0, 30.0],[40, 50],[60,70],[80,90]])

#locations is mxk
locations = np.array([[0.01, 0.02], [0.03,0.04], [0.05, 0.06]])

dev = np.empty([n,k,m])

for samples in range(n):
    for inputs in range(k):
        for gaussians in range(m):
            dev[samples,inputs,gaussians]=X[samples,inputs]-locations[gaussians,inputs]

output = np.empty([n,m]) 

for samples in range(n):
    for gaussians in range(m):
        output[samples,gaussians]=np.sum(dev[samples,:,gaussians]*dev[samples,:,gaussians])

I know that Numpy is able to make operations using arrays of different dimensions (Broadcast), but I am not able to use this concept here. Note that I am basically doing is removing the average of a sample of vectors and calculating its squared norm.

Answer 1

You could vectorize your for loops like this; The dev is basically an outer operation of X and locations with respect to the first dimension, so you can insert a new axis in locations (or X, this will only affect how you will transpose the result) and the subtraction will trigger numpy broadcasting and give back a cartesian/outer subtraction; For the second, you need to multiply dev with itself and sum along the second dimension (axis = 1):

mydev = np.transpose(X - locations[:,None], (1,2,0))

(mydev == dev).all()
# True

myoutput = (mydev**2).sum(axis=1)

(myoutput == output).all()
# True

---

Or to put together:

((X[:,None] - locations) ** 2).sum(axis=-1)

#array([[    99.6005,     99.2025,     98.8061],
#       [  1298.4005,   1296.4025,   1294.4061],
#       [  4097.2005,   4093.6025,   4090.0061],
#       [  8496.0005,   8490.8025,   8485.6061],
#       [ 14494.8005,  14488.0025,  14481.2061]])