Unnecessary elementary product of a 3d array

I have two 3D arrays A and B with the shape (N, 2, 2), which I would like to multiply by the elements along the N axis with the matrix product by each of the 2x2 matrices. With loop implementation, it looks like

C[i] = dot(A[i], B[i]) 

Is there a way to do this without using a loop? I looked at the tenderdot, but could not get it to work. I think I might need something like tensordot(a, b, axes=([1,2], [2,1])) , but that gives me the NxN matrix.

+8
python arrays vectorization numpy matrix
source share
2 answers

It seems you are doing matrix multiplications for each slice along the first axis. For this you can use np.einsum like this -

 np.einsum('ijk,ikl->ijl',A,B) 

Checking runtime and checking results -

 In [179]: N = 10000 ...: A = np.random.rand(N,2,2) ...: B = np.random.rand(N,2,2) ...: ...: def einsum_based(A,B): ...: return np.einsum('ijk,ikl->ijl',A,B) ...: ...: def forloop(A,B): ...: N = A.shape[0] ...: C = np.zeros((N,2,2)) ...: for i in range(N): ...: C[i] = np.dot(A[i], B[i]) ...: return C ...: In [180]: np.allclose(einsum_based(A,B),forloop(A,B)) Out[180]: True In [181]: %timeit forloop(A,B) 10 loops, best of 3: 54.9 ms per loop In [182]: %timeit einsum_based(A,B) 100 loops, best of 3: 5.92 ms per loop 
+9
source share

You just need to perform the operation on the first measurement of your tensors, which is labeled 0 :

 c = tensordot(a, b, axes=(0,0)) 

This will work as you wish. Also, you do not need a list of axes, because this is just one dimension that you take. With axes([1,2],[2,1]) you cross with the 2nd and 3rd dimensions. If you write it in index notation (Einstein summation convention), this corresponds to c[i,j] = a[i,k,l]*b[j,k,l] , so you shorten the indexes you want to keep .

EDIT: Well, the problem is that the tensor product of two three-dimensional objects is 6d. Because abbreviations are associated with pairs of indices, you cannot get a 3D object using the tensordot operation. The trick is to split your calculation into two parts: first you make tensordot by index to perform the matrix operation, and then you take the diagonal of the tensor to reduce your 4d object to 3d. In one team:

 d = np.diagonal(np.tensordot(a,b,axes=()), axis1=0, axis2=2) 

In the tensor notation, d[i,j,k] = c[i,j,i,k] = a[i,j,l]*b[i,l,k] .

0
source share

All Articles