numpyで3ウェイの外側の製品はどうですか？

Question

numpy.outerについて [link]について

a = [a0, a1, ..., aM]およびb = [b0, b1, ..., bN]の2つのベクトルが与えられると、外積はM * N行列になります。

1. しかし、どの手段3アレイ外積実装する： 第ベクターc = [c0, c1, ..., cP]与え、3つのnumpyのアレイ間の外積 を取得する方法を。

2. 及び方法"n"を処理する'i,j,k->ijk'を変更する方法、einsumの方法のために、numpyの中のn-アレイのNウェイ外積を取得します。

Answer 1

あなたはこのためにuse einsumことができます：あなたはeinsumを使用してベクトルの任意の金額を越えることはできません

>>> 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]]]) 

---

。インデックスはだけなのでat most 52 vectors are allowed in this form（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は残念ながらonly support up to 26 vectors数値インデックスの代わりの文字を使用する別の形式をサポートしています。私はしません。は、上記のバグが修正されるまでこれを使用することをお勧めします。

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）を容易にN値の場合にスケールアップすることができます。 33のベクターと

>>> 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 

：26のベクターで

>>> 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 

：3つのベクターと

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

---

はtimeit：まだ、32のベクターの実装限界がある

>>> 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 

Answer 2

numpy.ix_を使用できます。軸をオペランドに追加して、開いたグリッドを形成します。 A、B、Cの形状の下は(2, 1, 1), (1, 3, 1)と(1, 1, 4)なので、それらを掛け合わせるだけで外側の製品になります。放送をフルに活用して

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]]]) 

Answer 3

これを行うための直接的な方法は、次のとおりです。

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

np.ix_はスピードで控えめなコストで、あなたのため

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

     [[1]], 

     [[2]], 

     [[3]], 

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

と、この整形を行います結果として生じる配列は、乗算することができます

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

einsumバージョンはシンプルで、かつ高速

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

それはすべての3つの方法を学ぶことをお勧めします。

入れ子の問題outerは、入力を1dにするか、または平らにすることを想定しています。それは使用することができますが、いくつかの形を変える必要があります

np.outer(a,np.outer(b,c)).reshape(a.shape[0],b.shape[0],c.shape[0])