CODEWITHSUNDEEP
numpymatrix-multiplication
tejaswi
tejaswi

Reputation: 33

numpy outerproduct of sequence of arrays

I have a matrix A (nXm) . My ultimate goal is to get Z of dimension (nXmXm) Currently I am doing it using this but can it be done without using for loop using some matrix.tensordot or matrix.multiply.outer

 for i in range(0,A.shape[0]):
       Z[i,:,:] = np.outer(A[i,:],A[i,:])

Upvotes: 3

Views: 42

Answers (2)

Oliver W.
Oliver W.

Reputation: 13459

You could use numpy's Einstein summation, like this:

np.einsum('ij, ik -> ijk', a, a)

Just for completeness, the timing comparison with the also excellent answer (+1) from unutbu:

In [39]: A = np.random.random((1000,50))

In [40]: %timeit using_einsum(A)
100 loops, best of 3: 11.6 ms per loop

In [41]: %timeit using_broadcasting(A)
100 loops, best of 3: 10.2 ms per loop

In [42]: %timeit orig(A)
10 loops, best of 3: 27.8 ms per loop

Which teaches me that

  1. unutbu's machine is faster than mine
  2. broadcasting would be slightly faster than np.einsum

Upvotes: 3

unutbu
unutbu

Reputation: 880489

for i in range(0,A.shape[0]):
    Z[i,:,:] = np.outer(A[i,:],A[i,:])

means

Z_ijk = A_ij * A_ik

which can be computed using NumPy broadcasting:

Z = A[:, :, np.newaxis] * A[:, np.newaxis, :]

A[:, :, np.newaxis] has shape (n, m, 1) and A[:, np.newaxis, :] has shape (n, 1, m). Multiplying the two causes both arrays to be broadcasted up to shape (n, m, m).

NumPy multiplication is always performed elementwise. The values along the broadcasted axis are the same everywhere, so elementwise multiplication results in Z_ijk = A_ij * A_ik.

import numpy as np

def orig(A):
    Z = np.empty(A.shape+(A.shape[-1],), dtype=A.dtype)
    for i in range(0,A.shape[0]):
        Z[i,:,:] = np.outer(A[i,:],A[i,:])
    return Z

def using_broadcasting(A):
    return A[:, :, np.newaxis] * A[:, np.newaxis, :]

Here is a sanity check showing this produces the correct result:

A = np.random.random((1000,50))
assert np.allclose(using_broadcasting(A), orig(A))

By choosing A.shape[0] to be large we get an example which shows off the advantage of broadcasting over looping in Python:

In [107]: %timeit using_broadcasting(A)
10 loops, best of 3: 6.12 ms per loop

In [108]: %timeit orig(A)
100 loops, best of 3: 16.9 ms per loop

Upvotes: 3

Related Questions