Vectorized Evaluation of a Function, Broadcasting and Element Wise Operations

Question

Given this...

I have to explain what this code does, knowing that it performs the vectorized evaluation of F, using broadcasting and element wise operations concepts...

def F(x_pos, alpha):
    D = x_pos.reshape(1,-1) - x_pos.reshape(-1,1)
    return (1./alpha) * (alpha.reshape(1,-1) * R(D)).sum(axis=1)

My explanation is:

In the first line of the function F receives x_pos and alpha as parameters (both numpy arrays), in the second line the matrix D is calculated by means of broadcasting (basic operations such as addition in arrays numpy are performed elementwise, ie, element by element, but it is also possible with arranys of different size if numpy can transform them into others of the same size, this conversion is called broadcasting), subtracting an array of order 1xN with another of order Nx1, resulting in the matrix D of order NxN containing x_j - x_1, x_j - x_2, etc. as elements, finally, in the last line the reciprocal of alpha is calculated (which clearly is an arrangement), where each element is multiplied by the sum of the R evaluation of each cell of the matrix D multiplied by alpha_j horizontally (due to axis = 1 in the argument)

Questions:

1. Considering I'm new to Python, is my explanation OK?
2. The code has an error or not? Because I don't see that the "j must be different from 1, 2, ..., n" in each sum is taken into consideration in the code... and If it's in fact wrong... How can I fix the code so it do exactly the same thing as stated as in the image?

Answer 1

Few comments/improvements/fixes could be suggested here.

1] The first step could be alternatively done with just introducing a new axis and subtracting with itself, like so -

D = x_pos[:,None] - x_pos

In my opinion, this is a cleaner option. The performance benefit might be just marginal.

2] In the second line, I think it needs a fix as we need to avoid computations for the diagonal elements of R(D). So, If I got that correctly, the corrected code would be -

vals = R(D)
np.fill_diagonal(vals,0)
out = (1./alpha) * (alpha.reshape(1,-1) * vals).sum(axis=1)

Now, let's make the code a bit more idiomatic/cleaner.

At that line, we could write : (alpha * vals) instead of alpha.reshape(1,-1) * vals. This is because the shapes are already aligned for broadcasting as shown in a schematic diagram below -

alpha :      n
vals  :  n x n

Thus, alpha would be automatically extended to 2D with its elements broadcasted along the first axis for the length of vals and then elementwise multiplications being generated with it. Again, this is meant as a cleaner code.

There's a further performance improvement possible here with (alpha.reshape(1,-1) * vals).sum(axis=1) being replaceable with a matrix-multiplicatiion using np.dot as alpha.dot(vals). The benefit on performance should be noticeable with this step.

So, the second step reduces to -

out = (1./alpha) * alpha.dot(vals)