多个阵列上的 Numpy 广播

Question

我有一个平面的 3 个维度的基础：(u, v)。

我想获得这个基础的所有线性组合，以基本遍历我的整个平面：

对于 [0, 512[ 中的 i 和 [0, 512[ 中的 j，获取所有 (i * u + j * v)。

我需要这个速度很快，所以 for 循环并不是一个真正的选择。如何使用 numpy 广播做到这一点？

看了https://docs.scipy.org/doc/numpy-1.13.0/user/basics.broadcasting.html，我的印象是不可能做到...

试过了：

# This is an orthonormal basis but there is no guarantee it is
u = np.array([1, 0, 0])
v = np.array([0, 1, 0])
tmp = np.arange(512)
factors = itertools.combinations(tmp, 2)
pixels = factors[0] * u + factors[1] + v

但显然它不起作用。

这个问题有解决办法吗？如果是，那怎么办？

Answer 1

将 (u, v) 与 2D 索引网格相乘：

ind = np.indices((512, 512))
pixels = ind[0, ..., np.newaxis] * u + ind[1, ..., np.newaxis] * v

>>> %timeit ind = np.indices((512, 512)); pixels = ind[0, ..., np.newaxis] * u + ind[1, ..., np.newaxis] * v
8.06 ms ± 69.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

将 u 乘以 1D 索引范围，将 v 乘以 1D 索引范围，广播并合并为 2D：

i512 = np.arange(512)[:, np.newaxis]
pixels = (i512 * u)[:, np.newaxis, :] + (i512 * v)[np.newaxis, :, :]

>>> %timeit i512 = np.arange(512)[:, np.newaxis]; pixels = (i512 * u)[:, np.newaxis, :] + (i512 * v)[np.newaxis, :, :]
4.06 ms ± 58.6 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Answer 2

对您的代码进行一些修改

factors = itertools.product(tmp, tmp)
factors = list(zip(*factors))
factors = np.array(factors)
pixels = factors[0][:, None] * u + factors[1][:, None] + v

首先，一个数学错误：您想要笛卡尔积，而不是组合。

现在，实际的语法错误：itertools 生成 [(i1, j1), (i2, j2)..] 列表，而不是 numpy 数组。所以，为了让你的最后一行起作用，你必须

1. zip(*) 将列表转换为(i1, i2), (j1, j2) 格式
2. 使它成为一个 numpy 数组
3. 在1D的factors[0], factors[1]向量中，添加[:, None]以将其转换为列，以便广播工作。