Multiplying every element of one array by every element of another array
Asked Answered
L

5

6

Say I have two arrays,

import numpy as np


x = np.array([1, 2, 3, 4])
y = np.array([5, 6, 7, 8])

What's the fastest, most Pythonic, etc., etc. way to get a new array, z, with a number of elements equal to x.size * y.size, in which the elements are the products of every pair of elements (x_i, y_j) from the two input arrays.

To rephrase, I'm looking for an array z in which z[k] is x[i] * y[j].

A simple but inefficient way to get this is as follows:

z = np.empty(x.size * y.size)
counter = 0
for i in x:
    for j in y:
        z[counter] = i * j
        counter += 1

Running the above code shows that z in this example is

In [3]: z
Out[3]: 
array([  5.,   6.,   7.,   8.,  10.,  12.,  14.,  16.,  15.,  18.,  21.,
        24.,  20.,  24.,  28.,  32.])
Lawn answered 2/6, 2015 at 4:0 Comment(1)
@dbliss: Could've been an outer product, could've been a convolution, could've been something else. Outer product looked most likely, but there wasn't enough information to be sure.Taxidermy
N
3

Two more approaches could be suggested here.

Using matrix-multiplication with np.dot:

np.dot(x[:,None],y[None]).ravel()

With np.einsum:

np.einsum('i,j->ij',x,y).ravel()

Runtime tests

In [31]: N = 10000
    ...: x = np.random.rand(N)
    ...: y = np.random.rand(N)
    ...: 

In [32]: %timeit np.dot(x[:,None],y[None]).ravel()
1 loops, best of 3: 302 ms per loop

In [33]: %timeit np.einsum('i,j->ij',x,y).ravel()
1 loops, best of 3: 274 ms per loop

Same as @BilalAkil's answer but with ravel() instead of flatten() as a faster alternative -

In [34]: %timeit np.multiply.outer(x, y).ravel() 
1 loops, best of 3: 211 ms per loop

@BilalAkil's answer:

In [35]: %timeit np.multiply.outer(x, y).flatten()
1 loops, best of 3: 451 ms per loop

@Tim Leathart's answer:

In [36]: %timeit np.array([y * a for a in x]).flatten()
1 loops, best of 3: 766 ms per loop
Newsstand answered 2/6, 2015 at 5:30 Comment(2)
odd that ravel is faster than flatten. (i'd expect them to be the same.) can you explain that time difference?Lawn
@dbliss It has to do with the fact that flatten() produces a copy, whereas ravel() is just a view into the input array, so it's cheaper. To verify this, you can create a toy setup : A = np.random.rand(3,4) and then check for view or not with np.may_share_memory(). So, np.may_share_memory(A.ravel(),A) would be True, whereas np.may_share_memory(A.flatten(),A) would result in False. More on checking view or not, here: https://mcmap.net/q/77285/-how-can-i-tell-if-numpy-creates-a-view-or-a-copyNewsstand
S
5

Well I haven't much experience with numpy, but a quick search gave me this: http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.outer.html

>>> np.multiply.outer([1, 2, 3], [4, 5, 6])
array([[ 4,  5,  6],
   [ 8, 10, 12],
   [12, 15, 18]])

You can then flatten that array to get the same output as you requested: http://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.flatten.html

EDIT: @Divakar's answer showed us that ravel will do the same thing as flatten, except faster o.O So use that instead.

So in your case, it'd look like this:

>>> np.multiply.outer(x, y).ravel()

BONUS: You can go multi-dimensional with this!

Shul answered 2/6, 2015 at 4:15 Comment(0)
L
4

Here's one way to do it:

import itertools


z = np.empty(x.size * y.size)
counter = 0
for i, j in itertools.product(x, y):
    z[counter] = i * j
    counter += 1

It'd be nice to get rid of that counter, though, as well as the for loop (but at least I got rid of one of the loops).

UPDATE

Being one-liners, the other provided answers are better than this one (according to my standards, which value brevity). The timing results below show that @BilalAkil's answer is faster than @TimLeathart's:

In [10]: %timeit np.array([x * j for j in y]).flatten()
The slowest run took 4.37 times longer than the fastest. This could mean that an intermediate result is being cached 
10000 loops, best of 3: 24.2 µs per loop

In [11]: %timeit np.multiply.outer(x, y).flatten()
The slowest run took 5.59 times longer than the fastest. This could mean that an intermediate result is being cached 
100000 loops, best of 3: 10.5 µs per loop
Lawn answered 2/6, 2015 at 4:15 Comment(3)
Although you got rid of one of the loops, you're still performing a cross multiplication (in this case with itertools.product), so it's still got the same complexity as having both loops: O(xy). It'll be difficult to avoid that.Shul
Try timing your basic loop implementations too. The one with your nested loops and the itertools.product one :PShul
@BilalAkil i'd rather not embarrass myself ;)Lawn
N
3

Two more approaches could be suggested here.

Using matrix-multiplication with np.dot:

np.dot(x[:,None],y[None]).ravel()

With np.einsum:

np.einsum('i,j->ij',x,y).ravel()

Runtime tests

In [31]: N = 10000
    ...: x = np.random.rand(N)
    ...: y = np.random.rand(N)
    ...: 

In [32]: %timeit np.dot(x[:,None],y[None]).ravel()
1 loops, best of 3: 302 ms per loop

In [33]: %timeit np.einsum('i,j->ij',x,y).ravel()
1 loops, best of 3: 274 ms per loop

Same as @BilalAkil's answer but with ravel() instead of flatten() as a faster alternative -

In [34]: %timeit np.multiply.outer(x, y).ravel() 
1 loops, best of 3: 211 ms per loop

@BilalAkil's answer:

In [35]: %timeit np.multiply.outer(x, y).flatten()
1 loops, best of 3: 451 ms per loop

@Tim Leathart's answer:

In [36]: %timeit np.array([y * a for a in x]).flatten()
1 loops, best of 3: 766 ms per loop
Newsstand answered 2/6, 2015 at 5:30 Comment(2)
odd that ravel is faster than flatten. (i'd expect them to be the same.) can you explain that time difference?Lawn
@dbliss It has to do with the fact that flatten() produces a copy, whereas ravel() is just a view into the input array, so it's cheaper. To verify this, you can create a toy setup : A = np.random.rand(3,4) and then check for view or not with np.may_share_memory(). So, np.may_share_memory(A.ravel(),A) would be True, whereas np.may_share_memory(A.flatten(),A) would result in False. More on checking view or not, here: https://mcmap.net/q/77285/-how-can-i-tell-if-numpy-creates-a-view-or-a-copyNewsstand
S
2

Here's a way to do it:

import numpy as np

x = np.array([1, 2, 3, 4])
y = np.array([5, 6, 7, 8])
z = np.array([y * a for a in x]).flatten()
Scare answered 2/6, 2015 at 4:15 Comment(0)
B
0

I know I'm super late to the party here, but I thought I'd throw my hat into the ring for anyone reading this question in the future. Using the same metric as @Divakar, I added what I consider to be a much more intuitive solution to the list (the first code snippet measured):

import numpy as np
N = 10000
x = np.random.rand(N)
y = np.random.rand(N)

%timeit np.ravel(x[:,None] * y[None])
635 ms ± 19.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit np.outer(x, y).ravel()
640 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit np.dot(x[:,None],y[None]).ravel()
853 ms ± 57.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit np.einsum('i,j->ij',x,y).ravel()
754 ms ± 19.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Based on the similarity in execution time, it seems likely that numpy.outer functions exactly the same way as my solution internally, although you should take observations like this with a hefty grain of salt.

The reason why I find it more intuitive is that, unlike all other solutions, its syntax isn't strictly limited to multiplication. For example, np.ravel(x[:,None] / y[None]) will give you a / b for every a in x and b in y.

Bonkers answered 8/6, 2020 at 0:15 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.