Advanced numpy array multiplication

We start with the 3 arrays of components x, y, and z. I will change the values from your example so that they have unique values:

x = np.array([[[1,2],[3,4]],
              [[5,6],[7,8]]])
y = x + 10
z = y + 10

Each of the above have shape (2,2,2), but they could be any (n, m, l). This shape will have little impact on our process.

We next combine the three component arrays into a new array p, the "position vector", creating a new dimension i will iterate over the three physical dimensions x, y, z,

p = np.array([x, y, z])

so p[0] is x and so on, and p has shape (d, n, m, l) (where d=3 is the physical dimensionality of the vectors).

Now we look at your list of vectors sv which presumably has shape (N, d). Let us use a small number for N:

N = 4
d = 3
sv = np.arange(d*N).reshape(N,d) # a list of N vectors in 3d

OK the above was a little repetive but I want to be clear (and please correct any misunderstandings I may have had from your question).

You want to make some difference, diff in which you take each of the n*m*l vectors buried in p and subtract from it each of the N vectors in sv. This will give you N*n*m*l vectors, which each have d components. We need to align each of these dimensions before we do subtractions.

Basically we want to take p - sv but we must make sure that their shapes match so that the d axis is aligned, and the n, m, l and N axes basically just add up. The way numpy broadcasts is to take the shapes of the array, and aligns them from the end, so the last axis of each is aligned, and so on. To broadcast, each size must match exactly, or must be empty (on the left) or 1. That is, if your shapes were (a, b, c) and (b, c), you would be fine, and the second array would be repeated ("broadcasted") a times to match the a different subarrays of shape (b, c) in the first array. You can use dimensions length 1 which will force the position, so normally two arrays of shape (a, b, c) and (a, b) will not align because the last axis does not match, but you can add a new placeholder axis at the end of the second to give it shape (a, b, 1) which will match to (a, b, c) no matter what the value of c is.

We give shape (N, d, 1, 1, 1) to sv which matches the shape (d, n, m, l) of p. This can be done several ways:

sv = sv.reshape(sv.shape + (1,1,1)])
#or
sv.shape += (1, 1, 1)
#or
sv = sv[..., None, None, None]

Then, we can do the difference:

diff = p - sv[..., None, None, None]

where we have that diff.shape is (N, d, n, m, l). Now we can square it and sum over the second (d) dimension to get the norm/magnitude of each vector:

m = (diff*diff).sum(1)

which of course will have shape (N, n, m, l), or in the example case (4, 2, 2, 2)

So, all together:

import numpy as np

x = np.array([[[1,2],[3,4]],
              [[5,6],[7,8]]])
y = x + 10
z = y + 10
p = np.array([x, y, z])
print p.shape
N = 4
d = 3
sv = np.arange(d*N).reshape(N,d) # a list of N vectors in 3d
print sv.shape
diff = p - sv[..., None, None, None]
print diff.shape
m = (diff*diff).sum(1)
print m.shape