Numpy 学习笔记——3. 数组的线性代数运算

NumPy 提供了比较丰富的线性代数运算功能，主要集中在 numpy.linalg 模块中，下面罗列一些基础的操作，主要针对一维和二维数组，对于其它高维数组的语义可能不太符合直觉。

重要函数

点积 np.dot

np.dot() 函数可以计算两个向量的点积，最基本的用法例如

a = np.array([1, 2, 3], dtype=np.float64)
b = np.array([4, 5, 6], dtype=np.float64)
np.dot(a, b)  # np.float64(32.0)
# 1*4 + 2*5 + 3*6 = 32

np.dot(a, b)的完整语义需要根据两个数组的维度来区分：

第一类情况：
- 如果其中一个是标量，那么就是进行（逐元素的）乘法，等效于 np.multiply(a, b) 或 a * b。
第二类情况：
- 如果 a 是 $N$ 维数组，b 是一维数组，那么会沿着两者的最后一个轴求和，得到 $N-1$ 维数组。
- （特例，标准情况）如果两者都是尺寸为 (m,) 的一维数组，那么就是标准的向量内积，结果是一个标量。
- （特例）如果 a 是尺寸为 (m, n) 的二维数组，b 是尺寸为 (n,) 的一维数组，那么就是标准的矩阵乘以向量，得到尺寸为 (m,) 的一维数组。
第三类情况：
- 如果 a 是 $N$ 维数组，b 是 $M$ 维数组（$M \ge 2$），那么会沿着 a 的最后一个轴和 b 的倒数第二个轴求和，得到 $N+M-2$ 维数组。
- （特例）如果 a 是尺寸为 (m, n) 的二维数组，b 是尺寸为 (n, k) 的二维数组，那么就是在进行标准的矩阵乘法，得到尺寸为 (m, k) 的二维数组，等效于 np.matmul(a, b) 或 a @ b；

注意：

对于不含标量的情形，要求两个数组在求和运算的轴有相同的维数，否则运算会抛出异常。
如果 a 和 b 都是标量或者都是一维数组，那么返回值是标量，否则返回 ndarray 数组。

涉及标量的情况例如

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

np.dot(2,a) # array([2., 4., 6.])
np.dot(a,2) # array([2., 4., 6.])

np.dot() 的处理实际上是希望用户忽略行向量（无论形状是(n,)还是(1,n)）和列向量（形状为(n,1)）的差异，
下面两种情况都是在进行矩阵乘以向量的操作：

形状为 (m,n) 的二维数组（矩阵）和形状为(n,)的一维数组（行向量）进行 dot 运算，得到的结果是一个形状为 (m,) 的一维数组（行向量）。

1
2
3

a = np.array([[1, 2], [3, 4]], dtype=np.float64)
b = np.array([10, 20], dtype=np.float64)
np.dot(a, b)  # array([ 50., 110.])

形状为 (m,n) 的二维数组（矩阵）和形状为(n,1)的二维数组（列向量）进行 dot 运算，得到的结果是一个形状 (m, 1) 的二维数组（列向量）。

a = np.array([[1, 2], [3, 4]], dtype=np.float64)
b = np.array([[10], [20]], dtype=np.float64)
np.dot(a, b)
# array([[ 50.],
#        [110.]])

除了 numpy.dot() 函数，还有几个语义非常类似的函数：

numpy.vdot()：如果第一个数组是复值，那么会对其进行共轭转置；
numpy.tensordot()：可以在任意指定的轴进行求和；
numpy.inner()：总是对两个数组的最后一个轴求和。

矩阵乘法 np.matmul

np.matmul() 函数可以计算两个矩阵的乘积，最基本的用法例如

A = np.array([[1, 2], [3, 4]], dtype=np.float64)
B = np.array([[5, 6], [7, 8]], dtype=np.float64)

C = np.matmul(A, B) # shape: (2, 2)
# array([[19., 22.],
#        [43., 50.]])

np.matmul(a, b) 可以使用符号 a @ b 代替，注意 @ 和 * 处于同一优先级。

np.matmul(a, b) 的完整语义需要根据两个数组的维度来区分：

（标准情况）如果 a 是尺寸为 (m, n) 的二维数组，b 是尺寸为 (n, k) 的二维数组，那么就是在进行标准的矩阵乘法，得到尺寸为 (m, k) 的二维数组；
（自动加轴）
- 如果 a 是尺寸为 (n,) 的一维数组，会将其形状变成 (1, n) 的二维数组（行向量）参与运算，在最终中会将添加的轴移除；
- 如果 b 是尺寸为 (n,) 的一维数组，会将其形状变成 (n, 1) 的二维数组（列向量）参与运算，在最终中会将添加的轴移除；
对于涉及更高维度的情况，np.matmul() 会将其视作矩阵的堆叠，见下文的补充说明。

注意：

要求第一个数组的最后一个轴和第二个数组的倒数第二个轴有相同的维数，否则运算会抛出异常。
如果 a 和 b 都是一维数组，那么返回值是标量，否则返回 ndarray 数组。

自动加轴的情况例如

a = np.array([1, 2], dtype=np.float64)
b = np.array([[3, 4], [5, 6]], dtype=np.float64)
np.matmul(a, b)  # shape: (2,)
# array([13., 16.])

c = np.array([[3, 4], [5, 6]], dtype=np.float64)
d = np.array([1, 2], dtype=np.float64)
np.matmul(c, d)  # shape: (2,)
# array([11., 17.])

np.matmul() 与 np.dot() 的语义大体相同，但是有几个不同点：

np.matmul() 不支持传入标量，会直接报错；
在处理高维数组时，np.dot() 只是沿着固定的轴求和，而 np.matmul() 会将其视作矩阵的堆叠，即一组矩阵与一组矩阵相乘。

下面通过几个例子可以说明两者对于高维情况的处理差异：np.dots() 始终只是沿着第一个数组的倒数第一个轴和第二个数组的倒数第二个轴求和

a = np.ones([3, 2, 7, 4])
c = np.ones([1, 5, 4, 3])
np.dot(a, c).shape # shape: (3, 2, 7, 1, 5, 3)

a = np.ones([3, 2, 7, 4])
c = np.ones([3, 2, 4, 3])
np.dot(a, c).shape # shape: (3, 2, 7, 3, 2, 3)

但是 np.matmul() 会将其视作两组矩阵的批量相乘

1
2
3

a = np.ones([3, 2, 7, 4])
c = np.ones([3, 2, 4, 3])
np.matmul(a, c).shape # shape: (3, 2, 7, 3)

这可以视作以矩阵为元素的形状为(3,2)的两个“矩阵数组”的对应项进行矩阵乘法。
下面这种情况则会报错，因为两个“矩阵数组”的形状不相等

1
2
3

a = np.ones([3, 2, 7, 4])
c = np.ones([1, 5, 4, 3])
np.matmul(a, c).shape # error

范数 np.linalg.norm

np.linalg.norm() 是用于计算向量或矩阵的范数的通用函数，分成向量和矩阵两种情况进行讨论，常见的情况包括：

如果输入一维数组：
- ord=1 向量 $L^1$ 范数（绝对值求和）
- ord=2 向量 $L^2$ 范数（平方和开根号）
- ord=p 向量 $L^p$ 范数
- ord=np.inf 向量 $L^\infty$ 范数（绝对值最大值）
如果输入二维数组：
- ord='fro' 矩阵的 Frobenius 范数（所有元素平方和开根号）
- ord=1 矩阵的 $L^1$ 范数（绝对值列和最大值）
- ord=2 矩阵的 $L^2$ 范数（谱范数，最大奇异值）
- ord=np.inf 矩阵的 $L^\infty$ 范数（绝对值行和最大值）

向量范数例如

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

np.linalg.norm(a, ord=1)  # np.float64(6.0)
np.linalg.norm(a, ord=2)  # np.float64(3.7416573867739413)
np.linalg.norm(a, ord=np.inf)  # np.float64(3.0)

矩阵范数例如

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

np.linalg.norm(a, ord="fro")  # np.float64(5.477225575051661)
np.linalg.norm(a, ord=1)  # np.float64(6.0)
np.linalg.norm(a, ord=2)  # np.float64(5.464985704219043)
np.linalg.norm(a, ord=np.inf)  # np.float64(7.0)

注意：

在没有指定ord时，会将数组展平为一维数组，然后计算向量的 $L^2$ 范数，此时对于矩阵也相当于 Frobenius 范数。
大部分情况只支持一维或二维数组，对于高维数组，要么使用ord=None，要么使用参数axis指定具体的轴进行运算。

其它函数

除了最常用的几个函数，其它线性代数运算的语义并没有那么复杂，因此只提供基本的例子即可，需要时再详细查阅文档即可。

线性方程组求解

求解线性方程组例如

A = np.array([[3, 1], [1, 2]], dtype=np.float64)
b = np.array([9, 8], dtype=np.float64)

x = np.linalg.solve(A, b)  # array([2., 3.])

这要求矩阵必须是方阵且满秩。

矩阵求逆

求逆矩阵/伪逆矩阵

1
2
3

A = np.array([[1, 2], [3, 4]], dtype=np.float64)
A_inv = np.linalg.inv(A)
A_pinv = np.linalg.pinv(A) # Moore-Penrose

显然前者要求矩阵可逆，如果矩阵可逆，那么伪逆矩阵也等于逆矩阵。

矩阵的秩、行列式、迹

直接给例子即可

A = np.array([[3, 1], [1, 2]], dtype=np.float64)

np.linalg.matrix_rank(A)  # np.int64(2)
np.linalg.det(A) # np.float64(5.000000000000001)
np.trace(A) # np.float64(5.0)

注意秩的返回结果是整数。

矩阵特征值分解

获取矩阵特征值与特征向量例如

A = np.array([[4, -2], [1, 1]], dtype=np.float64)

eigvals, eigvecs = np.linalg.eig(A)
# eigvals: array([3., 2.])
# eigvecs:
#   array([[0.89442719, 0.70710678],
#          [0.4472136 , 0.70710678]])

# A @ v = λ * v
for i in range(A.shape[0]):
    lmd = eigvals[i]
    v = eigvecs[:, i]
    print(f"{lmd=}, {v=}", np.allclose(A @ v, lmd * v))  # True

其中 eigvals 是特征值（可能为复数）组成的一维数组，eigvecs 是特征向量（按列拼接）组成的二维数组。

注：如果输入的是实对称矩阵（或复 Hermitian 矩阵），使用 np.linalg.eigh() 函数的效果更好且效率更高。

矩阵奇异值分解

获取矩阵的奇异值分解例如

A = np.random.rand(4, 3)

U, S, Vh = np.linalg.svd(A)
# U shape: (4, 4)
# S shape: (3,)
# Vh shape: (3, 3)

Sigma = np.zeros(A.shape)
Sigma[: len(S), : len(S)] = np.diag(S)
# Sigma shape: (4, 3)

np.allclose(A, U @ Sigma @ Vh)  # True

注意：

返回的 S 是奇异值组成的一维数组而不是对角阵。
支持full_matrices=False选项用于进行精简的SVD分解。

MATLAB vs Numpy

由于同时使用MATLAB和Numpy，它们之间有很多区别，有必要梳理一下。

一些关键的区别包括：

MATLAB 使用的索引从 1 开始；Numpy 使用的索引从 0 开始。
MATLAB 将绝大多数的数据都视作数组，标量都是(1,1)数组，行向量和列向量有明显区别；Numpy 的一维数组尺寸是 (n,)。
MATLAB 使用懒拷贝提高运算效率；Numpy 则使用视图提高运算效率，需要拷贝时最好使用 ndarray.copy() 方法。
MATLAB 数组使用列主序，Numpy 数组则主要使用行主序。

一些容易记混的区别：

MATLAB 通过 size() 函数获取形状，而 Numpy 数组则主要通过 ndarray.shape 属性获取形状。
MATLAB 的转置使用'记号（实际是共轭转置），而 Numpy 数组的 ndarray.T 不是共轭转置。

更详细的内容参考官方文档：numpy-for-matlab-users。