不同大小向量的 Numpy 乘法避免循环

Question

我有一个矩阵，比方说，P 大小 (X,Y)。另外，我有两个矩阵，例如 Kx 和 Ky 大小 (M,N)，一个矩阵 pk 大小 (M,N) 和两个向量 u 和 v 分别是 X 和 Y。例如，它们可以定义如下：

import numpy as np
P = np.zeros((X,Y));
pk = np.random.rand(M,N);
Kx = np.random.rand(M,N);
Ky = np.random.rand(M,N);
u = np.random.rand(X);
v = np.random.rand(Y);

当然，在实际代码中它们不是随机的，但这对于本例来说无关紧要。问题是，是否存在等同于以下内容的纯 numpy：

for m in range(0, M):
    for n in range(0, N):
        for i in range(0,X):
            for j in range(0,Y):
               Arg = Kx[m,n]*u[i] + Ky[m,n]*v[j];
               P[i,j] += pk[m,n]*np.cos(Arg);

所有M、N、X、Y都不同，但是X和Y可以相同，如果否则不存在解决方案。

Answer 1

在 NumPy 计算中消除 for-loop 的常见策略是使用高维数组。

例如考虑以下行

Arg = Kx[m,n]*u[i] + Ky[m,n]*v[j]

此行取决于索引 m、n、i 和 j。 因此 Arg 取决于 m、n、i 和 j。这意味着 Arg 可以被认为是一个由 m、n、i 和 j 索引的 4 维数组。所以我们可以消除 4 for-loops——就 Arg 而言——通过计算

Kxu = Kx[:,:,np.newaxis]*u
Kyv = Ky[:,:,np.newaxis]*v
Arg = Kxu[:,:,:,np.newaxis] + Kyv[:,:,np.newaxis,:]

Kx[:,:,np.newaxis] 的形状为 (M, N, 1)，u 的形状为 (X,)。将它们相乘使用 NumPy broadcasting to create an array of shape (M, N, X). Thus, above, new axes are used somewhat like placeholders，因此 Arg 最终得到按 m、n、i、j 索引的 4 个轴.

同理，P可以定义为

P = (pk[:,:,np.newaxis,np.newaxis]*np.cos(Arg)).sum(axis=0).sum(axis=0)

sum(axis=0)（调用两次）沿 m 和 n 轴求和，因此 P 最终成为由 [= 索引的二维数组仅 19=] 和 j。

通过使用这些 4 维数组，我们可以对整个 NumPy 数组应用 NumPy 操作。相反，当使用 4 for-loops 时，我们必须对标量进行逐值计算。例如，当 Arg 是一个 4 维数组时，考虑 np.cos(Arg) 正在做什么。这在一个 NumPy 函数调用中卸载了 all 余弦的计算，该函数调用在已编译的 C 代码中执行底层循环。这比为每个标量调用一次 np.cos 快得多。这就是为什么使用高维数组最终比基于 for-loop 的代码快得多的原因。

import numpy as np

def orig(Kx, Ky, u, v, pk):
    M, N = Kx.shape
    X = u.size
    Y = v.size
    P = np.empty((X, Y), dtype=pk.dtype)
    for m in range(0, M):
        for n in range(0, N):
            for i in range(0,X):
                for j in range(0,Y):
                   Arg = Kx[m,n]*u[i] + Ky[m,n]*v[j]
                   P[i,j] += pk[m,n]*np.cos(Arg)
    return P

def alt(Kx, Ky, u, v, pk):
    Kxu = Kx[:,:,np.newaxis]*u
    Kyv = Ky[:,:,np.newaxis]*v
    Arg = Kxu[:,:,:,np.newaxis] + Kyv[:,:,np.newaxis,:]
    P = (pk[:,:,np.newaxis,np.newaxis]*np.cos(Arg)).sum(axis=0).sum(axis=0)
    return P

M, N = 10, 20
X, Y = 5, 15
Kx = np.random.random((M, N))
Ky = np.random.random((M, N))
u = np.random.random(X)
v = np.random.random(Y)
pk = np.random.random((M, N))

完整性检查，（显示 alt 和 orig return 相同的结果）：

In [57]: P2 = alt(Kx, Ky, u, v, pk)

In [58]: P1 = orig(Kx, Ky, u, v, pk)

In [59]: np.allclose(P1, P2)
Out[59]: True

一个基准测试，显示 alt 明显快于 orig：

In [60]: %timeit orig(Kx, Ky, u, v, pk)
10 loops, best of 3: 33.6 ms per loop

In [61]: %timeit alt(Kx, Ky, u, v, pk)
1000 loops, best of 3: 349 µs per loop

不同大小向量的 Numpy 乘法避免循环

Numpy multiplication of vectors of different size avoiding for loops

python

numpy

matrix

matrix-multiplication