在numpy中计算矩阵，使每个单元格包含该行其他所有条目的乘积
我有一个矩阵
A = np.array([[0.2, 0.4, 0.6],
[0.5, 0.5, 0.5],
[0.6, 0.4, 0.2]])
我想要一个新的矩阵，这个新矩阵中第 i 行第 j 列的值是原矩阵 A 的第 i 行所有值的乘积，但不包括第 j 列的那个值。
array([[ 0.24, 0.12, 0.08],
[ 0.25, 0.25, 0.25],
[ 0.08, 0.12, 0.24]])
我最开始想到的解决办法是
np.repeat(np.prod(A, 1, keepdims = True), 3, axis = 1) / A
但是这个方法只适用于没有值为零的情况。
有什么想法吗？谢谢！
编辑：我已经开发了
B = np.zeros((3, 3))
for i in range(3):
for j in range(3):
B[i, j] = np.prod(i, A[[x for x in range(3) if x != j]])
但我相信还有更优雅的方法来实现这个，能利用 numpy 高效的 C 后端，而不是低效的 python 循环？
5 个回答
这里有一种O(n^2)的方法，不使用Python的循环或数值近似：
def double_cumprod(A):
B = np.empty((A.shape[0],A.shape[1]+1),A.dtype)
B[:,0] = 1
B[:,1:] = A
L = np.cumprod(B, axis=1)
B[:,1:] = A[:,::-1]
R = np.cumprod(B, axis=1)[:,::-1]
return L[:,:-1] * R[:,1:]
注意：这个方法的速度大约是数值近似方法的两倍慢，这符合我们的预期。
如果你能接受一些小错误的话，可以使用你最开始提的那个解决方案。
A += 1e-10
np.around(np.repeat(np.prod(A, 1, keepdims = True), 3, axis = 1) / A, 9)
这是你解决方案的一种变体，使用了 repeat 和 [:,None]。
np.prod(A,axis=1)[:,None]/A
我第一次尝试处理 0s 的方法是：
In [21]: B
array([[ 0.2, 0.4, 0.6],
[ 0. , 0.5, 0.5],
[ 0.6, 0.4, 0.2]])
In [22]: np.prod(B,axis=1)[:,None]/(B+np.where(B==0,1,0))
array([[ 0.24, 0.12, 0.08],
[ 0. , 0. , 0. ],
[ 0.08, 0.12, 0.24]])
但是正如评论所指出的；[0,1] 这个单元格应该是 0.25。
这个修改解决了那个问题，但现在在有多个连续的 0 时又出现了新问题。
In [30]: I=B==0
In [31]: B1=B+np.where(I,1,0)
In [32]: B2=np.prod(B1,axis=1)[:,None]/B1
In [33]: B3=np.prod(B,axis=1)[:,None]/B1
In [34]: np.where(I,B2,B3)
Out[34]:
array([[ 0.24, 0.12, 0.08],
[ 0.25, 0. , 0. ],
[ 0.08, 0.12, 0.24]])
In [55]: C
array([[ 0.2, 0.4, 0.6],
[ 0. , 0.5, 0. ],
[ 0.6, 0.4, 0.2]])
In [64]: np.where(I,sum1[:,None],sum[:,None])/C1
array([[ 0.24, 0.12, 0.08],
[ 0.5 , 0. , 0.5 ],
[ 0.08, 0.12, 0.24]])
Blaz Bratanic 的 epsilon 方法是目前为止最好的非迭代解决方案：
In [74]: np.prod(C+eps,axis=1)[:,None]/(C+eps)
还有一种不同的解决方案是对列进行迭代：
def paulj(A):
P = np.ones_like(A)
for i in range(1,A.shape[1]):
P *= np.roll(A, i, axis=1)
return P
In [130]: paulj(A)
array([[ 0.24, 0.12, 0.08],
[ 0.25, 0.25, 0.25],
[ 0.08, 0.12, 0.24]])
In [131]: paulj(B)
array([[ 0.24, 0.12, 0.08],
[ 0.25, 0. , 0. ],
[ 0.08, 0.12, 0.24]])
In [132]: paulj(C)
array([[ 0.24, 0.12, 0.08],
[ 0. , 0. , 0. ],
[ 0.08, 0.12, 0.24]])
我在一个大矩阵上做了一些时间测试。
In [13]: A=np.random.randint(0,100,(1000,1000))*0.01
In [14]: timeit paulj(A)
1 loops, best of 3: 23.2 s per loop
In [15]: timeit blaz(A)
10 loops, best of 3: 80.7 ms per loop
In [16]: timeit zwinck1(A)
1 loops, best of 3: 15.3 s per loop
In [17]: timeit zwinck2(A)
1 loops, best of 3: 65.3 s per loop
epsilon 近似方法可能是我们能期待的最佳速度，但它有一些四舍五入的问题。需要对很多列进行迭代会影响速度。我不太明白为什么 np.prod(A[:,mask], 1) 的方法是最慢的。
eeclo https://stackoverflow.com/a/22441825/901925 建议使用 as_strided。我想这就是他所想的（改编自一个重叠块的问题，https://stackoverflow.com/a/8070716/901925）
def strided(A):
h,w = A.shape
A2 = np.hstack([A,A])
x,y = A2.strides
strides = (y,x,y)
shape = (w, h, w-1)
blocks = np.lib.stride_tricks.as_strided(A2[:,1:], shape=shape, strides=strides)
P = blocks.prod(2).T # faster to prod on last dim
# alt: shape = (w-1, h, w), and P=blocks.prod(0)
return P
对于 (1000,1000) 的数组，时间测试比列迭代有了很大改善，尽管仍然比 epsilon 方法慢很多。
In [153]: timeit strided(A)
1 loops, best of 3: 2.51 s per loop
另一种索引方法虽然相对简单，但速度较慢，并且更容易出现内存错误。
def foo(A):
h,w = A.shape
I = (np.arange(w)[:,None]+np.arange(1,w))
I1 = np.array(I)%w
P = A[:,I1].prod(2)
return P
我现在很忙，没时间去解决这个问题；不过我会这样做：在最后一个轴上创建一个连续的循环视图，也就是把数组自己拼接起来，然后用np.lib.index_tricks.as_strided来选择合适的元素进行np.prod运算。这样做不需要用到Python的循环，也不需要数值近似。
补充说明：给你看看：
import numpy as np
A = np.array([[0.2, 0.4, 0.6],
[0.5, 0.5, 0.5],
[0.5, 0.0, 0.5],
[0.6, 0.4, 0.2]])
B = np.concatenate((A,A),axis=1)
C = np.lib.index_tricks.as_strided(
B,
A.shape +A.shape[1:],
B.strides+B.strides[1:])
D = np.prod(C[...,1:], axis=-1)
print D
注意：这种方法并不是最理想的，因为它的复杂度是O(n^3)。你可以看看我另外一个发布的解决方案，复杂度是O(n^2)。
如果你能接受只用一个循环：
B = np.empty_like(A)
for col in range(A.shape[1]):
B[:,col] = np.prod(np.delete(A, col, 1), 1)
这个方法会一次计算你需要的一个列。虽然它的效率没有理论上那么高，因为 np.delete() 会创建一个副本；如果你非常在意内存的使用，建议用掩码来代替：
B = np.empty_like(A)
mask = np.ones(A.shape[1], dtype=bool)
for col in range(A.shape[1]):
mask[col] = False
B[:,col] = np.prod(A[:,mask], 1)
mask[col] = True