我们在前面的章节中学习了NumPy的通用函数,它们用来对数组进行 向量化 操作,从而抛弃了性能低下的Python循环。 还有一种对NumPy数组进行向量化操作的方式我们称为 广播 。 广播简单来说就是一整套用于在不同尺寸或形状的数组之间进行二元 ufuncs 运算(如加法、减法、乘法等)的规则。
import numpy as np
a = np.array([0, 1, 2])
b = np.array([5, 5, 5])
a + b
array([5, 6, 7])
广播机制允许这样的二元运算能够在不同尺寸和形状的数组之间进行。 例如,我们可以用数组和一个标量相加(标量可以认为是一个零维数组):
a + 5
array([5, 6, 7])
我们可以认为上面的运算首先将标量扩展成了一个一维的数组[5, 5, 5]
,然后在和 a
进行了加法运算。
NumPy的广播方式并不是真的需要将元素复制然后扩展,但是这对于理解广播的运行方式很有帮助。
我们可以很简单的将上面的情形推广到更高纬度的数组上。下面我们使用广播将一个一维数组和一个二维数组进行加法运算:
M = np.ones((3, 3))
M
array([[1., 1., 1.], [1., 1., 1.], [1., 1., 1.]])
M + a
array([[1., 2., 3.], [1., 2., 3.], [1., 2., 3.]])
上例中一维数组a
在第二个维度上进行了扩展或者广播,这样才能符合M
的形状。
上面两个例子相对来说非常容易理解,但是当参与运算的两个数组都需要广播时,情况就相对复杂一些了。看下面的例子:
a = np.arange(3)
b = np.arange(3)[:, np.newaxis]
print(a)
print(b)
[0 1 2] [[0] [1] [2]]
a + b
array([[0, 1, 2], [1, 2, 3], [2, 3, 4]])
浅色格子代表的是广播后的值:再次说明,这些广播的值不会真正占用内存,只是为了辅助我们理解广播的机制。
M = np.ones((2, 3))
a = np.arange(3)
我们先看一下两个数组的形状:
M.shape = (2, 3)
a.shape = (3,)
依据规则1,数组a
的维度较少,因此首先对其进行维度扩增,我们在其最前面(最左边)增加一个维度,长度为1。此时两个数组的形状变为:
M.shape -> (2, 3)
a.shape -> (1, 3)
依据规则2,我们可以看到双方在第一维度上不相同,因此我们将第一维度具有长度1的a
的第一维度扩展为2。此时双方的形状变为:
M.shape -> (2, 3)
a.shape -> (2, 3)
经过变换之后,双方形状一致,可以进行加法运算了,我们可以预知最终结果的形状为(2, 3)
:
M + a
array([[1., 2., 3.], [1., 2., 3.]])
a = np.arange(3).reshape((3, 1))
b = np.arange(3)
开始时双方的形状为:
a.shape = (3, 1)
b.shape = (3,)
由规则1我们需要将数组b
扩增第一维度,长度为1:
a.shape -> (3, 1)
b.shape -> (1, 3)
由规则2我们需要将数组a
的第二维度扩展为3,还需要将数组b
的第一维度扩展为3,得到:
a.shape -> (3, 3)
b.shape -> (3, 3)
双方形状相同,可以进行运算:
a + b
array([[0, 1, 2], [1, 2, 3], [2, 3, 4]])
M = np.ones((3, 2))
a = np.arange(3)
这个例子和例子1有一点点区别,那就是本例中的M
是例子1中M
的转置矩阵。它们的形状是:
M.shape = (3, 2)
a.shape = (3,)
由规则1我们需要在数组a
上扩增第一维度,长度为1:
M.shape -> (3, 2)
a.shape -> (1, 3)
由规则2我们需要将数组a
的第一维度扩展为3才能与数组M
保持一致,除此之外双方都没有长度为1的维度了:
M.shape -> (3, 2)
a.shape -> (3, 3)
观察得到的形状,你可以发现这个结果满足规则3,双方的各维度长度不完全一致且不为1,因此无法完成广播,最终会产生错误:
M + a
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In [14], line 1 ----> 1 M + a ValueError: operands could not be broadcast together with shapes (3,2) (3,)
这里你可能会发现一个问题:如果广播的时候不一定按照最前面(最左边)维度的原则进行扩增维度的话,那不是很多的数组都可以进行广播计算吗?这样处理不是更灵活吗?例如上例中如果我们在数组a
的第二维度上扩增的话,那广播就能正确进行了。很可惜,广播并不会支持这种处理方式,虽然这种方法在某些情况下会更加灵活,但是在部分情况下会带来不确定性。如果你确实希望进行右维度扩增的话,你必须明确指定。利用我们在NumPy数组基础中介绍的np.newaxis
属性可以进行这个操作:
a[:, np.newaxis].shape
(3, 1)
M + a[:, np.newaxis]
array([[1., 1.], [2., 2.], [3., 3.]])
还要说明的是,上面的例子中我们都是使用加法进行说明,实际上广播可以应用到任何的二元ufunc上。例如下面我们采用logaddexp(a, b)
函数求值,这个函数计算的是$\log(e^a + e^b)$的值,使用这个函数能比采用原始的exp和log函数进行计算得到更高的精度:
np.logaddexp(M, a[:, np.newaxis])
array([[1.31326169, 1.31326169], [1.69314718, 1.69314718], [2.31326169, 2.31326169]])
更多关于通用函数的介绍,请复习使用Numpy计算:通用函数。
广播操作在本书后面很多例子中都会见到。因此这里我们看一些简单的例子,更好的说明它。
在前一节中,我们看到了ufuncs提供了我们可以避免使用Python循环的低效方式,而广播则大大扩展了这种能力。 一个常见的例子就是我们需要将数据集进行中心化。 例如我们我们进行了10次采样观测,每次都会得到3个数据值。 按照惯例(参见Scikit-Learn数据表现方式),我们可以将这些数据存成一个$10 \times 3$的数组:
X = np.random.random((10, 3))
我们使用mean
函数沿着第一维度求出每个特征的平均值:
Xmean = X.mean(0)
Xmean
array([0.63651576, 0.6188268 , 0.51705641])
下面我们就可以将数组X
减去它的各维度平均值就可以将其中心化(这里就是一个广播操作):
X_centered = X - Xmean
我们来检查一下结果的正确性,我们可以通过查看中心化后的数组在各特征上的平均值是够接近于0来进行判断:
X_centered.mean(0)
array([-1.44328993e-16, 0.00000000e+00, 8.88178420e-17])
考虑到机器精度情况,平均值已经等于0了。
广播还有一个很有用的场景,就是当你需要绘制一个二维函数的图像时。如果我们希望定义一个函数$z = f(x, y)$,广播可以被用来计算二维平面上每个网格的数值:
# x和y都是0~5范围平均分的50个点
x = np.linspace(0, 5, 50)
y = np.linspace(0, 5, 50)[:, np.newaxis]
z = np.sin(x) ** 10 + np.cos(10 + y * x) * np.cos(x)
算出z后,我们使用Matplotlib来画出这个二维数组(我们将在密度和轮廓图中详细介绍):
%matplotlib inline
import matplotlib.pyplot as plt
plt.imshow(z, origin='lower', extent=[0, 5, 0, 5],
cmap='viridis')
plt.colorbar();
上面的图形以一种极其吸引人的方式为我们展现了二维函数的分布情况。