用广播替换轴上的循环，第 2 部分

Question

之前我问过一个类似的，其中的答案使用了 np.dot，利用了点积涉及乘积总和的事实。（据我了解。）

现在我有一个类似的问题，我认为 dot 不适用，因为我想取元素对角线来代替总和。如果是的话，我还没能正确应用它。

给定矩阵 x 和数组 err:

x = np.matrix([[ 0.02984406, -0.00257266],
               [-0.00257266,  0.00320312]])

err = np.array([  7.6363226 ,  13.16548267])

我当前的循环实现是：

res = np.array([np.sqrt(np.diagonal(x * err[i])) for i in range(err.shape[0])])

print(res)
[[ 0.47738755  0.15639712]
 [ 0.62682649  0.20535487]]

对于err中的每个i取x.dot(i)的对角线。这可以矢量化吗？换句话说，x * err 的输出是否可以是 3 维的，然后 np.diagonal 产生一个二维数组，每条对角线有一个元素？

Answer 1

节目：

import numpy as np
x = np.matrix([[ 0.02984406, -0.00257266],
               [-0.00257266,  0.00320312]])

err = np.array([  7.6363226 ,  13.16548267])
diag = np.diagonal(x)
ans = np.sqrt(diag*err[:,np.newaxis])  # sqrt of outer product
print(ans)

# use out keyword to avoid making new numpy array for many times.
ans = np.empty(x.shape, dtype=x.dtype)
for i in range(100):
    ans = np.multiply(diag, err, out=ans)
    ans = np.sqrt(ans, out=ans)

结果：

[[ 0.47738755  0.15639712]
 [ 0.62682649  0.20535487]]

Answer 2

这是一种使用 diagonal-view 和 ndarray.flat 进入 x 然后使用 broadcasting 进行元素乘法的方法，就像这样 -

np.sqrt(x.flat[::x.shape[1]+1].A1 * err[:,None])

样本运行-

In [108]: x = np.matrix([[ 0.02984406, -0.00257266],
     ...:                [-0.00257266,  0.00320312]])
     ...: 
     ...: err = np.array([  7.6363226 ,  13.16548267])
     ...: 

In [109]: np.sqrt(x.flat[::x.shape[1]+1].A1 * err[:,None])
Out[109]: 
array([[ 0.47738755,  0.15639712],
       [ 0.62682649,  0.20535487]])

运行时测试以查看 view 如何帮助创建副本的 np.diagonal -

In [104]: x = np.matrix(np.random.rand(5000,5000))

In [105]: err = np.random.rand(5000)

In [106]: %timeit np.diagonal(x)*err[:,np.newaxis]
10 loops, best of 3: 66.8 ms per loop

In [107]: %timeit x.flat[::x.shape[1]+1].A1 * err[:,None]
10 loops, best of 3: 37.7 ms per loop

用广播替换轴上的循环，第 2 部分

Replace looping-over-axes with broadcasting, pt 2

python

numpy

python-3.x

array-broadcasting