在 numpy 中，计算一个矩阵，其中每个单元格包含该行中所有其他条目的乘积答案

【问题标题】：In numpy, calculating a matrix where each cell contains the product of all the other entries in that row在 numpy 中，计算一个矩阵，其中每个单元格包含该行中所有其他条目的乘积
【发布时间】：2014-04-21 12:36:40
【问题描述】：

我有一个矩阵

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 循环？

【问题讨论】：

标签： python math numpy matrix

【解决方案1】：

如果你愿意容忍一个循环：

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

【讨论】：

+1，但在内存方面，两者是等价的。使用布尔掩码进行索引也会创建一个副本。第二个版本没有优势。（布尔索引是“花式”索引的一种形式，它总是会生成一个副本。与普通切片不同，没有办法用偏移量和步幅来描述结果。）
@JoeKington：你说得对，谢谢。你知道是否有一种更聪明的方法可以按照我们需要的方式对列进行子集化？我们可能希望避免复制，但我们需要为任意列表达“除一个之外的所有列”的想法。我不知道偏移量和步幅如何实现，但也许您知道其他可以节省内存的方法？或者也许今天在 NumPy 中没有办法。

【解决方案2】：

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

【讨论】：

我认为如果有零，它确实工作。第一行的第 0 个元素应该是 0.25，不是吗？
我已经纠正了这个问题，但解决方案比我想要的要复杂。
不幸的是，我认为这也行不通。假设一行中有两个零：那么每个输出行都应该是零，因为每个产品都包含一个零。但是您的新方法会将零转换为非零元素的乘积（因为它们会采用零替换为 1 的情况。）如果您使用 4x4 矩阵进行测试，这更容易看出，它有具有两个非零元素和两个零的行。
[“每个输出行”当然是指“输出行的每个元素”，而不是整个输出数组应该为零。]
hpaulj：没有注意到您已经包含了我的想法的一个版本，但我也对我的帖子进行了编辑。你的版本不完全正确；但无论如何，巨大的速度差异让我有些惊讶。确实 prod 被接管了一个小步幅轴，因此我们应该获得良好的缓存行为，因此唯一的额外成本是连接数组；然而，epsilon 方法除了实际的 prod 之外，还对数组进行了三遍传递......我不太明白。

【解决方案3】：

我正在运行，所以我没有时间制定这个解决方案；但是 id 所做的是在最后一个轴上创建一个连续的圆形视图，方法是将数组沿最后一个轴连接到自身，然后使用 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)

【讨论】：

这与我的跨步版本运行相同。
但是np.allclose(strided(A),double_cumprod(A)) 返回True。
确实如此，但大矩阵会导致到处都是零。 ~0.5**1000 始终使您的浮点数下溢。尝试打印一个小矩阵的结果。

【解决方案4】：

如果您愿意容忍小错误，您可以使用您最初提出的解决方案。

A += 1e-10
np.around(np.repeat(np.prod(A, 1, keepdims = True), 3, axis = 1) / A, 9)

【讨论】：

【解决方案5】：

这是一个没有 python 循环或数值近似的 O(n^2) 方法：

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

注意：它似乎比数值逼近法慢两倍左右，符合预期。

【讨论】：

在我的大矩阵上，这个测试比 Blaz 的近似值快一点。
我也对其进行了测试，但使用的是浮点 1000x1000 矩阵。在这种情况下，近似值在我的机器上会领先，但对于整数，它们的速度大致相同。
当问题被表述为prod(A[j,:i-1])*prod(A[j,i+1:]) 时，这种 cumprod 方法变得更加明显。也就是说，2 个序列的乘积。` 如最初所述，问题的重点是从单个 prod() 中消除 A[j,i]。