3

About the numpy.outer [link] .

Given two vectors, a = [a0, a1, ..., aM] and b = [b0, b1, ..., bN], the outer product will be M*N matrix.

  1. But how to implement a 3-array outer product, which means : given third vector c = [c0, c1, ..., cP], how to get the outer product between the 3 numpy arrays.

  2. and how to get n-way outer product for n-array in numpy, for the method of einsum, how to change 'i,j,k->ijk' to process "n" .

Alice Smith
  • 41
  • 5

3 Answers3

11

The direct way of doing this, taking full advantage of broadcasting is:

a[:,None,None] * b[None,:,None] * c[None,None,:]

np.ix_ does this reshaping for you, at a modest cost in speed

In [919]: np.ix_(a,b,c)
Out[919]: 
(array([[[0]],

        [[1]],

        [[2]],

        [[3]],

        [[4]]]), array([[[10],
         [11],
         [12],
         [13]]]), array([[[20, 21, 22]]]))

and the resulting arrays can be multiplied with

np.prod(np.ix_(a,b,c))

The einsum version is simple, and fast

np.einsum('i,j,k',a,b,c)

It's a good idea to learn all 3 methods.

The problem with nesting outer is that expects the inputs to be 1d, or it flattens them. It can be used, but needs some reshaping

np.outer(a,np.outer(b,c)).reshape(a.shape[0],b.shape[0],c.shape[0])
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • `np.prod(np.ix_(a,b,c))` returns `ValueError ... The requested array has an inhomogeneous shape` as of Python 3.10.9 and Numpy 1.24.2. Using `math.prod` instead of `np.prod` works. – Giancarlo Sportelli Mar 10 '23 at 16:25
  • @GiancarloSportelli, that's because `np.prod` is first making an array from the tuple, and new numpy versions have tighten up that process. `np.prod(np.array((a,b,c)). dtype=object))` should run ok. The change was in the `np.array...` step, not the `np.prod`. `math.prod` works because it iterates on the tuple directly. Both do a repeated product, `(A*B) * C`. – hpaulj Mar 10 '23 at 16:33
5

You can use numpy.ix_; it adds axes to its operands such that they form an open grid. Below the shapes of A,B,C are (2, 1, 1), (1, 3, 1) and (1, 1, 4), so just multiplying them together results in the outer product.

a = np.arange(1, 3)
b = np.arange(1, 4)
c = np.arange(1, 5)
A,B,C = np.ix_(a,b,c)
A*B*C
# array([[[ 1,  2,  3,  4],
#         [ 2,  4,  6,  8],
#         [ 3,  6,  9, 12]],
#
#        [[ 2,  4,  6,  8],
#         [ 4,  8, 12, 16],
#         [ 6, 12, 18, 24]]])
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
4

You could use einsum for this:

>>> numpy.einsum('i,j,k->ijk', [1, 2], [3, 4, 5], [6,7,8,9])
array([[[18, 21, 24, 27],
        [24, 28, 32, 36],
        [30, 35, 40, 45]],

       [[36, 42, 48, 54],
        [48, 56, 64, 72],
        [60, 70, 80, 90]]])

You cannot cross arbitrary amount of vectors using einsum. The indices can only be letters, so at most 52 vectors are allowed in this form (although it will raise "too many operands" when 32 vectors are provided):

def nary_outer_einsum_52(*vectors):
    import string
    subscripts = string.ascii_letters[:len(vectors)]
    subscripts = ','.join(subscripts) + '->' + subscripts
    # expands to `numpy.einsum('a,b,c,d,e->abcde', v[0], v[1], v[2], v[3], v[4])`
    return numpy.einsum(subscripts, *vectors)

einsum supports another form which uses numeric indices instead of letters, which unfortunately only support up to 26 vectors because it just translates to the letter version internally. I do not recommend using this until the bug I mentioned is fixed.

def nary_outer_einsum_26(*vectors):
    operations = (x for i, v in enumerate(vectors) for x in (v, [i]))
    # expands to `numpy.einsum(v[0], [0], v[1], [1], v[2], [2], v[3], [3], v[4], [4])`
    return numpy.einsum(*operations)

numpy.ix_ (@PaulPanzer's answer) can easily be scaled up to the N-ary case. Still, there is an implementation-limit of 32 vectors:

def nary_outer_ix(*vectors):
    return numpy.prod(numpy.ix_(*vectors), axis=0)

timeit with 3 vectors:

>>> timeit.Timer('nary_outer_einsum_52([1,2], [3,4,5], [6,7,8,9])', globals=globals()).autorange()
(100000, 0.8991979530110257)   # 9 µs / iteration
>>> timeit.Timer('nary_outer_einsum_26([1,2], [3,4,5], [6,7,8,9])', globals=globals()).autorange()
(100000, 1.0089023629989242)   # 10 µs / iteration
>>> timeit.Timer('nary_outer_ix([1,2], [3,4,5], [6,7,8,9])', globals=globals()).autorange()
(10000, 0.30978472300921567)   # 31 µs / iteration

with 26 vectors:

>>> timeit.Timer('nary_outer_einsum_52(*([[1]] * 26))', globals=globals()).autorange()
(10000, 0.6589978060073918)   # 66 µs / iteration
>>> timeit.Timer('nary_outer_einsum_26(*([[1]] * 26))', globals=globals()).autorange()
(10000, 0.7502327310066903)   # 75 µs / iteration
>>> timeit.Timer('nary_outer_ix(*([[1]] * 26))', globals=globals()).autorange()
(1000, 0.2026848039968172)    # 203 µs / iteration

with 33 vectors:

>>> nary_outer_einsum_52(*([[1]] * 33))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 6, in nary_outer_einsum_52
  File "/usr/local/lib/python3.6/site-packages/numpy/core/einsumfunc.py", line 948, in einsum
    return c_einsum(*operands, **kwargs)
ValueError: too many operands
>>> nary_outer_ix(*([[1]] * 33))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 2, in nary_outer_ix
  File "/usr/local/lib/python3.6/site-packages/numpy/lib/index_tricks.py", line 83, in ix_
    new = new.reshape((1,)*k + (new.size,) + (1,)*(nd-k-1))
ValueError: sequence too large; cannot be greater than 32
Community
  • 1
  • 1
kennytm
  • 510,854
  • 105
  • 1,084
  • 1,005