1

Suppose array_1 and array_2 are two arrays of matrices of the same sizes. Is there any vectorised way of multiplying element-wise, the elements of these two arrays(which their elements' multiplication is well defined)?

The dummy code:

def mat_multiply(array_1,array_2): size=np.shape(array_1)[0] result=np.array([]) for i in range(size): result=np.append(result,np.dot(array_1[i],array_2[i]),axis=0) return np.reshape(result,(size,2))

example input:

a=[[[1,2],[3,4]],[[1,2],[3,4]]] b=[[1,3],[4,5]]

output:

[[ 7. 15.] [ 14. 32.]]
Improve this question
2

2 Answers 2

Reset to default
4

Contrary to your first sentence, a and b are not the same size. But let's focus on your example.

So you want this - 2 dot products, one for each row of a and b

np.array([np.dot(x,y) for x,y in zip(a,b)])

or to avoid appending

X = np.zeros((2,2)) for i in range(2): X[i,...] = np.dot(a[i],b[i])

the dot product can be expressed with einsum (matrix index notation) as

[np.einsum('ij,j->i',x,y) for x,y in zip(a,b)]

so the next step is to index that first dimension:

np.einsum('kij,kj->ki',a,b)

I'm quite familiar with einsum, but it still took a bit of trial and error to figure out what you want. Now that the problem is clear I can compute it in several other ways

A, B = np.array(a), np.array(b) np.multiply(A,B[:,np.newaxis,:]).sum(axis=2) (A*B[:,None,:]).sum(2) np.dot(A,B.T)[0,...] np.tensordot(b,a,(-1,-1))[:,0,:]

I find it helpful to work with arrays that have different sizes. For example if A were (2,3,4) and B (2,4), it would be more obvious the dot sum has to be on the last dimension.

Another numpy iteration tool is np.nditer. einsum uses this (in C). http://docs.scipy.org/doc/numpy/reference/arrays.nditer.html

it = np.nditer([A, B, None],flags=['external_loop'], op_axes=[[0,1,2], [0,-1,1], None]) for x,y,w in it: # x, y are shape (2,) w[...] = np.dot(x,y) it.operands[2][...,0]

Avoiding that [...,0] step, requires a more elaborate setup.

C = np.zeros((2,2)) it = np.nditer([A, B, C],flags=['external_loop','reduce_ok'], op_axes=[[0,1,2], [0,-1,1], [0,1,-1]], op_flags=[['readonly'],['readonly'],['readwrite']]) for x,y,w in it: w[...] = np.dot(x,y) # w[...] += x*y print C # array([[ 7., 15.],[ 14., 32.]])
Improve this answer
Sign up to request clarification or add additional context in comments.

5 Comments

thanks. vote up! But they are the same size! a consists of two matrices: [[1,2],[3,4]] and [[1,2],[3,4]]. And b consists of two matrices:[1,3] and [4,5].
I think you can replace the k in your np.einsum' call with ...` and get it to work for broadcastable arrays of matrices and vectors, not just 1D arrays.
Naji - your a and b are lists of lists. A is a 3d array, B 2d. If you want to vectorize operations, you need to think in terms of these higher dimensional arrays. Otherwise you are stuck with iterating over lower dimension arrays.
@hpaulj, I understand your point. You mentioned that I was being wrong by saying that a nd b are of the same size. The point is I was talking in terms of mathematics and in that terms a and b are of the same size because they include the same size of matrices. And what I am expecting is to get element-wise multiplication between elements(each matrix in this array) of these arrays(whenever the multiplication between the paired elements are well defined mathematically). My only point was to clarify that I was not talking about data-structures but mathematical k-tumples of matrices.
A concept that is used in numpy documentation is compatible shapes, docs.scipy.org/doc/numpy/user/basics.broadcasting.html
2

There's one more option that @hpaulj left out in his extensive and comprehensive list of options:

>>> a = np.array(a) >>> b = np.array(b) >>> from numpy.core.umath_tests import matrix_multiply >>> matrix_multiply.signature '(m,n),(n,p)->(m,p)' >>> matrix_multiply(a, b[..., np.newaxis]) array([[[ 7], [15]], [[14], [32]]]) >>> matrix_multiply(a, b[..., np.newaxis]).shape (2L, 2L, 1L) >>> np.squeeze(matrix_multiply(a, b[..., np.newaxis]), axis=-1) array([[ 7, 15], [14, 32]])

The nice thing about matrix_multiply is that, it being a gufunc, it will work not only with 1D arrays of matrices, but also with broadcastable arrays. As an example, if instead of multiplying the first matrix with the first vector, and the second matrix with the second vector, you wanted to compute all possible multiplications, you could simply do:

>>> a = np.arange(8).reshape(2, 2, 2) # to have different matrices >>> np.squeeze(matrix_multiply(a[...,np.newaxis, :, :], ... b[..., np.newaxis]), axis=-1) array([[[ 3, 11], [ 5, 23]], [[19, 27], [41, 59]]])
Improve this answer

Comments

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.