数组的拷贝

Python 本身有深拷贝和浅拷贝等复杂的概念，对于多层列表等数据结构，对应的行为差异非常明显。
但是 Numpy 的 ndarray 数组无论形状如何，在内存中始终是连续存储的，高维数组并不是数组的数组，
因此并没有深拷贝和浅拷贝的明显差异。

和其它 Python 数据一样，在函数传参过程中对 ndarray 数组的传递并不会触发任何的拷贝行为。

为了获取一个 ndarray 数组的完整拷贝，可以使用 ndarray.copy() 方法，例如

1
2
3
4
5
a = np.array([1, 2, 3])
b = a.copy()

b[0] = 100
a  # array([1, 2, 3])

上述测试代码表明，拷贝得到的是完全独立的数组，修改不会相互影响。

也可以使用numpy.copy()函数，例如

1
2
a = np.array([1, 2, 3])
b = np.copy(a)

ndarray.copy() 方法和 numpy.copy() 函数的行为是基本一致的，都是获取原数组的拷贝，但是两者关于内存排布参数的默认值有细微差异，可以参考具体文档。

数组的视图

为了优化计算效率，避免不必要的拷贝，Numpy 提供了视图机制。

一个数组对象其实只需要对外提供访问任意位置元素的接口，为了实现这个需求，在内部只需要持有一个指向内存中实际数据的指针，数组的元素类型，数组的形状和一些偏移量参数即可。
因此两个数组对象实际上完全可以共享内存数据，但是允许两者的指针有一定的初始偏移，并且使用完全不同的形状，这就是视图机制的本质。
这两个数组通常不是地位对等的，习惯上将后创建的数组称为前者的视图。

视图虽然节约了很多不必要的拷贝操作，但是由于它们实际共享的是同一块内存，只是解读方式不同，因此对其进行的修改会相互影响，这是尤其需要注意的。

最基础的创建视图的做法是使用 ndarray.view() 方法，例如

1
2
3
4
5
a = np.array([1, 2, 3])
b = a.view()

b[0] = 100
a  # array([100, 2, 3])

上述测试代码表明，对视图的修改会影响原始数组。在实践中很少直接使用 ndarray.view() 方法，很多具有实际意义的操作都会获得原始数组的视图。

关于视图和拷贝，有两类检查手段，第一种是检查共享内存是否重叠：使用 np.shares_memory 函数判断两个数组是否共享内存，例如

1
2
3
4
5
6
a = np.array([1, 2, 3])
b = a.copy()
c = a.view()

np.shares_memory(a,b) # False
np.shares_memory(a,c) # True

第二种是通过ndarray.base属性检查数组是否是派生的视图：

  • 如果base属性为空，那么表示该数组是独立的，可能来自于直接创建或拷贝；
  • 如果base属性不为空，那么表示该数组是视图，base 属性指向源数组（但是如果存在多次派生，只会指向最近的源数组）。
1
2
3
4
5
6
a = np.array([1, 2, 3])
b = a.copy()
c = a.view()

b.base is None # True
c.base is a # True

注意下面这种写法：

1
2
a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]).view()
a.base is None # False

此时我们直接获取了一个新的 ndarray 数组的视图a，此时a.base属性非空，但是实际上我们只有一种方式去访问对应的底层数据，因此即使是视图，也不需要考虑修改的相互干扰问题。

Python & Numpy 涉及这部分内容的底层机制非常复杂，很多函数可能返回的是视图或拷贝（如果视图不可用），判断也不是一定可靠的，这些都是 Python 的底层实现所决定的。与之不同的是，MATLAB 采用了深拷贝+懒拷贝的机制，既保证了运算效率，也降低了编程时的心智困难。

几个基本的原则：

  • 为了提供效率，很多针对数组的操作的语义都是返回视图；
  • 为了确保获得拷贝，可以使用 ndarray.copy() 方法。
  • 如果语义是获得视图，那么就会尽量返回视图，在尝试构造视图失败时也会返回拷贝；
  • 如果语义是获得拷贝，那么保证不可能返回视图。

索引获得视图？拷贝？

Numpy 对数组的各种索引操作可以用来提取数组的部分元素，此时可以直接对其赋值，会对数组进行修改，
也可以将其赋值给一个新的变量，此时存在一个问题：赋值产生的新变量是原数组的视图还是拷贝？

通常情况下，使用整数或切片索引进行赋值会创建视图，例如

1
2
3
4
5
6
7
a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

b = a[0]
np.shares_memory(a, b)  # True

c = a[0, :]
np.shares_memory(a, c)  # True

使用高级索引进行赋值则会创建新数组，例如

1
2
3
4
5
6
7
a = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

b = a[[0]]
np.shares_memory(a, b)  # False

c = a[a>3]
np.shares_memory(a, c)  # False

这只是最基础的两类情况，实际还存在很多复杂的情况，例如整数或切片索引与高级索引混合使用，这里不作讨论。

数组的形状变换（一）

下面介绍一下常用的 ndarray 数组的方法，这些方法通常都有对应且同名的函数版本，两者的语义基本相同，但是可能在参数处理上存在细微区别，我们主要关注方法的行为。

reshape

一个最常见的需求是修改数组的形状，使用 ndarray.reshape() 方法可以返回修改后形状的视图，例如

1
2
3
4
a = np.array([[1, 2, 3], [4, 5, 6]])
b = a.reshape((3,2))

np.shares_memory(a,b) # True

在指定新的形状时可以使用一个 -1 占位，Numpy 在可行的情况下会自动推导这个维度对应的尺寸。

1
2
3
4
a = np.array([[1, 2, 3], [4, 5, 6]])
b = a.reshape((-1,1)) # shape = (6,1)

np.shares_memory(a,b) # True

为了便于使用，ndarray.reshape() 不仅接收一个完整的形状参数，还支持依次传递形状参数（numpy.reshape() 不支持），例如

1
a.reshape((2, 3)) == a.reshape(2, 3) # ok

虽然某些情况下直接修改 ndarray.shape 是可行的，Numpy 为这个属性提供了 setter，但是这并不是推荐的做法，在后续版本可能被废弃。

resize

与 ndarray.reshape() 不同，使用 ndarray.resize() 方法可以就地修改数组形状，没有返回值。
由于新的形状可以改变总元素个数（截断或延拓），因此不允许使用 -1 占位。

例如

1
2
3
4
5
6
a = np.array([[1, 2, 3], [4, 5, 6]])
a.resize((1,4))
# array([[1, 2, 3, 4]])

a.resize((1,6))
# array([[1, 2, 3, 4, 0, 0]])

通过实验可知，如果 resize 使得总长度发生改变，对应的行为就是直接丢弃最后的元素，或者使用 0 填充作为新的元素。

flatten / ravel

考虑这样的需求：将高维数组按照实际内存顺序展平为一维数组进行访问，此时有两类具体的行为：返回原数组的视图或拷贝。

使用 ndarray.flatten() 可以获取原数组展平为一维数组所得到的拷贝，例如

1
2
3
4
a = np.array([[1, 2], [3, 4]])
b = a.flatten() # array([1, 2, 3, 4])

np.shares_memory(a,b) # False

使用 ndarray.ravel() 则可以获取原数组展平为一维数组所得到的视图，例如

1
2
3
4
a = np.array([[1, 2], [3, 4]])
b = a.ravel() # array([1, 2, 3, 4])

np.shares_memory(a,b) # True

如果需要返回展平后的一维数组视图，官方文档更推荐的做法其实是 a.reshape(-1)，因为常用方法的可读性更高。

transpose

ndarray.transpose() 方法可以返回二维数组转置后的视图，例如

1
2
3
4
5
6
7
a = np.array([[1, 2, 3], [4, 5, 6]])
b = a.transpose()
# array([[1, 4],
#        [2, 5],
#        [3, 6]])

np.shares_memory(a,b) # True

由于转置操作非常普遍，Numpy 直接提供了一个 ndarray.T 属性，相当于调用 ndarray.transpose() 方法

1
a.T == a.transpose()

转置操作对于非二维数组也是合法的，但是效果可能并不是所期望的，默认情况下，它只是将所有的维度翻转，也就是之前的第一个轴翻转为最后一个，第二个轴翻转为倒数第二个，以此类推。
例如 shape = (2,3,4) 的数组，使用 ndarray.transpose() 方法会返回 shape = (4,3,2) 的数组。

1
2
a = np.zeros((2,3,4))
a.transpose().shape # (4, 3, 2)

按照这个规则，对于一维数组，由于 shape = (n,)，调用 ndarray.transpose() 方法仍然返回原数组的视图，而不是通常所期望的 shape = (n, 1) 的列向量，这一点需要特别注意。

1
2
a = np.array([1,2,3])
a.transpose() # array([1, 2, 3])

当然，对于多维数组，我们可以传递参数以指定操作前后轴的对应关系，例如

1
2
3
4
5
6
a = np.zeros((2, 3, 4))
a.transpose().shape  # (4, 3, 2)

a.transpose((2, 1, 0)).shape  # (4, 3, 2)
a.transpose((1, 2, 0)).shape  # (3, 4, 2)
a.transpose((2, 0, 1)).shape  # (4, 2, 3)

swapaxes

除了通过 ndarray.transpose() 方法一次性调整所有轴，
还可以使用 ndarray.swapaxes() 方法对指定的两个轴进行交换，同样也是返回视图，例如

1
2
3
4
5
6
a = np.zeros((2, 3, 4))
b = a.swapaxes(0, 1)
np.shares_memory(a, b)  # True

a.swapaxes(0, 1).shape  # (3, 2, 4)
a.swapaxes(1, 2).shape  # (2, 4, 3)

squeeze

有时我们需要处理类似 (1,3,1) 这种“尺寸冗余”的数据，可以使用 ndarray.squeeze() 方法将其压缩为一维数组，可以简化后续对数组的操作，注意返回的仍然是视图。

1
2
3
4
5
6
7
8
9
a = np.array([[[1], [2], [3]]])  # shape: (1, 3, 1)
# array([[[1],
#         [2],
#         [3]]])

b = a.squeeze()  # shape: (3,)
# array([1, 2, 3])

np.shares_memory(a, b)  # True

也可以指定需要压缩的轴，此时其它轴即使长度为1也会保持，但是必须保证指定的轴维数为1，否则会抛出错误。

1
2
3
a = np.array([[[1], [2], [3]]])  # shape: (1, 3, 1)

c = a.squeeze(axis=2) # shape: (1, 3)

数组的形状变换（二）

下面这几个是 Numpy 函数或特殊对象，并没有对应的 ndarray 方法。

np.expand_dims

作为 squeeze 的逆操作，有时我们需要增加数组的轴来匹配某些运算要求，可以使用 numpy.expand_dims() 函数（注意没有对应的方法），，注意返回的仍然是视图。

1
2
3
4
5
6
a = np.array([1, 2, 3])  # shape: (3,)

b = np.expand_dims(a, axis=0)  # shape: (1, 3)
np.shares_memory(a, b)  # True

c = np.expand_dims(a, axis=(0, 1))  # shape: (1, 1, 3)

通过axis参数指定额外的长度为1的轴的位置，注意这个参数必须是合理的，否则会报错。

np.newaxis

我们可以在索引操作中直接增加数组的轴，这需要用到一个特殊的 np.newaxis 对象来占位，表示在对应位置增加一个轴，例如

1
2
3
4
5
a = np.array([1, 2, 3])  # shape: (3,)

a[:, np.newaxis].shape  # (3, 1)
a[np.newaxis, :].shape  # (1, 3)
a[np.newaxis, :, np.newaxis].shape  # (1, 3, 1)

np.newaxis 在本质上只是 None 的别名，在一些代码中会直接使用 None。

1
np.newaxis is None # True

np.atleast_Xd

与 np.expand_dims() 函数类似，np.atleast_Xd() 函数也是增加数组的轴（其中 X=1,2,3），可以用于保证数组的维度不小于 X，返回满足要求的视图：

  • 对于满足要求的数组，不进行任何操作；
  • 对于不满足要求的数组，会在末尾增加若干个长度为1的轴，使其达到要求。

例如

1
2
3
4
a = np.array([2, 3])

b = np.atleast_2d(a)  # shape: (1, 2)
np.shares_memory(a, b)  # True

显然，np.atleast_1d() 函数只在输入标量时有意义，此时会将其转换为 ndarray 数组

1
2
3
a = 1
b = np.atleast_1d(a)  # shape: (1,)
# array([1])

这对于实现兼容标量和 ndarray 数组作为参数的函数很有用。

np.atleast_2d() 函数只在输入标量和一维数组时有意义

1
2
3
a = np.array([1,2,3])
b = np.atleast_2d(a) # shape: (3,1)
# array([[1, 2, 3]])

在使用 np.atleast_1d() 的同时，可以搭配使用 np.isscalar() 判断输入是否是标量，
从而对于不同的输入类型提供不同的输出类型（如果输入标量，仍然返回标量）

1
2
3
4
5
np.isscalar(1) # True
np.isscalar(2.0) # True

np.isscalar(np.array([1])) # False
np.isscalar([2.0]) # False

数组的拼接

有时我们需要将两个数组拼接成一个大数组，此时显然不可能返回一个视图，必然是一个独立的新数组，
这里关注最通用的 np.concatenate() 函数，它可以将若干个数组沿着已经存在的某个轴进行拼接。

这个函数在使用时提供的参数通常有两个：

  • 第一个是位置参数，需要提供待合并的多个数组组成的元组或列表；
  • 第二个则是关键字参数axis，用于指定合并的轴。（默认值axis=0）

函数的语义为：

  • 如果 axis 被指定为非负整数，则沿着 axis 对应的轴进行拼接，要求所有数组除了 axis 轴之外的其它轴的维度必须一致，拼接得到的数组的 shape 对于其它轴保持不变，对于执行拼接的轴，维度为所有数组这个轴的维数之和；
  • 如果 axis=None，则将所有数组展平为一维数组，然后顺序拼接。

例如

1
2
3
4
5
a = np.array([[1, 2], [3, 4]])  # shape (2, 2)
b = np.array([[5, 6]])  # shape (1, 2)

np.concatenate([a, b], axis=0)  # shape (3, 2)
# np.concatenate([a, b], axis=1)  # error

Numpy 实际上根据各种具体情景还提供了 np.stack(), np.hstack(), np.vstack() 等函数，但是适用情景比较特殊，而且操作可以被 np.concatenate() 函数替代。

数组的堆叠

numpy.tile() 函数可以将一个已知的小数组重复堆叠来构造大矩阵，例如

1
2
3
4
5
6
7
A = np.array([[1, 2], [3, 4]])

B = np.tile(A, (2, 3))
# array([[1, 2, 1, 2, 1, 2],
#        [3, 4, 3, 4, 3, 4],
#        [1, 2, 1, 2, 1, 2],
#        [3, 4, 3, 4, 3, 4]])

在理想情况下，堆叠参数指定了数组每一个维度的重复次数，例如形状为 (m1, ..., mk) 的数组，堆叠次数为 (r1, ..., rk)，那么最终的结果是形状为 (m1*r1, ..., mk*rk) 的数组。
例如

1
2
A = np.zeros((2, 3, 4))
np.tile(A, (5, 6, 7)).shape # (10, 18, 28)

在两者的轴数不匹配的情况下：若数组的轴数过少，会在前面加 1；若堆叠参数的轴数过少，同样也在前面加 1，以保证两者匹配。
例如

1
2
3
4
5
A = np.zeros((2, 3))
np.tile(A, (4, 5, 6)).shape  # (4, 10, 18)

A = np.zeros((2, 3, 4))
np.tile(A, (5, 6)).shape  # (2, 15, 24)

其它

重复索引更新问题

考虑下面一段代码

1
2
3
4
5
6
b = np.array([0.0])
idx = np.array([0, 0, 0])

b[idx] += np.array([1.0, 2.0, 3.0])

print(b)

按照直觉来说，这段代码会对 b[0] 进行三次就地更新：

1
2
3
b[0] += 1
b[0] += 2
b[0] += 3

最终结果应该是 b[0]=6.0。但是实际执行的结果是 b[0]=3.0。

这是因为在这段代码中使用的 NumPy 高级索引会生成拷贝的临时数组，而不是视图。

b[idx] 实际等价于

1
tmp = np.array([b[0], b[0], b[0]])

b[idx] += v 的实际过程为：对生成的临时数组进行 +=

1
2
3
tmp = np.array([b[0], b[0], b[0]])
tmp += np.array([1.0, 2.0, 3.0])
# tmp: [1.0, 2.0, 3.0]

然后写回原数组

1
2
3
b[0] = 1
b[0] = 2
b[0] = 3

由于重复索引会被覆盖写回，最终结果为 b[0]=3.0。

numpy 提供了 np.add.at 函数来解决这类需求，具体细节不再展开，编程时注意重复索引可能出现问题即可。事实上，重复索引即使在语义明确的情况下，也会因为并行计算中的数据竞争而出现结果随机的问题。

np.where

numpy 数组支持对 array 数组进行条件判断和条件赋值，都通过numpy.where()实现：

  • 条件判断：np.where(condition)，返回满足条件的索引分量元组，只是np.asarray(condition).nonzero()的快捷方式，注意与之不同的是，直接 condition 返回的是布尔数组。
  • 条件赋值：np.where(condition, x, y)，在满足条件的位置使用数组 x 赋值，不满足条件的位置使用数组 y 赋值，也就是对每一个位置应用三元表达式

条件判断的用法例如

1
2
3
4
a = np.array([10, 20, 30, 40])

np.where(a > 15)  # (array([1, 2, 3]),)
np.nonzero(a > 15)  # (array([1, 2, 3]),)

条件赋值则更加常用，例如

1
2
3
4
a = np.array([10, 20, 30, 40])
b = np.array([20, 20, 20, 20])

np.where(a > 20, a, b) # array([20, 20, 30, 40])

在操作中会自动使用广播机制，因此直接提供一个标量也可以

1
2
a = np.array([10, 20, 30, 40])
np.where(a > 20, 20, a)  # array([10, 20, 20, 20])

遍历数组

array 数组可以直接放在 for 语句中进行遍历，此时的行为是按照第一个轴进行遍历。
还可以使用 ndarray.flat 属性获取一个对所有元素遍历的迭代器，例如

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
a = np.array([[1,2,3],[4,5,6],[7,8,9]])

for item in a:
    print(item)

# [1 2 3]
# [4 5 6]
# [7 8 9]

for item in a.flat:
    print(item)

# 1
# 2
# 3
# 4
# 5
# 6
# 7
# 8
# 9

保存和加载数组

Numpy 直接提供了 .npy 和 .npz 格式用于保存和加载 array 数组。
这是一种 NumPy 专用的二进制格式，速度快、精度不变，推荐用于保存和加载中间数据。

可以使用 np.save() 和 np.load() 来保存和加载 array 数组。

1
2
3
4
5
a = np.array([[1, 2, 3], [4, 5, 6]])
np.save("array1", a)

a_loaded = np.load("array1.npy")
print(a_loaded)

其中 np.save() 函数对于文件名可能会自动添加 .npy 后缀，但是 np.load() 函数不会补全后缀。

对于多个数组，可以使用 np.savez() 将其保存为 .npz 文件。

1
2
3
4
5
6
7
8
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

np.savez("data.npz", array1=a, array2=b)

data = np.load("data.npz")
print(data["array1"])  # [1 2 3]
print(data["array2"])  # [4 5 6]