Internal vector arrangement and numpy automatic broadcasting

1

I wanted to interpolate an expression of the form:

f\left(\mathbf{x}\right) = \left\Vert \mathbf{x} - \mathbf{a} \right\Vert

which meant evaluating

\left\Vert\mathbf{x}_i-\mathbf{a}\right\Vert,\quad i=1,2,\dots,N_\mathrm{points}

Me being used to numpy automatic broadcasting, wrote:

f.interpolate( lambda x : np.linalg.norm(x - a))

where a is a 3 elements array.

However that does not work, as x is a 3 \times N_\mathrm{points} array, and numpy does automatic broadcasting from the right.

I can easily fix it by calling x.transpose() inside the expression, however the code is less clean. (my actual expression is a little bit more cumbersome)

Is there a reason for choosing the d \times N_\mathrm{points} arrangement instead of the transposed one?

2

You can give numpy.linalg.norm the axis to evaluate over: numpy.linalg.norm — NumPy v2.1 Manual

See the following example in the docs:

c = np.array([[ 1, 2, 3],
              [-1, 1, 4]])
LA.norm(c, axis=0)
array([ 1.41421356,  2.23606798,  5.        ])
LA.norm(c, axis=1)
array([ 3.74165739,  4.24264069])
LA.norm(c, ord=1, axis=1)
array([ 6.,  6.])

The reason for choosing this is that it allows for vectorized operations in Python, rather than looping over each point, which would be inefficient. (It also gives nice memory alignment).

3

Following up, the reason for shape=(3, npoints) is that it allows, for example, f = y \cos(x) to be interpolated using

f.interpolate(lambda x : x[1] * np.cos(x[0]))
1 Like
4

nice! dokken answer was already convincing, but this one is also beautiful