Numpy数组的计算/操作

Question

希望能够尽快进行此计算。我有一个n x m的numpy数组X。我想定义Y如下：

Y_11 = 1 / (exp(X_11-X_11) + exp(X_11-X_12) + ... exp(X_11 - X_1N) ).

或者对于Y_00

1/np.sum(np.exp(X[0,0]-X[0,:]))

所以基本上，Y也是n x m的，其中i，j元素为1 / sum_j' exp(X_ij - X_ij')

任何提示都将很棒！谢谢。

根据要求提供示例代码：

np.random.seed(111)
J,K = 110,120
X = np.random.rand(J,K)
Y = np.zeros((J,K))
for j in range(J):
    for k in range(K):
        Y[j,k] = 1/np.sum(np.exp(X[j,k]-X[j,:]))

# note each row will sum to 1 under this operation
np.sum(Y,axis=1)

- Kevin

---

你的第二个表达式中 exp 的调用去哪了？ - user2357112

谢谢你发现我漏掉了exp()。我还添加了一些代码。 - Kevin

这是一个好问题。将来，如果您能像现在这样提前写好东西，那就太棒了。 - Ami Tavory

Answer 1

这里有几种全矢量化的方法：

def vectorized_app1(X):
    return 1/np.exp(X[:,None] - X[...,None]).sum(1)

def vectorized_app2(X):
    exp_vals = np.exp(X)
    return 1/(exp_vals[:,None]/exp_vals[...,None]).sum(1)

def vectorized_app3(X):
    exp_vals = np.exp(X)
    return 1/np.einsum('ij,ik->ij',exp_vals,1/exp_vals)

使用的技巧和学到的经验

• 通过使用None/np.newaxis扩展维度，引入广播并以向量化的方式执行所有操作。

• np.exp是一项昂贵的操作。因此，在广播庞大数组上使用它将代价高昂。所以，所使用的技巧是：exp(A-B) = exp(A)/exp(B)。因此，我们先执行np.exp(X)，然后执行广播除法。

• 这些求和缩减可以使用np.einsum实现。这带来了内存效率，因为我们不必创建巨大的广播数组，这会极大地提高性能。

运行时间测试-

In [111]: # Setup input, X
     ...: np.random.seed(111)
     ...: J,K = 110,120
     ...: X = np.random.rand(J,K)
     ...: 

In [112]: %timeit original_approach(X)
1 loop, best of 3: 322 ms per loop

In [113]: %timeit vectorized_app1(X)
10 loops, best of 3: 124 ms per loop

In [114]: %timeit vectorized_app2(X)
100 loops, best of 3: 14.6 ms per loop

In [115]: %timeit vectorized_app3(X)
100 loops, best of 3: 3.01 ms per loop

看起来，einsum 再次以 100x+ 的速度展示了它的魔力！

- Divakar

---

非常有用的答案。感谢您提供的所有解释。 - Ami Tavory

谢谢！不过有一个问题，有时候A或B的值很大，所以我会得到一个错误。我认为采用exp（A-B）的方法是我将它们缩小的好办法，并且可以放入np.min（A-B，200）或类似的数据来避免错误。还有其他想法吗？ - Kevin

@Kevin 究竟是什么错误？它说了什么？ - Divakar

指数溢出。我刚刚发现了bigfloat，但需要花时间研究如何在我的代码中实现 - 并查看它是否与knitro兼容。 - Kevin

@Kevin 嗯，我猜你可以在这些情况下继续使用 vectorized_app1。 - Divakar

Answer 2

这里是减少双重循环的第一步：

def foo2(X):
    Y = np.zeros_like(X)
    for k in range(X.shape[1]):
        Y[:,k]=1/np.exp(X[:,[k]]-X[:,:]).sum(axis=1)
    return Y

我猜我也可以去掉 k 循环，但我需要更多的时间来弄清楚它具体的功能。在整个程序中，X[:,[k]]-X[:,:] 不是很明显。

另一个步骤：

Z = np.stack([X[:,[k]]-X for k in range(X.shape[1])],2)
Y = 1/np.exp(Z).sum(axis=1)

这可以进一步细化（避免过多的试错），使用

Z = X[:,None,:]-X[:,:,None]

- hpaulj

---

很好的回答。非常小的评论：在您的最后一个括号中，您可能是指“没有太多的试错”。 - Ami Tavory

不，我尝试了许多替代方案。X可以以多种方式扩展为3D，并且可以在3个轴上进行求和。我没有像我本可以做的那样系统化。 - hpaulj