认识广播

广播，英文 broadcasting，有些读者是在使用 NumPy 时，从报错信息中第一次见到 broadcasting。

那么，什么是广播？广播的规则又是怎样的？

在做 10 道 NumPy 练习题前，请先了解这个重要的规则。

在 NumPy 中，下列操作是有效的：

In [1]: import numpy as np

In [2]: v1 = np.arange(10).reshape(2,5) # v1 shape: (2,5)

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

In [4]: v1 + v2
Out[4]:
array([[ 2,  3,  4,  5,  6],
       [ 7,  8,  9, 10, 11]])

v1 的 shape 为 (2,5)，v2 的 shape 为 (1,)，一个为二维，一个为一维。如果没有广播机制，一定会抛(出，维数不等无法相加的异常。

但是，因为广播机制的存在，v2 数组会适配 v1 数组，按照第 0、1 维度，分别发生一次广播，广播后的 v2 变为：

v2 = np.tile(2,(2,5))
Out[10]:
array([[2, 2, 2, 2, 2],
       [2, 2, 2, 2, 2]])

然后，执行再执行加法操作时，因为 v1、v2 的 shape 变得完全一致，所以就能实现相加操作了。

但是，如果 v2 的 shape 为 (2,)，如下，是否 v2 广播后，能实现 v1 + v2 操作？

v2 = np.array([1,2])

执行 v1 + v2 后，抛出 shapes (2,5) (2,) 无法广播到一起的异常。

In [11]: v2 = np.array([1,2])

In [12]: v1 + v2

ValueError: operands could not be broadcast together with shapes (2,5) (2,)

因为 v1、v2 按照广播的规则，无法达成一致的 shape，所以抛出异常。下面了解广播的具体规则。

广播规则

以上看到，不是任意 shape 的多个数组，操作时都能广播到一起，必须满足一定的约束条件。

• NumPy 首先会比较最靠右的维度，如果最靠右的维度相等或其中一个为 1，则认为此维度相等；
• 那么，再继续向左比较，如果一直满足，则认为两者兼容；
• 最后，分别在对应维度上发生广播，以此补齐直到维度一致。

如下，两个数组 a、b，shape 分别为 (2,1,3)、(4,3), 它们能否广播兼容？我们来分析下。

a = np.arange(6).reshape(2,1,3) # shape: (2,1,3)
b = np.arange(12).reshape(4,3) # shape: (4,3)

• 按照规则，从最右侧维度开始比较，数组 a, b 在此维度上的长度都为 3，相等；
• 继续向左比较，a 在此维度上长度为 1，b 长度为 4，根据规则，也认为此维度是兼容的；
• 继续比较，但是数组 b 已到维度终点，停止比较。

结论，数组 a 和 b 兼容，通过广播能实现 shape 一致。

下面看看，数组 a 和 b 广播操作实施的具体步骤。

In [36]: a # 初始 a
Out[36]:
array([[[0, 1, 2]],

       [[3, 4, 5]]])

In [39]: b # 初始 b
Out[39]:
array([[ 0,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8],
       [ 9, 10, 11]])

维度编号从 0 开始，数组 a 在维度 1 上发生广播，复制 4 次：

a = np.repeat(a,4,axis=1)

打印 a ：

Out[35]:
array([[[0, 1, 2],
        [0, 1, 2],
        [0, 1, 2],
        [0, 1, 2]],

       [[3, 4, 5],
        [3, 4, 5],
        [3, 4, 5],
        [3, 4, 5]]])

此时，数组 a 和 b 在后两个维度一致，但是数组 b 维度缺少一维，所以 b 也会广播一次：

b = b[np.newaxis,:,:]  # 首先增加一个维度
b = np.repeat(b,2,axis=0) # 在维度 0 上复制 2 次

经过以上操作，数组 a 和 b 维度都变为 (2,4,3)，至此广播完成，做个加法操作：

In [55]: a + b
Out[55]:
array([[[ 0,  2,  4],
        [ 3,  5,  7],
        [ 6,  8, 10],
        [ 9, 11, 13]],

       [[ 3,  5,  7],
        [ 6,  8, 10],
        [ 9, 11, 13],
        [12, 14, 16]]])

验证我们自己实现的广播操作，是否与 NumPy 中的广播操作一致，直接使用原始的 a 和 b 数组相加，看到与上面得到的结果一致。

a = np.arange(6).reshape(2,1,3) # shape: (2,1,3)
b = np.arange(12).reshape(4,3) # shape: (4,3)
a + b
Out[56]:
array([[[ 0,  2,  4],
        [ 3,  5,  7],
        [ 6,  8, 10],
        [ 9, 11, 13]],

       [[ 3,  5,  7],
        [ 6,  8, 10],
        [ 9, 11, 13],
        [12, 14, 16]]])

至此，广播规则总结完毕。

建议大家都好好理解广播机制，因为接下来使用 NumPy 函数，或者查看文档时，再遇到 broadcast，就知道它的规则，加快对函数的理解。

并且，即便遇到广播不兼容的 Bug 时，相信也能很快解决。

NumPy 练习题

返回有规律的数组

已知数组：

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

返回如下数组：

array([1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3])

分析：

数组前半部分 1, 1, 1, 2, 2, 2, 3, 3, 3 通过 repeat 函数复制 3 次，后面部分通过 tile 函数复制 3 次，体会二者区别。然后合并数据。

In [59]: np.hstack((np.repeat(a,3), np.tile(a,3)))
Out[59]: array([1, 1, 1, 2, 2, 2, 3, 3, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3])

Python 实现向量化

借助 NumPy 的 vectorize 实现操作向量化。

原生的 Python 列表不支持向量化操作，两个列表相加默认不是逐个元素相加：

a = [1,3,5]
b = [2,4,6]
a + b # 默认实现的不是逐个元素相加操作

但是，借助 vectorize 能实现矢量相加：

def add(x,y):
    return x+y
addv = np.vectorize(add)

addv 函数就能实现两个数组相加：

In [67]: addv(a,b)
Out[67]: array([ 3,  7, 11])

限制打印元素的个数

使用 set_printoptions 限制打印元素的个数：

np.set_printoptions(threshold=5)

求中位数

求如下三维数组 a，沿 axis = 1 的中位数。

使用 median 方法，因为 axis 为 1 的数组元素长度为 4，所以中位数为中间两个数的平均数。

如切片 a[0,:,0] 为 [4,8,5,3]，排序后的中间两个元素为 [4,5]，平均值为 4.5。

In [75]: a
array([[[4, 2, 4],
        [8, 2, 7],
        [5, 3, 6],
        [3, 2, 3]],

       [[2, 6, 1],
        [5, 9, 8],
        [9, 7, 1],
        [2, 1, 1]]])
In [73]: ma = np.median(a,axis = 1)
Out[74]:
array([[4.5, 2. , 5. ],
       [3.5, 6.5, 1. ]])

计算 softmax 得分值

已知数组 a，求 softmax 得分值。

In [81]: a
Out[81]:
array([0.07810512, 0.12083313, 0.23554504, 0.62057901, 0.3437597 ,
       0.10876455, 0.08338525, 0.28873765, 0.54033942, 0.71941148])

定义 softmax 函数：

def softmax(a):
    e_a = np.exp(a - np.max(a))
    return e_a / e_a.sum(axis=0)

调用 softmax，得到每个元素的得分，因为 softmax 单调递增函数，所以输入值越大，得分值越高。

In [85]: sm
Out[85]:
array([0.07694574, 0.08030473, 0.0900658 , 0.13236648, 0.10035914,
       0.07934139, 0.0773531 , 0.09498634, 0.12216039, 0.14611689])

sum(sm) 等于 1。

求任意分位数

已知数组 a，求 20 分位数，80 分位数：

a = np.arange(11)

使用 percentile 函数，q 为分位数列表：

In [95]: np.percentile(a,q=[20,80])
Out[95]: array([2., 8.])

找到 NumPy 中缺失值

NumPy 使用 np.nan 标记缺失值，给定如下数组 a，求出缺失值的索引。

如下使用 where 函数，返回满足条件的位置索引：

In [119]: a = np.array([ 0.,  1., np.nan,  3., np.nan, np.nan,  6.,  7.,  8.,  9.])
In [123]: np.where(np.isnan(a))
Out[123]: (array([2, 4, 5], dtype=int64),)   

返回无缺失值的行

给定数组，找出没有任何缺失值的行：

a = np.array([[ 0., np.nan,  2.,  3.],
       [ 4.,  5., np.nan,  7.],
       [ 8.,  9., 10., 11.],
       [12., 13., np.nan, 15.],
       [16., 17., np.nan, 19.],
       [20., 21., 22., 23.]])

求解方法：

In [135]: m = np.sum(np.isnan(a), axis = 1) == 0

In [136]: m
Out[136]: array([False, False,  True, False, False,  True])

In [137]: a[m]
Out[137]:
array([[ 8.,  9., 10., 11.],
       [20., 21., 22., 23.]])

求相关系数

求如下二维数组 a 的相关系数：

a = array([[ 2, 12, 21, 10],
       [ 1, 20,  8, 22],
       [ 7,  1,  5,  1],
       [ 7, 10, 14, 14],
       [12, 13, 13, 14],
       [ 0, 12, 21,  2]])

如下使用 corrcoef 方法，求得两列的相关系数，相关系数为 0.242：

In [149]: np.corrcoef(a[:,1],a[:,2])
Out[149]:
array([[1.        , 0.24230838],
       [0.24230838, 1.        ]])

缺失值默认填充为 0

如下数组，含有缺失值，使用 0 填充：

a = np.array([[ 0., np.nan,  2.,  3.],
       [ 4.,  5., np.nan,  7.],
       [ 8.,  9., 10., 11.],
       [12., 13., np.nan, 15.],
       [16., 17., np.nan, 19.],
       [20., 21., 22., 23.]])

一行代码，np.isnan(a) 逐元素检查，若为空则为 True，否则为 False，得到一个与原来 shape 相同的值为 True 和 False 的数组。

In [157]: a[np.isnan(a)] = 0

In [158]: a
Out[158]:
array([[ 0.,  0.,  2.,  3.],
       [ 4.,  5.,  0.,  7.],
       [ 8.,  9., 10., 11.],
       [12., 13.,  0., 15.],
       [16., 17.,  0., 19.],
       [20., 21., 22., 23.]])

使用 NumPy 处理 fashion-mnist 数据集

fashion-mnist 是一个与手写字一样经典的数据集，与服饰相关。

导入数据特征数 785，我们先提取 0 到 784，然后 reshape 为 28*28 的二维数组。

首先导入 NumPy，截取前 784 个元素，reshape 为 28*28 的元素。大家自行下载此数据集，下载并导入后，train_data 的 shape 为 (*,785)。

import numpy as np
train_data = fashion_mnist_train.to_numpy() # Pandas DataFrame 转 numpy 对象
row0 = train_data[0,:784].reshape(28,-1) 

导入 Matplotlib，使用 imshow 绘制 784 个像素（取值为 0~255）。

import matplotlib.pyplot as plt
plt.imshow(row0)
plt.show()

展示的图像，如下所示：

依次展示前 10 幅图：

for i in range(10):
    print('图%d'%(i,))
    plt.imshow(train_data[i,:784].reshape(28,-1))
    plt.colorbar()
    plt.show()

小结

今天学习了NumPy 中最重要的机制之一：广播，解答了广播是什么，规则又是什么，推荐大家理解透。

后面介绍 10 道 NumPy 练习题，进一步巩固 NumPy 中常用的函数。

最后使用一个实际数据集，图形背后也是数值型的多维数组，NumPy 也会发挥重大作用。