# NumPy 的应用 - 2

# 数组的运算

使用 NumPy 最为方便的是当需要对数组元素进行运算时,不用编写循环代码遍历每个元素,所有的运算都会自动的矢量化(使用高效的、提前编译的底层代码来对数据序列进行数学操作)。简单的说就是,NumPy 中的数学运算和数学函数会自动作用于数组中的每个成员。

# 数组跟标量的运算

代码:

array35 = np.arange(1, 10)
print(array35 + 10)
print(array35 * 10)

输出:

[11 12 13 14 15 16 17 18 19]
[10 20 30 40 50 60 70 80 90]

# 数组跟数组的运算

代码:

array36 = np.array([1, 1, 1, 2, 2, 2, 3, 3, 3])
print(array35 + array36)
print(array35 * array36)
print(array35 ** array36)

输出:

[ 2  3  4  6  7  8 10 11 12]
[ 1  2  3  8 10 12 21 24 27]
[  1   2   3  16  25  36 343 512 729]

# 通用一元函数

通用函数是对 ndarray 中的数据执行元素级运算的函数。你可以将其看做普通函数(接收一个标量值作为参数,返回一个标量值)的矢量化包装器,如下所示。

代码:

print(np.sqrt(array35))
print(np.log2(array35))

输出:

[1.         1.41421356 1.73205081 2.         2.23606798 2.44948974
 2.64575131 2.82842712 3.        ]
[0.         1.         1.5849625  2.         2.32192809 2.5849625
 2.80735492 3.         3.169925  ]

表 1:通用一元函数

函数说明
abs / fabs求绝对值的函数
sqrt求平方根的函数,相当于 array ** 0.5
square求平方的函数,相当于 array ** 2
exp计算exe^x 的函数
log / log10 / log2对数函数( e 为底 / 10 为底 / 2 为底)
sign符号函数( 1 - 正数; 0 - 零; -1 - 负数)
ceil / floor上取整 / 下取整
isnan返回布尔数组,NaN 对应 True ,非 NaN 对应 False
isfinite / isinf判断数值是否为无穷大的函数
cos / cosh / sin三角函数
sinh / tan / tanh三角函数
arccos / arccosh / arcsin反三角函数
arcsinh / arctan / arctanh反三角函数
rint / round四舍五入函数

# 通用二元函数

代码:

array37 = np.array([[4, 5, 6], [7, 8, 9]])
array38 = np.array([[1, 2, 3], [3, 2, 1]])
print(array37 ** array38)
print(np.power(array37, array38))

输出:

[[  4  25 216]
 [343  64   9]]
[[  4  25 216]
 [343  64   9]]

表 2:通用二元函数

函数说明
add(x, y) / substract(x, y)加法函数 / 减法函数
multiply(x, y) / divide(x, y)乘法函数 / 除法函数
floor_divide(x, y) / mod(x, y)整除函数 / 求模函数
allclose(x, y)检查数组 xy 元素是否几乎相等
power(x, y)数组xx 的元素xix_i 和数组yy 的元素yiy_i,计算x_i^
maximum(x, y) / fmax(x, y)两两比较元素获取最大值 / 获取最大值(忽略 NaN)
minimum(x, y) / fmin(x, y)两两比较元素获取最小值 / 获取最小值(忽略 NaN)
dot(x, y)点积运算(数量积,通常记为\cdots,用于欧几里得空间(Euclidean space))
inner(x, y)内积运算(内积的含义要高于点积,点积相当于是内积在欧几里得空间 $$ 的特例,而内积可以推广到赋范向量空间,只要它满足平行四边形法则即可)
cross(x, y)叉积运算(向量积,通常记为×\times,运算结果是一个向量)
outer(x, y)外积运算(张量积,通常记为\bigotimes,运算结果通常是一个矩阵)
intersect1d(x, y)计算 xy 的交集,返回这些元素构成的有序数组
union1d(x, y)计算 xy 的并集,返回这些元素构成的有序数组
in1d(x, y)返回由判断 x 的元素是否在 y 中得到的布尔值构成的数组
setdiff1d(x, y)计算 xy 的差集,返回这些元素构成的数组
setxor1d(x, y)计算 xy 的对称差,返回这些元素构成的数组

补充说明:在二维空间内,两个向量A=[a1a2]\boldsymbol{A}=\begin{bmatrix} a_1 \\ a_2 \end{bmatrix}B=[b1b2]\boldsymbol{B}=\begin{bmatrix} b_1 \\ b_2 \end{bmatrix} 的叉积是这样定义的:A×B=a1a2b1b2=a1b2a2b1\boldsymbol{A}\times \boldsymbol{B}=\begin{vmatrix} a_1 \quad a_2 \\ b_1 \quad b_2 \end{vmatrix}=a_1b_2 - a_2b_1,其中a1a2b1b2\begin{vmatrix} a_1 \quad a_2 \\ b_1 \quad b_2 \end{vmatrix} 称为行列式。但是一定要注意,叉积并不等同于行列式,行列式的运算结果是一个标量,而叉积运算的结果是一个向量。如果不明白,我们可以看看三维空间两个向量,A=[a1a2a3]\boldsymbol{A}=\begin{bmatrix} a_1 \\ a_2 \\ a_3 \end{bmatrix}B=[b1b2b3]\boldsymbol{B}=\begin{bmatrix} b_1 \\ b_2 \\ b_3 \end{bmatrix} 的叉积是<i^a2a3b2b3,j^a1a3b1b3,k^a1a2b1b2>\left< \hat{i} \begin{vmatrix} a_2 \quad a_3 \\ b_2 \quad b_3 \end{vmatrix}, -\hat{j} \begin{vmatrix} a_1 \quad a_3 \\ b_1 \quad b_3 \end{vmatrix}, \hat{k} \begin{vmatrix} a_1 \quad a_2 \\ b_1 \quad b_2 \end{vmatrix} \right>,其中i^,j^,k^\hat{i}, \hat{j}, \hat{k} 代表每个维度的单位向量。

# 广播机制

上面的例子中,两个二元运算的数组形状是完全相同的,我们再来研究一下,两个形状不同的数组是否可以直接做二元运算或使用二元函数进行运算,请看下面的例子。

代码:

array39 = np.array([[0, 0, 0], [1, 1, 1], [2, 2, 2], [3, 3, 3]])
array40 = np.array([1, 2, 3])
array39 + array40

输出:

array([[1, 2, 3],
       [2, 3, 4],
       [3, 4, 5],
       [4, 5, 6]])

代码:

array41 = np.array([[1], [2], [3], [4]])
array39 + array41

输出:

array([[1, 1, 1],
       [3, 3, 3],
       [5, 5, 5],
       [7, 7, 7]])

通过上面的例子,我们发现形状不同的数组仍然有机会进行二元运算,但也绝对不是任意的数组都可以进行二元运算。简单的说,只有两个数组后缘维度相同或者其中一个数组后缘维度为 1 时,广播机制会被触发,而通过广播机制如果能够使两个数组的形状一致,才能进行二元运算。所谓后缘维度,指的是数组 shape 属性对应的元组中最后一个元素的值(从后往前数最后一个维度的值),例如,我们之前打开的图像对应的数组后缘维度为 3,3 行 4 列的二维数组后缘维度为 4,而有 5 个元素的一维数组后缘维度为 5。简单的说就是,后缘维度相同或者其中一个数组的后缘维度为 1,就可以应用广播机制;而广播机制如果能够使得数组的形状一致,就满足了两个数组对应元素做运算的需求,如下图所示。

# 其他常用函数

除了上面讲到的函数外,NumPy 中还提供了很多用于处理数组的函数, ndarray 对象的很多方法也可以通过直接调用函数来实现,下表给出了一些常用的函数。

表 3:NumPy 其他常用函数

函数说明
unique去除数组重复元素,返回唯一元素构成的有序数组
copy返回拷贝数组得到的数组
sort返回数组元素排序后的拷贝
split / hsplit / vsplit将数组拆成若干个子数组
stack / hstack / vstack将多个数组堆叠成新数组
concatenate沿着指定的轴连接多个数组构成新数组
append / insert向数组末尾追加元素 / 在数组指定位置插入元素
argwhere找出数组中非 0 元素的位置
extract / select / where按照指定的条件从数组中抽取或处理数组元素
flip沿指定的轴翻转数组中的元素
fromiter通过迭代器创建数组对象
fromregex通过读取文件和正则表达式解析获取数据创建数组对象
repeat / tile通过对元素的重复来创建新数组
roll沿指定轴对数组元素进行移位
resize重新调整数组的大小
place / put将数组中满足条件的元素 / 指定的元素替换为指定的值
partition用选定的元素对数组进行一次划分并返回划分后的数组

提示:上面的 resize 函数和 ndarray 对象的 resize 方法是有区别的, resize 函数在调整数组大小时会重复数组中的元素作为填补多出来的元素的值,而 ndarry 对象的 resize 方法是用 0 来填补多出来的元素。这些小细节不清楚暂时也不要紧,但是如果用到对应的功能了就要引起注意。

代码:

array42 = np.array([[1, 1, 1], [2, 2, 2], [3, 3, 3]])
array43 = np.array([[4, 4, 4], [5, 5, 5], [6, 6, 6]])
np.hstack((array42, array43))

输出:

array([[1, 1, 1, 4, 4, 4],
       [2, 2, 2, 5, 5, 5],
       [3, 3, 3, 6, 6, 6]])

代码:

np.vstack((array42, array43))

输出:

array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3],
       [4, 4, 4],
       [5, 5, 5],
       [6, 6, 6]])

代码:

np.concatenate((array42, array43))

输出:

array([[1, 1, 1],
       [2, 2, 2],
       [3, 3, 3],
       [4, 4, 4],
       [5, 5, 5],
       [6, 6, 6]])

代码:

np.concatenate((array42, array43), axis=1)

输出:

array([[1, 1, 1, 4, 4, 4],
       [2, 2, 2, 5, 5, 5],
       [3, 3, 3, 6, 6, 6]])

# 矩阵运算

NumPy 中提供了专门用于线性代数(linear algebra)的模块和表示矩阵的类型 matrix ,当然我们通过二维数组也可以表示一个矩阵,官方并不推荐使用 matrix 类而是建议使用二维数组,而且有可能在将来的版本中会移除 matrix 类。无论如何,利用这些已经封装好的类和函数,我们可以轻松愉快的实现线性代数中很多的操作。

# 线性代数快速回顾

  1. 向量也叫矢量,是一个同时具有大小和方向,且满足平行四边形法则的几何对象。与向量相对的概念叫标量数量,标量只有大小、绝大多数情况下没有方向。
  2. 向量可以进行数乘点积叉积等运算。
  3. 行列式由向量组成,它的性质可以由向量解释。
  4. 行列式可以使用行列式公式计算:det(A)=n!±a1αa2βanωdet(\boldsymbol{A})=\sum_{n!} \pm {a_{1\alpha}a_{2\beta} \cdots a_{n\omega}}
  5. 高阶行列式可以用代数余子式展开成多个低阶行列式,如:det(A)=a11C11+a12C12++a1nC1ndet(\boldsymbol{A})=a_{11}C_{11}+a_{12}C_{12}+ \cdots +a_{1n}C_{1n}
  6. 矩阵是由一系列元素排成的矩形阵列,矩阵里的元素可以是数字、符号或数学公式。
  7. 矩阵可以进行加法减法数乘乘法转置等运算。
  8. 逆矩阵A1\boldsymbol{A^{-1}} 表示,AA1=A1A=I\boldsymbol{A}\boldsymbol{A^{-1}}=\boldsymbol{A^{-1}}\boldsymbol{A}=\boldsymbol{I};没有逆矩阵的方阵是奇异矩阵
  9. 如果一个方阵是满秩矩阵 (矩阵的秩等于矩阵的阶数),该方阵对应的线性方程有唯一解。

说明矩阵的秩是指矩阵中线性无关的行 / 列向量的最大个数,同时也是矩阵对应的线性变换的像空间的维度。

# NumPy 中矩阵相关函数

  1. 创建矩阵对象。

    代码:

    # matrix构造函数可以传入类数组对象也可以传入字符串
    m1 = np.matrix('1 2 3; 4 5 6')
    m1
    

    输出:

    matrix([[1, 2, 3],
            [4, 5, 6]])
    

    代码:

    # asmatrix函数也可以写成mat函数,它们其实是同一个函数
    m2 = np.asmatrix(np.array([[1, 1], [2, 2], [3, 3]]))
    m2
    

    输出:

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

    代码:

    m1 * m2
    

    输出:

    matrix([[14, 14],
            [32, 32]])
    

    说明:注意 matrix 对象和 ndarray 对象乘法运算的差别,如果两个二维数组要做矩阵乘法运算,应该使用 @ 运算符或 matmul 函数,而不是 * 运算符。

  2. 矩阵对象的属性。

    | 属性 | 说明 |
    | ------- | ----------------------------------------- |
    | A | 获取矩阵对象对应的 ndarray 对象 |
    | A1 | 获取矩阵对象对应的扁平化后的 ndarray 对象 |
    | I | 可逆矩阵的逆矩阵 |
    | T | 矩阵的转置 |
    | H | 矩阵的共轭转置 |
    | shape | 矩阵的形状 |
    | size | 矩阵元素的个数 |

  3. 矩阵对象的方法。

矩阵对象的方法跟之前讲过的 ndarray 数组对象的方法基本差不多,此处不再进行赘述。

# NumPy 的线性代数模块

NumPy 的 linalg 模块中有一组标准的矩阵分解运算以及诸如求逆和行列式之类的函数,它们跟 MATLAB 和 R 等语言所使用的是相同的行业标准线性代数库,下面的表格列出了 numpy 以及 linalg 模块中常用的跟线性代数相关的函数。

函数说明
diag以一维数组的形式返回方阵的对角线元素或将一维数组转换为方阵(非对角元素元素为 0)
vdot向量的点积
dot数组的点积
inner数组的内积
outer数组的叉积
trace计算对角线元素的和
norm求模(范数)运算
det计算行列式的值(在方阵上计算会得到一个标量)
matrix_rank计算矩阵的秩
eig计算矩阵的特征值(eigenvalue)和特征向量(eigenvector)
inv计算非奇异矩阵(nn 阶方阵)的逆矩阵
pinv计算矩阵的摩尔 - 彭若斯(Moore-Penrose)广义逆
qrQR 分解(把矩阵分解成一个正交矩阵与一个上三角矩阵的积)
svd计算奇异值分解(singular value decomposition)
solve解线性方程组Ax=b\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b},其中A\boldsymbol{A} 是一个方阵
lstsq计算Ax=b\boldsymbol{A}\boldsymbol{x}=\boldsymbol{b} 的最小二乘解

大家如果有兴趣可以用下面的代码验证上面的函数。

代码:

m3 = np.array([[1., 2.], [3., 4.]])
np.linalg.inv(m3)

输出:

array([[-2. ,  1. ],
       [ 1.5, -0.5]])

代码:

m4 = np.array([[1, 3, 5], [2, 4, 6], [4, 7, 9]])
np.linalg.det(m4)

输出:

2

代码:

# 解线性方程组ax=b
# 3*x1 + x2= 9,x1 + 2*x2 = 8
a = np.array([[3,1], [1,2]])
b = np.array([9, 8])
np.linalg.solve(a, b)

输出:

array([2., 3.])