Using numpy for the summation where i != j

Viewed 147

I need to calculate the sum below where X is a matrix 2x5 and i/j is a selected column.

sum

I've solved this with the following code, but I would like to know if there is an easier solution with numpy to solve this sum (without the for..range)

sum = 0
size = 5
for i in range(0, size):
    for j in range(0, size):
        if i != j:
            sum = sum + np.power(np.linalg.norm(X[:, i] - X[:, j]), 2)
2 Answers

You can optimize quite a bit here. First, you only need to compute for off-diagonal elements of the full i-j matrix. Second, you can compute for just one half of the matrix and double the sum, since norm(a - b) == norm(b - a).

Use np.triu_indices (or np.tril_indices):

i, j = np.triu_indices(size, 1)
sum = 2 * (np.linalg.norm(X[:, i] - X[:, j], axis=0)**2).sum()

You can avoid the square roots (assuming you're using 2-norm):

sum = 2 * (X[:, i] - X[:, j])**2).sum(None)

(np.array([np.subtract.outer(X,X)[0,:,0,:],np.subtract.outer(X,X)[1,:,1,:]])**2).sum()

seems to work .... I'd definitely stick with the for loops (or if this were time-critical production code, I'd probably use Cython ... This "solution" does involve more computations than it needs to).

However, removing one for loop can be done quite a bit more cleanly:

sum = 0
size = 5
XT = X.T
for i in range(0, size):
    sum = sum + ((XT[i, :] - XT[:, :])**2).sum()

... Also it's not a good idea to use sum as a variable name. It overrides the built-in sum function!

Related