如何在numpy中使用3向外部產品？

Question

關於numpy.outer [link]。

給定兩個向量，a = [a0, a1, ..., aM]和b = [b0, b1, ..., bN]，外積是M * N矩陣。

1. 但是如何實現一個三陣列外積，這意味着：給定 第三矢量c = [c0, c1, ..., cP]，如何讓外積 3個numpy的陣列之間。

2. 以及如何在numpy中獲得n路陣列的n路外積，方法爲einsum，如何將'i,j,k->ijk'更改爲處理"n"。

Answer 1

你可以爲這個use 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]]]) 

---

你不能越過使用einsum載體的任意量。該指數只能是字母，所以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進制的情況。儘管如此，仍有32個向量的實現限制：

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

---

timeit用3個載體：

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

26個載體：

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

用33個載體：

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