NumPy 是一个为 Python 提供高性能向量、矩阵和高维数据结构的科学计算库。它通过 C 和 Fortran 实现,因此用向量和矩阵建立方程并实现数值计算有非常好的性能。NumPy 基本上是所有使用 Python 进行数值计算的框架和包的基础,例如 TensorFlow 和 PyTorch,构建机器学习模型最基础的内容就是学会使用 NumPy 搭建计算过程。

深入理解 NumPy

广播机制

广播操作是 NumPy 非常重要的一个特点,它允许 NumPy 扩展矩阵间的运算。例如它会隐式地把一个数组的异常维度调整到与另一个算子相匹配的维度以实现维度兼容。例如将一个维度为 [3,2] 的矩阵与另一个维度为 [3,1] 的矩阵相加是合法的,NumPy 会自动将第二个矩阵扩展到等同的维度。

为了定义两个形状是否是可兼容的,NumPy 从最后开始往前逐个比较它们的维度大小。在这个过程中,如果两者的对应维度相同,或者其一(或者全是)等于 1,则继续进行比较,直到最前面的维度。若不满足这两个条件,程序就会报错。

如下展示了一个广播操作:

>>>a = np.array([1.0,2.0,3.0,4.0, 5.0, 6.0]).reshape(3,2)
>>>b = np.array([3.0])
>>>a * b

array([[  3.,   6.],
       [  9.,  12.],
       [ 15.,  18.]])

高级索引

NumPy 比一般的 Python 序列提供更多的索引方式。除了之前看到的用整数和截取的索引,数组可以由整数数组和布尔数组 indexed。

通过数组索引

如下我们可以根据数组 i 和 j 索引数组 a 中间的元素,其中输出数组保持索引的 shape。

>>> a = np.arange(12)**2                       # the first 12 square numbers
>>> i = np.array( [ 1,1,3,8,5 ] )              # an array of indices
>>> a[i]                                       # the elements of a at the positions i
array([ 1,  1,  9, 64, 25])

>>> j = np.array( [ [ 3, 4], [ 9, 7 ] ] )      # a bidimensional array of indices
>>> a[j]                                       # the same shape as j
array([[ 9, 16],
       [81, 49]])

当使用多维数组作为索引时,每一个维度就会索引一次原数组,并按索引的 shape 排列。下面的代码展示了这种索引方式,palette 可以视为简单的调色板,而数组 image 中的元素则表示索引对应颜色的像素点。

>>> palette = np.array( [ [0,0,0],                # black
...                       [255,0,0],              # red
...                       [0,255,0],              # green
...                       [0,0,255],              # blue
...                       [255,255,255] ] )       # white
>>> image = np.array( [ [ 0, 1, 2, 0 ],           # each value corresponds to a color in the palette
...                     [ 0, 3, 4, 0 ]  ] )
>>> palette[image]                            # the (2,4,3) color image
array([[[  0,   0,   0],
        [255,   0,   0],
        [  0, 255,   0],
        [  0,   0,   0]],
       [[  0,   0,   0],
        [  0,   0, 255],
        [255, 255, 255],
        [  0,   0,   0]]])
       [81, 49]])

我们也可以使用多维索引获取数组中的元素,多维索引的每个维度都必须有相同的形状。如下多维数组 i 和 j 可以分别作为索引 a 中第一个维度和第二个维度的参数,例如 a[i, j] 分别从 i 和 j 中抽取一个元素作为索引 a 中元素的参数。

>>> a = np.arange(12).reshape(3,4)
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>> i = np.array( [ [0,1],                        # indices for the first dim of a
...                 [1,2] ] )
>>> j = np.array( [ [2,1],                        # indices for the second dim
...                 [3,3] ] )
>>>
>>> a[i,j]                                     # i and j must have equal shape
array([[ 2,  5],
       [ 7, 11]])
>>>
>>> a[i,2]
array([[ 2,  6],
       [ 6, 10]])
>>>
>>> a[:,j]                                     # i.e., a[ : , j]
array([[[ 2,  1],
        [ 3,  3]],
       [[ 6,  5],
        [ 7,  7]],
       [[10,  9],
        [11, 11]]])

同样,我们把 i 和 j 放在一个序列中,然后用它作为索引:

>>> l = [i,j]
>>> a[l]                                       # equivalent to a[i,j]
array([[ 2,  5],
       [ 7, 11]])

然而,我们不能如上把 i 和 j 放在一个数组中作为索引,因为数组会被理解为索引 a 的第一维度。

>>> s = np.array( [i,j] )
>>> a[s]                                       # not what we want
Traceback (most recent call last):
  File "<stdin>", line 1, in ?
IndexError: index (3) out of range (0<=index<=2) in dimension 0
>>>
>>> a[tuple(s)]                                # same as a[i,j]
array([[ 2,  5],
       [ 7, 11]])

另一个将数组作为索引的常用方法是搜索时间序列的最大值:

>>> time = np.linspace(20, 145, 5)                 # time scale
>>> data = np.sin(np.arange(20)).reshape(5,4)      # 4 time-dependent series
>>> time
array([  20.  ,   51.25,   82.5 ,  113.75,  145.  ])
>>> data
array([[ 0.        ,  0.84147098,  0.90929743,  0.14112001],
       [-0.7568025 , -0.95892427, -0.2794155 ,  0.6569866 ],
       [ 0.98935825,  0.41211849, -0.54402111, -0.99999021],
       [-0.53657292,  0.42016704,  0.99060736,  0.65028784],
       [-0.28790332, -0.96139749, -0.75098725,  0.14987721]])
>>>
>>> ind = data.argmax(axis=0)                  # index of the maxima for each series
>>> ind
array([2, 0, 3, 1])
>>>
>>> time_max = time[ind]                       # times corresponding to the maxima
>>>
>>> data_max = data[ind, range(data.shape[1])] # => data[ind[0],0], data[ind[1],1]...
>>>
>>> time_max
array([  82.5 ,   20.  ,  113.75,   51.25])
>>> data_max
array([ 0.98935825,  0.84147098,  0.99060736,  0.6569866 ])
>>>
>>> np.all(data_max == data.max(axis=0))
True

你也可以用数组索引作为一个分配目标:

>>> a = np.arange(5)
>>> a
array([0, 1, 2, 3, 4])
>>> a[[1,3,4]] = 0
>>> a
array([0, 0, 2, 0, 0])

然而,当索引列表中有重复时,赋值任务会执行多次,并保留最后一次结果。

>>> a = np.arange(5)
>>> a[[0,0,2]]=[1,2,3]
>>> a
array([2, 1, 3, 3, 4])

这是合理的,但注意如果你使用 Python 的 +=创建,可能不会得出预期的结果:

>>> a = np.arange(5)
>>> a[[0,0,2]]+=1
>>> a
array([1, 1, 3, 3, 4])

虽然 0 在索引列表中出现两次,第 0 个元素只会增加一次。这是因为 Python 中「a+=1」等于「a = a + 1」.

用布尔数组做索引

当我们索引数组元素时,我们在提供索引列表。但布尔值索引是不同的,我们需要清楚地选择被索引数组中哪个元素是我们想要的哪个是不想要的。

布尔索引需要用和原数组相同 shape 的布尔值数组,如下只有在大于 4 的情况下才输出 True,而得出来的布尔值数组可作为索引。

>>> a = np.arange(12).reshape(3,4)
>>> b = a > 4
>>> b                                          # b is a boolean with a s shape
array([[False, False, False, False],
       [False,  True,  True,  True],
       [ True,  True,  True,  True]])
>>> a[b]                                       # 1d array with the selected elements
array([ 5,  6,  7,  8,  9, 10, 11])

这个性质在任务中非常有用,例如在 ReLu 激活函数中,只有大于 0 才输出激活值,因此我们就能使用这种方式实现 ReLU 激活函数。

>>> a[b] = 0                                   # All elements of  a  higher than 4 become 0
>>> a
array([[0, 1, 2, 3],
       [4, 0, 0, 0],
       [0, 0, 0, 0]])

第二种使用布尔索引的方法与整数索引更加相似的;在数组的每个维度中,我们使用一维布尔数组选择我们想要的截取部分:

>>> a = np.arange(12).reshape(3,4)
>>> b1 = np.array([False,True,True])             # first dim selection
>>> b2 = np.array([True,False,True,False])       # second dim selection
>>>
>>> a[b1,:]                                   # selecting rows
array([[ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>>
>>> a[b1]                                     # same thing
array([[ 4,  5,  6,  7],
       [ 8,  9, 10, 11]])
>>>
>>> a[:,b2]                                   # selecting columns
array([[ 0,  2],
       [ 4,  6],
       [ 8, 10]])
>>>
>>> a[b1,b2]                                  # a weird thing to do
array([ 4, 10])

注意一维布尔数组的长度必须和想截取轴的长度相同。在上面的例子中,b1 的长度 3、b2 的长度为 4,它们分别对应于 a 的第一个维度与第二个维度。

线性代数

简单的数组运算

如下仅展示了简单的矩阵运算更多详细的方法可在实践中遇到在查找 API。如下展示了矩阵的转置、求逆、单位矩阵、矩阵乘法、矩阵的迹、解线性方程和求特征向量等基本运算:

 

>>> import numpy as np
>>> a = np.array([[1.0, 2.0], [3.0, 4.0]])
>>> print(a)
[[ 1.  2.]
 [ 3.  4.]]

>>> a.transpose()
array([[ 1.,  3.],
       [ 2.,  4.]])

>>> np.linalg.inv(a)
array([[-2. ,  1. ],
       [ 1.5, -0.5]])

>>> u = np.eye(2) # unit 2x2 matrix; "eye" represents "I"
>>> u
array([[ 1.,  0.],
       [ 0.,  1.]])
>>> j = np.array([[0.0, -1.0], [1.0, 0.0]])

>>> np.dot (j, j) # matrix product
array([[-1.,  0.],
       [ 0., -1.]])

>>> np.trace(u)  # trace
2.0

>>> y = np.array([[5.], [7.]])
>>> np.linalg.solve(a, y)
array([[-3.],
       [ 4.]])

>>> np.linalg.eig(j)
(array([ 0.+1.j,  0.-1.j]), array([[ 0.70710678+0.j        ,  0.70710678-0.j        ],
       [ 0.00000000-0.70710678j,  0.00000000+0.70710678j]]))

Parameters:
    square matrix
Returns
    The eigenvalues, each repeated according to its multiplicity.
    The normalized (unit "length") eigenvectors, such that the
    column ``v[:,i]`` is the eigenvector corresponding to the
    eigenvalue ``w[i]`` .