python - NumPy 中序列的乘积

Question

我需要用 NumPy 实现以下功能 -

其中 F_l(x) 是我需要计算的 N 个数组，它们依赖于数组 G(x) ，我得到的，A_j 是也给出的 N 系数。我想在 NumPy 中实现它，因为我必须为程序的每次迭代计算 F_l(x)。执行此操作的虚拟方法是使用 for 循环和 ifs:

import numpy as np
A = np.arange(1.,5.,1)
G = np.array([[1.,2.],[3.,4.]])

def calcF(G,A):
    N = A.size
    print A
    print N
    F = []
    for l in range(N):
        F.append(G/A[l])
        print F[l]
        for j in range(N):
            if j != l:
                F[l]*=((G - A[l])/(G + A[j]))*((A[l] - A[j])/(A[l] + A[j]))
    return F

F= calcF(G,A)
print F

至于循环和 if 语句相对较慢，我正在寻找一种 NumPy 机智的方法来做同样的事情。有人有想法吗？

Answer 1

本文中列出的是大量使用 NumPy's powerful broadcasting feature 的矢量化解决方案使用 np.newaxis/None 将输入数组的维度扩展到 3D 和 4D 情况后根据所涉及的计算在不同的地方。这是实现-

# Get size of A
N = A.size

# Perform "(G - A[l])/(G + A[j]))" in a vectorized manner
p1 = (G - A[:,None,None,None])/(G + A[:,None,None])

# Perform "((A[l] - A[j])/(A[l] + A[j]))" in a vectorized manner
p2 = ((A[:,None] - A)/(A[:,None] + A))

# Elementwise multiplications between the previously calculated parts
p3 = p1*p2[...,None,None]

# Set the escaped portion "j != l" output as "G/A[l]"
p3[np.eye(N,dtype=bool)] = G/A[:,None,None]
Fout = p3.prod(1)

# If you need separate arrays just like in the question, split it
Fout_split = np.array_split(Fout,N)

sample 运行-

In [284]: # Original inputs
     ...: A = np.arange(1.,5.,1)
     ...: G = np.array([[1.,2.],[3.,4.]])
     ...: 

In [285]: calcF(G,A)
Out[285]: 
[array([[-0.        , -0.00166667],
        [-0.01142857, -0.03214286]]), array([[-0.00027778,  0.        ],
        [ 0.00019841,  0.00126984]]), array([[  1.26984127e-03,   1.32275132e-04],
        [ -0.00000000e+00,  -7.93650794e-05]]), array([[-0.00803571, -0.00190476],
        [-0.00017857,  0.        ]])]

In [286]: vectorized_calcF(G,A) # Posted solution 
Out[286]: 
[array([[[-0.        , -0.00166667],
         [-0.01142857, -0.03214286]]]), array([[[-0.00027778,  0.        ],
         [ 0.00019841,  0.00126984]]]), array([[[  1.26984127e-03,   1.32275132e-04],
         [ -0.00000000e+00,  -7.93650794e-05]]]), array([[[-0.00803571, -0.00190476],
         [-0.00017857,  0.        ]]])]

运行时测试-

In [289]: # Larger inputs
     ...: A = np.random.randint(1,500,(400))
     ...: G = np.random.randint(1,400,(20,20))
     ...: 

In [290]: %timeit calcF(G,A)
1 loops, best of 3: 4.46 s per loop

In [291]: %timeit vectorized_calcF(G,A)  # Posted solution 
1 loops, best of 3: 1.87 s per loop

---

使用 NumPy/MATLAB 进行矢量化:一般方法

感觉我可以在我的一般方法上投入两分钱，我认为其他人在尝试向量化代码时会遵循类似的策略，尤其是在像 NumPy 或 MATLAB 这样的高级平台上。所以，这里有一份快速 list ，列出了可以考虑用于 Vectorization 的事项。 -

关于扩展维度的想法:要为输入数组扩展维度，以便新维度包含本应在嵌套循环中迭代生成的结果。

从哪里开始向量化？从计算的最深层(代码迭代次数最多的循环)阶段开始，看看如何扩展输入并引入相关计算。仔细跟踪涉及的迭代器并相应地扩展维度。向外移动到外部循环，直到您对完成的矢量化感到满意。

如何处理条件语句？对于简单的情况，暴力计算所有内容，然后查看如何处理 IF/ELSE 部分。这将是高度特定于上下文的。

是否存在依赖关系？如果存在，请查看是否可以相应地跟踪和实现依赖关系。这可能会形成另一个讨论主题，但这里有 few examples我参与其中。