Broadcasting - 3D field of coefficients to 3D field of matrices given matrix basis

0
votes

I have a (large) 4D array, consisting of the 5 coefficients in a given basis for a matrix field. Given the 5 basis matrices, I want to efficiently calculate the matrix field.

The coefficient field c[x,y,z,i] being the value of i-th coefficient at position x,y,z

And the matrix field M[x,y,z,a,b] being the (3,3) matrix at position x,y,z

And the basis matrices T_1,...T_5, being the (3,3) basis matrices

I could loop over each position in space:

M[x,y,z,:,:] = T_1[:,:]*c[x,y,z,0] + T_2[:,:]*c[x,y,z,1]...T_5[:,:]*c[x,y,z,4]

But this is very inefficient. My attempts at using np.multiply,np.sum result in broadcasting errors due to the ambiguity of the desired product being a field of 3x3 matrices.

1
pythonnumpymatrix
So in your text you mention that M[x,y,z,a,b] and then later when you calculate it M[x,y,z] Which one is the dimension? Can you provide a dummy example? – anishtain4
Quick comment, it appears that np.apply_along_axis() is a solution when supplied with a coefficient-->matrix conversion function. – Smectic
That is a very inefficient way to attach your problem, since you are not having any laborious function that cannot be done with tensors. look at the einsum, you can pack T_i into one matrix, and then use tensor notation: docs.scipy.org/doc/numpy-1.14.0/reference/generated/… – anishtain4

1 Answers

0
votes

Keep in mind that to numpy, these 4 and 5d arrays are just that, not 3d arrays containing 2d matrices, etc.

Let's try to write your calculation in a way that clarifies dimensions:

M[x,y,z] = T_1*c[x,y,z,0] + T_2*c[x,y,z,1]...T_5*c[x,y,z,4]

M[x,y,z,:,:] = T_1[:,:]*c[x,y,z,0] + T_2[:,:]*c[x,y,z,1]...T_5[:,:]*c[x,y,z,4]

c[x,y,z,i] is a coefficient, right? So M is a weighted sum of the T_n arrays?

One way of expressing this is:

T = np.stack([T_1, T_2, ...T_5], axis=0)   # 3d  (nab)
M = np.einsum('nab,xyzn->xyzab', T, c)

We could alternatively stack T_i on a new last axis

T = np.stack([T_1, T_2 ...T_5], axis=2)   # (abn)
M = np.einsum('abn,xyzn->xyzab', T, c)

or as broadcasted multiplication plus sum:

M = (T[None,None,None,:,:,:] * c[:,:,:,None,None,:]).sum(axis=-1)

I'm writing this code without testing, so there may be errors, but I think the basic outline is right.

It could also be written as a dot, if I can put the n dimension last in one argument, and 2nd to the last in the other. Or with tensordot. But there's less control over broadcasting of the other dimensions.

For test calculations you could also reshape these arrays so that the x,y,z are rolled into one, and the a,b into another, e.g

 M[xyz,:] = T_n[ab]*c[xyz,n]   # etc