• 本文是NumPy快速入门系列的最后一篇。
• 本文提到的矩阵均指 np.array 对象，而非 np.matrix 对象，后者已经不再推荐使用。
• NumPy中，与线性代数有关的模块是 np.linalg （linear algebra），我们通常会简写为 LA，即

from numpy import linalg as LA

一、@ 算子

Numpy中，@ 算子可用来计算两个矩阵的乘积：

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A @ B)
# [[19 22]
#  [43 50]]

np.matmul 也可用来计算两个矩阵的乘积，但这里更推荐使用 @。

二、向量乘积与矩阵乘积

函数	作用
np.inner(a, b)	a和b均是一维数组（向量）；计算两个向量的内积
np.outer(a, b)	计算两个向量的外积
np.kron(a, b)	计算两个数组的克罗内克积
LA.matrix_power(A, n)	$A$ 是方阵；计算 $A^n$

2.1 np.inner()

设 $a$ 和 $b$ 是两个（列）向量，则 np.inner(a, b) 相当于 $a^{\mathrm{T}}b$。

a = np.array([1, 2, 3])
b = np.array([3, 2, 1])
print(np.inner(a, b))
# 10

2.2 np.outer()

设 $a$ 和 $b$ 是两个（列）向量，则 np.outer(a, b) 相当于 $a{b}^{\mathrm{T}}$。

a = np.array([1, 2, 3])
b = np.array([3, 2, 1])
print(np.outer(a, b))
# [[3 2 1]
#  [6 4 2]
#  [9 6 3]]

2.3 np.kron()

有关克罗内克积的定义可参考Wikipedia. Kronecker product。

大多数情况下我们会计算两个矩阵的克罗内克积，设 $A$ 和 $B$ 是两个矩阵，则 np.kron(a, b) 相当于 $A\otimes B$。

A = np.array([[1, 2], [3, 4]])
B = np.array([[0, 5], [6, 7]])
print(np.kron(A, B))
# [[ 0  5  0 10]
#  [ 6  7 12 14]
#  [ 0 15  0 20]
#  [18 21 24 28]]

2.4 LA.matrix_power()

我们使用幂零矩阵来观察效果：

A = np.eye(4, k=1)
print(A)
# [[0. 1. 0. 0.]
#  [0. 0. 1. 0.]
#  [0. 0. 0. 1.]
#  [0. 0. 0. 0.]]
print(LA.matrix_power(A, 2))
# [[0. 0. 1. 0.]
#  [0. 0. 0. 1.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]]
print(LA.matrix_power(A, 3))
# [[0. 0. 0. 1.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]]
print(LA.matrix_power(A, 4))
# [[0. 0. 0. 0.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]
#  [0. 0. 0. 0.]]

三、行列式、秩、迹、范数

函数	作用
LA.det(A)	计算矩阵 A 的行列式
LA.matrix_rank(A)	计算矩阵 A 的秩
np.trace(A)	计算矩阵 A 的迹
LA.norm(a, ord=None)	ord控制范数；计算向量或矩阵的范数

3.1 LA.det()

A = np.array([[1, 2], [3, 4]])
print(LA.det(A))
# -2.0000000000000004

3.2 LA.matrix_rank()

A = np.array([[1, 2], [3, 4]])
B = np.array([[1, 1], [1, 1]])
C = np.array([[0, 0], [0, 0]])
print(LA.matrix_rank(A))
print(LA.matrix_rank(B))
print(LA.matrix_rank(C))
# 2
# 1
# 0

3.3 np.trace()

A = np.identity(4)
print(np.trace(A))
# 4.0

3.4 LA.norm()

在进一步讲解用法之前，我们先来回顾一下矩阵和向量的常见范数：

矩阵范数	表达式
Frobenius 范数	$\Vert A{\Vert }_F=(\sum _{i,j}\mid a_{ij}{\mid }^2{)}^{1/2}$
$1$ - 范数	$\Vert A{\Vert }_1=\underset{j}{\max }\sum _i\mid a_{ij}\mid$
$2$ - 范数	$\Vert A{\Vert }_2=\sqrt{{\lambda }_{max}(A^{H}A)}$
$\infty$ - 范数	$\Vert A{\Vert }_{\infty }=\underset{i}{\max }\sum _j\mid a_{ij}\mid$
核范数	$\Vert A{\Vert }_{\ast }=\mathrm{t}\mathrm{r}(\sqrt{A^{H}A})$

向量范数	表达式
$p$ 范数（$1\le p<+\infty$）	$\Vert \mathbf{x}{\Vert }_p=(\sum _i\mid x_i{\mid }^p{)}^{1/p}$
$\infty$ - 范数	$\Vert \mathbf{x}{\Vert }_{\infty }=\underset{i}{\max }\mid x_i\mid$

再来看一下 ord 的常用取值：

ord	矩阵范数	向量范数
None	Frobenius范数	$2$ - 范数
‘nuc’	核范数	——
inf	max(sum(abs(a), axis=1))	max(abs(a))
1	max(sum(abs(a), axis=0))	sum(abs(a))
2	$2$ - 范数	sum(abs(a) $\ast \ast$ 2) $\ast \ast$ (1./2)

具体代码示例不再给出，读者若有需要可自行练习。

四、特征值与特征向量

函数	作用
LA.eig(A)	计算方阵A的特征值和特征向量
LA.eigh(A)	计算矩阵A的特征值和特征向量（A是实对称矩阵或厄米特矩阵）
LA.eigvals(A)	计算方阵A的特征值
LA.eigvalsh(A)	计算方阵A的特征值（A是实对称矩阵或厄米特矩阵）

4.1 LA.eig()

我们先来回顾一下相关的定义。

这里只考虑实数域。设 $A$ 是 $n$ 阶方阵，如果存在常数 $\lambda$ 和非零向量 $x$ 使得

$$Ax=\lambda x$$

则称 $\lambda$ 是 $A$ 的特征值，$x$ 是 $A$ 属于特征值 $\lambda$ 的特征向量。

不考虑重复的情形，$A$ 一定有 $n$ 个特征值，记为 ${\lambda }_1,\cdots ,{\lambda }_n$，相应的特征向量记为 ${\mathbf{p}}_1,\cdots ,{\mathbf{p}}_n$。

令 $\mathbf{\lambda }=[{\lambda }_1,\cdots ,{\lambda }_n]$，$\mathbf{P}=[{\mathbf{p}}_1,\cdots ,{\mathbf{p}}_n]$，可知 $\mathbf{\lambda }$ 是一个一维数组，$\mathbf{P}$ 是一个二维数组。对 $A$ 使用 LA.eig(A) 会返回一个元组 $(\mathbf{\lambda },\mathbf{P})$，我们通常会使用两个参数 w, v 进行接收：

A = np.diag((3, 2, 1))
w, v = LA.eig(A)
print(w)
# [3. 2. 1.]
print(v)
# [[1. 0. 0.]
#  [0. 1. 0.]
#  [0. 0. 1.]]

结合上述定义，我们可以得知：

v[:, i] 是 $A$ 属于特征值 w[i] 的特征向量。

我们再来看一个例子

A = np.array([[1, -1, 3],[0, 1, 2],[0, 0, 2]
])
w, v = LA.eig(A)
print(w)
# [1. 1. 2.]
print(v)
# [[1.00000000e+00 1.00000000e+00 4.08248290e-01]
#  [0.00000000e+00 2.22044605e-16 8.16496581e-01]
#  [0.00000000e+00 0.00000000e+00 4.08248290e-01]]

从结果中不难看出有 $\Vert {\mathbf{p}}_i\Vert =1$，即返回后的每个特征向量 v[:, i] 都是标准化了的。此外，从第一个例子也能看出，我们的特征值列表 w 并不是从小到大进行排列的。

4.2 LA.eigh()

该函数返回的特征值列表 w 是从小到大进行排列的（实对称/厄米特矩阵的特征值总是实数）。

对于实对称矩阵和厄米特矩阵，我们当然也可以采用 LA.eig 进行计算，但毫无疑问使用 LA.eigh 是更好的，因为它利用了实对称/厄米特矩阵的特性从而加快了计算速度。

具体代码示例不再给出，读者若有需要可自行练习。

4.3 LA.eigvals()

该函数只返回特征值列表 w，且 w 不一定是从小到大进行排列的。

A = np.diag((3, 2, 1))
print(LA.eigvals(A))
# [3. 2. 1.]

4.4 LA.eigvalsh()

该函数返回的特征值列表 w 是从小到大进行排列的。

具体代码示例不再给出，读者若有需要可自行练习。