numpy的/ Python的：高效矩阵作为输入矩阵的乘积的乘法

Question

问题：

输入是（I，J）矩阵M的期望的输出是（ⅰ^ N，J^n）的矩阵K，其中n是采取的产品数量。冗长的方式来获得所需的输出是如下

• 生成N个行置换I的所有阵列（总的i ** N N阵列）
• 生成N个列置换的所有阵列Ĵ（总的j ** n个数组）
• K [i，j] = m [I [0]，J [0]] * ... * m [I [n]，J [n]]对于所有n （len（J））

我这样做的直接方法是通过生成范围内所有数字n置换的标签列表（len（np.shape（m）[0] ））和范围（len（np.shape（m）[1]））分别为行和列。之后，您可以像上面最后一个项目符号点那样放大它们。然而，这对于大型输入矩阵并不实用 - 所以我正在寻找优化上述方法。预先感谢您

实施例：

输入

np.array（[[1,2,3]，[4,5,6]]）

输出对于n = 3

[[1. 2. 3. 2. 4. 6. 3. 6. 9. 2. 4. 6. 4. 8. 12. 6. 12. 18. 3 6. 9. 6. 12. 18. 9. 18. 27.]

[4. 5. 6. 8. 10. 12. 12. 15. 18. 8. 10. 12. 16. 20. 24. 24. 30. 36。 12. 15. 18. 24. 30. 36. 36. 45. 54.]

[4. 8. 12. 5. 10. 15. 6. 12. 18. 8. 16. 24. 10. 20. 30. 12. 24. 36. 12. 24. 36. 15. 30. 45. 18. 36. 54.]

[16. 20. 24. 20. 25. 30. 24 30. 36. 32. 40. 48. 40. 50. 60. 48. 60. 72. 48. 60. 72. 60. 75. 90. 72. 90. 108.]

[4. 8. 12. 8. 16. 24. 12. 24. 36. 5. 10. 15. 10. 20. 30. 15. 30. 45. 6. 12. 18. 12. 24. 36. 18. 36. 54.]

[16. 20. 24. 32. 48. 48. 60. 72. 20. 25. 30. 40. 50. 60. 60. 75。 90. 24. 30. 36. 48. 60. 72. 72. 90. 108.]

[16. 32. 48. 20. 40. 60. 24. 48. 72. 20. 40. 60 25. 50. 75. 30. 60. 90. 24. 48. 72. 30. 60. 90. 36. 72. 108.]

[64. 80. 96. 80. 100. 120 96. 120. 144. 80. 100. 120. 100. 125. 150. 120. 150. 180. 96. 120. 144. 120. 150. 180. 144. 180. 216。]

偏溶液：

我发现的最好的是一个函数来创建矩阵的笛卡尔乘积这里提出：https://*.com/a/1235363/4003747 的问题是，该输出不是矩阵而是阵列阵列。乘以每个数组的元素给出我后面的值，但以无序的方式。我尝试了一段时间，但我不知道如何明智地重新排列它们。

对于n低效溶液= 3：

import numpy as np 
import itertools 

m=np.array([[1,2,3],[4,5,6]]) 

def f(m): 
    labels_i = [list(p) for p in itertools.product(range(np.shape(m)[0]),repeat=3)] 
    labels_j = [list(p) for p in itertools.product(range(np.shape(m)[1]),repeat=3)] 


    out = np.zeros([len(labels_i),len(labels_j)]) 
    for i in range(len(labels_i)): 
     for j in range(len(labels_j)): 
      out[i,j] = m[labels_i[i][0],labels_j[j][0]] * m[labels_i[i][1],labels_j[j][1]] * m[labels_i[i][2],labels_j[j][2]] 
    return out 

Answer 1

下面是使用和的broadcasting组合linear indexing向量化的方法 -

from itertools import product 

# Get input array's shape 
r,c = A.shape 

# Setup arrays corresponding to labels i and j 
arr_i = np.array(list(product(range(r), repeat=n))) 
arr_j = np.array(list(product(range(c), repeat=n))) 

# Use linear indexing with ".ravel()" to extract elements. 
# Perform elementwise product along the rows for the final output 
out = A.ravel()[(arr_i*c)[:,None,:] + arr_j].prod(2) 

运行时测试并输出验证 -

In [167]: # Inputs   
    ...: n = 4 
    ...: A = np.array([[1,2,3],[4,5,6]]) 
    ...: 
    ...: def f(m): 
    ...: labels_i = [list(p) for p in product(range(np.shape(m)[0]),repeat=n)] 
    ...: labels_j = [list(p) for p in product(range(np.shape(m)[1]),repeat=n)] 
    ...: 
    ...: out = np.zeros([len(labels_i),len(labels_j)]) 
    ...: for i in range(len(labels_i)): 
    ...:  for j in range(len(labels_j)): 
    ...:   out[i,j] = m[labels_i[i][0],labels_j[j][0]] \ 
    ...:     * m[labels_i[i][1],labels_j[j][1]] \ 
    ...:     * m[labels_i[i][2],labels_j[j][2]] \ 
    ...:     * m[labels_i[i][3],labels_j[j][3]] 
    ...: return out 
    ...: 
    ...: def f_vectorized(A,n): 
    ...: r,c = A.shape 
    ...: arr_i = np.array(list(product(range(r), repeat=n))) 
    ...: arr_j = np.array(list(product(range(c), repeat=n))) 
    ...: return A.ravel()[(arr_i*c)[:,None,:] + arr_j].prod(2) 
    ...: 

In [168]: np.allclose(f_vectorized(A,n),f(A)) 
Out[168]: True 

In [169]: %timeit f(A) 
100 loops, best of 3: 2.37 ms per loop 

In [170]: %timeit f_vectorized(A,n) 
1000 loops, best of 3: 202 µs per loop 

Answer 2

这应该工作：

import numpy as np 
import itertools 

m=np.array([[1,2,3],[4,5,6]]) 
n=3 # change your n here 


def f(m): 
    labels_i = [list(p) for p in itertools.product(range(np.shape(m)[0]),repeat=n)] 
    labels_j = [list(p) for p in itertools.product(range(np.shape(m)[1]),repeat=n)] 


    out = np.zeros([len(labels_i),len(labels_j)]) 
    for i in range(len(labels_i)): 
     for j in range(len(labels_j)): 
      out[i,j] = np.prod([m[labels_i[i][k],labels_j[j][k]] for k in range(n)]) 
    return out