Numpy: outer product of n vectors

Question

I'm trying to do something simple in numpy, and I'm sure there should be an easy way of doing it.

Basically, I have a list of n vectors with various lengths. If v1[i] is the i'th entry of the first vector then I want to find a n-dimensional array, A, such that

A[i,j,k...] = v1[i] v2[j] v3[k] ...

My problem is that:

1. outer only takes two vector arguments.

2. einsum requires a parameter like "abcd..." which seems unnecessary.

3. kron requires what seems like rather complex reshaping, and takes only two arguments.

I'd like to avoid as much complexity as possible, so as to avoid introducing bugs. So preferably I would like a single command.

So far, the best I have some up with is:

 vs = [v1, v2, v3 ...]
 shape = map(len, vs)

 # specify the orientation of each vector
 newshapes = diag(array(shape)-1)+1
 reshaped = [x.reshape(y) for x,y in zip(vs, newshapes)]

 # direct product
 A = reduce(lambda a,b: a*b, reshaped, 1)

Answer 1

You use use following one line code:

reduce(np.multiply, np.ix_(*vs))

np.ix_() will do the outer broadcast, you need reduce, but you can pass the ufunc np.multiply without lambda function.

Here is the comparing:

import numpy as np
vs = [np.r_[1,2,3.0],np.r_[4,5.0],np.r_[6,7,8.0]]
shape = map(len, vs)

 # specify the orientation of each vector
newshapes = np.diag(np.array(shape)-1)+1
reshaped = [x.reshape(y) for x,y in zip(vs, newshapes)]

# direct product
A = reduce(lambda a,b: a*b, reshaped, 1)
B = reduce(np.multiply, np.ix_(*vs))

np.all(A==B)

The reuslt:

True

Answer 2

There is an alternative line of code:

reduce(np.multiply.outer, vs)

It is more transparent for me than np.ix_(*vs) construction and support multidimensional arrays like in this question.

Timings are the same within a tolerance:

import numpy as np
from functools import reduce

def outer1(*vs):
    return np.multiply.reduce(np.ix_(*vs))
def outer2(*vs):
    return reduce(np.multiply.outer, vs)

v1 = np.random.randn(100)
v2 = np.random.randn(200)
v3 = np.random.randn(300)
v4 = np.random.randn(50)

%timeit outer1(v1, v2, v3, v4)
# 1 loop, best of 3: 796 ms per loop

%timeit outer2(v1, v2, v3, v4)
# 1 loop, best of 3: 795 ms per loop

np.all(outer1(v1, v2, v3, v4) == outer2(v1, v2, v3, v4))
# True