Numpy多维数组中索引的顺序

Question

例如，假设我正在模拟一堆粒子随时间做某些事情，我有一个名为particles的多维数组，其中包含这些索引：

• 粒子的x/y/z坐标（长度为a，对于3d空间来说是3）
• 单个粒子的索引（长度为b）
• 它所在的时间步骤的索引（长度为c）

构建数组时，particles.shape == (a, b, c)或particles.shape == (c, b, a)哪种更好？

我更关心惯例而不是效率：Numpy数组可以按C风格（最后一个索引变化最快）或Fortran风格（第一个索引）设置，因此它可以高效地支持任何一种设置。我也意识到我可以使用transpose将索引以任何顺序排列，但我想将其最小化。

我开始自己研究并找到了支持两种方式的论据：

赞成(c,b,a)：

• 默认情况下，Numpy使用C风格的数组，其中最后一个索引是变化最快的。
• 大多数向量代数函数（inner、cross等）作用于最后一个索引。（dot在一个数组的最后一个索引和另一个数组的倒数第二个索引之间进行计算。）
• matplotlib集合对象（LineCollection、PolyCollection）期望在最后一个轴中具有空间坐标的数组。

Pro-(a,b,c):

• 如果我想使用meshgrid和mgrid来生成一组点，它会将空间轴放在第一位。例如，np.mgrid[0:5,0:5,0:5].shape == (3,5,5,5)。我知道这些函数主要是用于整数数组索引，但通常也可以用它们来生成点的网格。
• matplotlib的scatter和plot函数将其参数分离，因此它对数组的形状不感知，但是ax.plot3d(particles[0], particles[1], particles[2])比带有particles[..., 0]的版本更短。

总的来说，似乎存在两种不同的约定（可能是由于C和Fortran之间的历史差异造成的），并且不清楚哪种约定在Numpy社区中更常见，或者更适合我的工作。

- Apocheir

Answer 1

根据我的经验，这类问题的常规做法更多地与特定的文件格式有关，而不是其他方面。然而，有一种快速的方法可以回答你正在做什么最好的问题：

如果你需要迭代一个轴，你最有可能迭代哪一个轴？换句话说，这些中哪一个最有可能：

# a first
for dimension in particles:
    ...

# b first
for particle in particles:
    ...

# c first
for timestep in particles:
    ...

就效率而言，这假设是C顺序，但这实际上在这里是无关紧要的。在Python层面上，无论内存布局如何，访问NumPy数组都被视为C顺序。(您总是遍历第一个轴，即使它不是内存中“最连续”的轴。)

当然，在许多情况下，应避免以这种方式直接迭代NumPy数组。尽管如此，这是您应该考虑的方式，特别是涉及到磁盘文件结构时。使您最常见的用例成为最快/最容易的用例。

如果没有其他内容，希望这能给您提供有用的思考方式。

- Joe Kington

Answer 2

另一个偏见是，在需要添加新维度时，numpy 的偏好是在左侧添加。也就是说，x[None,...] 是自动的。

np.array([x,y,z])   # produces a (3,...) array

np.ones((3,2)) + np.ones((1,2,10)) # error
np.ones((3,2,1)) + np.ones((2,10))  # (3,2,10)

但我不明白这种先前广播如何有利于 x/y/z 坐标的任何一种位置。

虽然 np.dot 使用最后两个轴的约定，但 np.tensordot 和 np.einsum 更加通用。

---

Apocheir 指出，在最后一个轴上进行缩减可能需要添加一个 newaxis。

 x / np.linalg.norm(x,axis=0)   # automatic newaxis at beginning
 x / np.linalg.norm(x,axis=-1)[...,np.newaxis]  # explicit newaxis

对于小的x，这个明确的newaxis会增加可测量的执行时间。但对于大的x，第二个计算更快。我认为这是因为在最后一个轴上进行缩减更快——这是改变速度更快的轴（对于order='C'）。许多内置的缩减方法都有一个keepdims参数，以便在像这样的用途中进行广播（例如sum、mean）。

- hpaulj

---

实际上，它确实有影响...将x = np.ones((3,4,5));y = np.linalg.norm(x,axis=0)与x = np.ones((5,4,3));y = np.linalg.norm(x,axis=-1)进行比较。当空间索引在第一位时，x/y可以在没有任何索引调整的情况下对x进行标准化。 当空间索引在最后时，您必须执行类似于x/y[..., np.newaxis]的操作。 - Apocheir

许多“缩减”方法（例如sum，mean）都有一个keepdims参数，以消除需要将此“newaxis”添加回去的需要。 - hpaulj