NumPy线性代数
数据分析专题 · 矩阵运算与线性方程组
专题：Python数据分析系统学习
关键词：数据分析, NumPy, 线性代数, linalg, 矩阵乘法, SVD, 特征值, 矩阵分解, 范数
一、概述
NumPy 的 linalg 模块是 Python 生态中最核心的线性代数工具库。它构建在 BLAS/LAPACK 等高性能 Fortran 库之上，提供了从基础矩阵运算到高级矩阵分解的全套函数接口。掌握 np.linalg 的用法，是从事数据科学、机器学习、计算机视觉和科学计算的基础。
线性代数在数据分析中无处不在：主成分分析(PCA)依赖特征值分解或 SVD；线性回归的闭式解需要求解线性方程组；推荐系统中的协同过滤用到 SVD 分解；图神经网络的消息传递本质上是矩阵乘法。理解这些运算不仅要知道"怎么调函数"，更要理解"为什么这样算"。
子模块概览：NumPy 线性代数相关操作分散在三个层面——基础数组运算层（np.dot、np.transpose 等）、np.linalg 标准模块、以及 np.matmul 等运算符层面。建议优先使用 np.linalg 中的统一接口，但也要理解各方式的性能差异与语义区别。
二、矩阵乘法：dot、@ 与 matmul
矩阵乘法是线性代数中最基础也最频繁使用的运算。NumPy 提供了三种等价方式：np.dot、np.matmul 和 @ 运算符（Python 3.5+）。三者在二维数组相乘时结果完全一致，但在高维数组（Tensor）和多维广播行为上存在差异。
2.1 基础矩阵乘法
import numpy as np
# 定义两个二维矩阵
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
# 三种等价写法（二维情况）
C1 = np.dot(A, B)
C2 = np.matmul(A, B)
C3 = A @ B
print(C1)
# [[19 22]
# [43 50]]
print(np.allclose(C1, C2) and np.allclose(C2, C3)) # True
2.2 dot、matmul 与 @ 的区别
| 特性 | np.dot | np.matmul / @ |
|---|
| 二维矩阵 | 标准矩阵乘 | 标准矩阵乘 |
| 一维向量 | 内积（标量） | 不做特殊处理（保持维度） |
| 多维数组（>2维） | 对最后两轴做矩阵乘，广播方式不同 | 将最后两轴视为矩阵，前面维度做广播 |
| 标量乘 | 支持（点乘标量） | 不支持（抛出异常） |
| 推荐场景 | 向后兼容 | 明确的矩阵乘法语义 |
# 一维向量行为对比
x = np.array([1, 2, 3])
y = np.array([4, 5, 6])
print(np.dot(x, y)) # 32（内积，返回标量）
print(np.matmul(x, y)) # 32（内积，返回标量）
print(x @ y) # 32
# 高维数组的广播行为
a = np.ones((2, 3, 4))
b = np.ones((2, 4, 5))
print(np.matmul(a, b).shape) # (2, 3, 5)
print((a @ b).shape) # (2, 3, 5)
# dot 在高维下行为不同
c = np.ones((3, 4))
d = np.ones((4, 5))
print(np.dot(c, d).shape) # (3, 5) —— 二维时相同
最佳实践：在数据分析和机器学习代码中，推荐统一使用 @ 运算符进行矩阵乘法。语法最简洁，语义最清晰，且符合数学记法。仅在需要向后兼容 Python 3.5 以下版本时才使用 np.dot，在高维张量运算中则注意 np.matmul 与 np.dot 的广播差异。
2.3 批量矩阵乘法
# 批量处理：同时计算多组矩阵的乘积
# 形状 (batch, m, n) @ (batch, n, p) -> (batch, m, p)
batched_A = np.random.randn(10, 3, 4) # 10 个 3x4 矩阵
batched_B = np.random.randn(10, 4, 5) # 10 个 4x5 矩阵
result = batched_A @ batched_B
print(result.shape) # (10, 3, 5)
# 利用广播做"单对多"乘法
single = np.random.randn(3, 4)
result2 = single @ batched_B # (3,4) @ (10,4,5) -> (10,3,5)
print(result2.shape) # (10, 3, 5)
三、矩阵运算：转置、逆、行列式、迹与秩
这些是线性代数中最基础的矩阵标量/矩阵属性运算，几乎出现在每一个数据分析实践中。NumPy 提供了高效且易用的接口。
3.1 转置（Transpose）
M = np.array([[1, 2, 3], [4, 5, 6]])
print(M.T)
# [[1 4]
# [2 5]
# [3 6]]
# 高维转置：指定轴顺序
T = np.arange(24).reshape(2, 3, 4)
T_T = T.transpose(2, 0, 1) # 交换轴顺序
print(T_T.shape) # (4, 2, 3)
3.2 逆矩阵（Inverse）与伪逆
A = np.array([[1, 2], [3, 4]])
A_inv = np.linalg.inv(A)
print(A_inv)
# [[-2. 1. ]
# [ 1.5 -0.5]]
# 验证：A @ A_inv = I
print(np.round(A @ A_inv, 10))
# [[1. 0.]
# [0. 1.]]
# 伪逆（Moore-Penrose）：对奇异/非方阵也有效
B = np.array([[1, 2, 3], [4, 5, 6]])
B_pinv = np.linalg.pinv(B)
print(B_pinv.shape) # (3, 2)
print(B @ B_pinv) # 近似单位阵
注意：np.linalg.inv 仅适用于方阵且满秩的情况。若矩阵奇异（行列式为零）或接近奇异，求解会抛出 LinAlgError。在不确定矩阵性质时，优先使用 np.linalg.pinv（伪逆）或使用 np.linalg.solve 求解线性方程组。
3.3 行列式（Determinant）
A = np.array([[1, 2], [3, 4]])
det = np.linalg.det(A)
print(det) # -2.0000000000000004（浮点误差）
# 行列式的几何意义：面积/体积放缩因子
# 二维 → 平行四边形面积；三维 → 平行六面体体积
# det = 0 说明矩阵奇异，线性相关
D = np.array([[1, 2], [2, 4]]) # 线性相关
print(np.linalg.det(D)) # 0.0（浮点误差范围内）
3.4 迹（Trace）与秩（Rank）
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
print(np.trace(A)) # 15（主对角线元素之和）
# 矩阵秩：线性无关的行（或列）的最大数量
print(np.linalg.matrix_rank(A)) # 2（第三行是前两行的线性组合）
# 秩的应用：判断线性系统是否可解
B = np.eye(4)
print(np.linalg.matrix_rank(B)) # 4
矩阵运算总结：转置（.T）改变视图不复制数据；逆（inv）只能用于满秩方阵；伪逆（pinv）通用但计算更慢；行列式（det）用于判断奇异性；迹（trace）等于特征值之和；秩（matrix_rank）基于 SVD 的容差判断，对数值误差鲁棒。
四、矩阵分解
矩阵分解是将一个矩阵拆解为多个结构简单矩阵的乘积。它是求解线性系统、降维、压缩、特征提取的核心技术。NumPy 的 linalg 模块直接支持 QR 分解、SVD 分解、特征值分解和 Cholesky 分解，LU 分解可通过 scipy.linalg 获取。
4.1 LU 分解
LU 分解将矩阵 A 分解为下三角矩阵 L 和上三角矩阵 U 的乘积：A = P L U。其中 P 是置换矩阵。LU 分解是计算机求解线性方程组的标准方法，复杂度 O(n^3)。
from scipy.linalg import lu
import numpy as np
A = np.array([[2, 1, 1], [4, 3, 3], [8, 7, 9]])
P, L, U = lu(A)
print("P:")
print(P)
print("L:")
print(L)
print("U:")
print(U)
# 验证
print(np.allclose(P @ L @ U, A)) # True
4.2 QR 分解
QR 分解将矩阵分解为正交矩阵 Q 和上三角矩阵 R：A = Q R。常用于求解最小二乘问题和特征值算法的基础步骤。
A = np.array([[1, 2], [3, 4], [5, 6]], dtype=float)
Q, R = np.linalg.qr(A)
print("Q (正交):")
print(Q)
print("R (上三角):")
print(R)
# Q 是正交矩阵：Q^T Q = I
print(np.round(Q.T @ Q, 10))
# [[1. 0.]
# [0. 1.]]
# 验证分解
print(np.allclose(Q @ R, A)) # True
4.3 特征值分解
对于方阵 A，特征值分解寻找满足 A v = λ v 的特征值 λ 和特征向量 v。只有可对角化矩阵才能进行完整的特征值分解。
A = np.array([[3, 1], [1, 2]])
eigenvalues, eigenvectors = np.linalg.eig(A)
print("特征值:", eigenvalues)
print("特征向量矩阵:")
print(eigenvectors)
# 验证：A @ v = λ v
v0 = eigenvectors[:, 0]
lam0 = eigenvalues[0]
print(np.allclose(A @ v0, lam0 * v0)) # True
# 对角化验证：A = V diag(λ) V^{-1}
V = eigenvectors
D = np.diag(eigenvalues)
print(np.allclose(V @ D @ np.linalg.inv(V), A)) # True
4.4 SVD 奇异值分解
SVD 是最通用、数值最稳定的矩阵分解方法，适用于任意形状的矩阵（不仅是方阵）。它将矩阵分解为 A = U Σ V^T。SVD 是 PCA、推荐系统、图像压缩、降维等技术的基础。
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]], dtype=float)
U, S, Vt = np.linalg.svd(A)
print("U 形状:", U.shape) # (3, 3)
print("奇异值:", S) # [1.6848e+01 1.0684e+00 3.3352e-16]
print("Vt 形状:", Vt.shape) # (3, 3)
# 验证分解：A = U @ diag(S) @ Vt
Sigma = np.zeros((3, 3))
np.fill_diagonal(Sigma, S)
print(np.allclose(U @ Sigma @ Vt, A)) # True
# 低秩近似（取前k个奇异值）
k = 2
U_k = U[:, :k]
S_k = S[:k]
Vt_k = Vt[:k, :]
A_approx = U_k @ np.diag(S_k) @ Vt_k
print(np.allclose(A_approx, A, atol=1e-1)) # True（近似成立）
SVD 的关键性质：奇异值 S 按从大到小排列，数值越大包含的信息越多。截断 SVD（只保留最大的 k 个奇异值）是降维和去噪的标准做法。相比特征值分解，SVD 的数值稳定性更好，且适用于任意形状矩阵。
4.5 Cholesky 分解
Cholesky 分解是针对对称正定矩阵的"平方根"分解：A = L L^T。它比 LU 分解快约 2 倍，且数值稳定性更好。常用于多元正态采样、卡尔曼滤波等。
# 构造对称正定矩阵
A = np.array([[4, 2], [2, 3]], dtype=float)
# 确保正定
print(np.all(np.linalg.eigvals(A) > 0)) # True
L = np.linalg.cholesky(A)
print("L (下三角):")
print(L)
# 验证：A = L @ L^T
print(np.allclose(L @ L.T, A)) # True
# 应用：从标准正态生成相关正态样本
Z = np.random.randn(2, 1000) # 独立标准正态
X = L @ Z # 相关正态（协方差为 A）
print(np.cov(X).round(1)) # 接近 [[4, 2], [2, 3]]
分解选择指南：求解方阵线性方程组用 LU 或 Cholesky（对称正定时）；最小二乘问题用 QR；数据分析中的降维与可视化用 SVD；可对角化方阵的特征分析用特征值分解。在数值稳定性要求高时，SVD 始终是最安全的选择。
五、线性方程组求解
求解线性系统 A x = b 是科学计算最经典的问题。NumPy 提供了 solve（精确求解方阵系统）和 lstsq（最小二乘求解，适合任意形状）两个核心函数。
5.1 np.linalg.solve（精确求解）
# 求解方程组：
# 2x + 3y = 8
# x - y = -1
A = np.array([[2, 3], [1, -1]])
b = np.array([8, -1])
x = np.linalg.solve(A, b)
print(x) # [1. 2.] → x=1, y=2
# 验证
print(np.allclose(A @ x, b)) # True
重要：solve 要求 A 是方阵且满秩。若 A 奇异或接近奇异，函数会抛出 np.linalg.LinAlgError。在求解之前，建议先用 np.linalg.matrix_rank(A) 检查秩，或直接使用 lstsq。
5.2 np.linalg.lstsq（最小二乘求解）
当方程数多于未知数（超定系统）或矩阵不满秩时，lstsq 返回最小化 ||A x - b||_2 的最小二乘解。这是线性回归的理论基础。
# 超定系统：3 个方程，2 个未知数
A = np.array([[1, 1], [1, 2], [1, 3]])
b = np.array([2, 3, 5])
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)
print("解:", x) # [0.66666667 1.33333333]
print("残差:", residuals) # [0.16666667]
print("秩:", rank) # 2
print("奇异值:", s) # [4.07914333 0.60049122]
# 预测值
y_pred = A @ x
print(y_pred) # [2. 3. 4.] vs 真实值 [2, 3, 5]
print("MSE:", np.mean((b - y_pred)**2)) # 0.0556
5.3 两种方法的选择策略
| 场景 | 推荐方法 | 说明 |
|---|
| 方阵、满秩、需精确解 | solve | 速度最快，直接使用 LU 分解 |
| 任意矩阵、允许近似 | lstsq | 基于 SVD，最通用、最稳定 |
| 对称正定方阵 | cholesky + solve_triangular | 利用 Cholesky 分解加速 |
| 超大规模稀疏系统 | scipy.sparse.linalg | 利用稀疏性大幅降低内存 |
六、特征值与特征向量
特征值与特征向量揭示了线性变换的内在结构。在数据科学中，它们直接关联到主成分分析（PCA）、谱聚类、PageRank 算法、马尔科夫链稳态分析等核心方法。
6.1 通用特征值：eig 与 eigvals
A = np.array([[4, -2], [1, 1]])
# 同时返回特征值和特征向量
eigvals, eigvecs = np.linalg.eig(A)
print("特征值:", eigvals) # [3. 2.]
print("特征向量:\n", eigvecs)
# 仅返回特征值（计算量更小）
vals_only = np.linalg.eigvals(A)
print("仅特征值:", vals_only) # [3. 2.]
# 验证：A @ v = λ v 对每个特征向量成立
for i in range(len(eigvals)):
v = eigvecs[:, i]
lam = eigvals[i]
print(f"λ={lam}: {np.allclose(A @ v, lam * v)}")
6.2 对称矩阵专用：eigh 与 eigvalsh
对于对称（Hermitian）矩阵，eigh 比 eig 快 2-4 倍，且保证返回实数特征值和正交特征向量。PCA 中的协方差矩阵就是对称矩阵，因此应优先使用 eigh。
# 构造对称矩阵（协方差矩阵的典型形式）
A = np.array([[3, 1, 0], [1, 4, 2], [0, 2, 5]], dtype=float)
# 确保对称
A = (A + A.T) / 2
# 使用 eigh（更高效）
eigvals_sym, eigvecs_sym = np.linalg.eigh(A)
print("特征值（eigh）:", eigvals_sym) # 已排序：升序
print("特征向量正交:", np.allclose(eigvecs_sym.T @ eigvecs_sym, np.eye(3)))
# 使用 eigvalsh（仅特征值，更高效）
vals = np.linalg.eigvalsh(A)
print("仅特征值（eigvalsh）:", vals)
# 对比：eig 的耗时（用大矩阵演示）
big = np.random.randn(500, 500)
big = big @ big.T # 对称正定
性能对比：对于 500x500 的对称矩阵，eigh 比 eig 快约 3 倍。在数据分析中，协方差矩阵、Gram 矩阵、拉普拉斯矩阵都是对称的，应始终使用 eigh/eigvalsh。注意 eigh 返回的特征值默认升序排列（最小在前），eig 则是无序的。
6.3 特征值分解的应用：PCA 实现
# 手动实现 PCA 以理解特征值的意义
def pca(X, n_components):
# 中心化
X_centered = X - X.mean(axis=0)
# 协方差矩阵
cov = X_centered.T @ X_centered / (X.shape[0] - 1)
# 特征值分解（对称矩阵用 eigh）
eigvals, eigvecs = np.linalg.eigh(cov)
# eigh 返回升序，取最大的 n_components 个
idx = np.argsort(eigvals)[::-1] # 降序排列
top_idx = idx[:n_components]
# 投影矩阵
W = eigvecs[:, top_idx]
# 降维
X_reduced = X_centered @ W
# 解释方差比
explained_ratio = eigvals[top_idx] / eigvals.sum()
return X_reduced, explained_ratio
# 测试
np.random.seed(42)
X = np.random.randn(100, 5)
X_reduced, ratio = pca(X, 2)
print("降维后形状:", X_reduced.shape) # (100, 2)
print("解释方差比:", ratio)
print("累计解释方差:", ratio.sum())
七、向量范数与矩阵范数
范数是衡量向量或矩阵"大小"的数学工具。在机器学习中，范数用于正则化（L1、L2）、误差度量、梯度裁剪等领域。np.linalg.norm 是统一的范数计算接口。
7.1 向量范数
v = np.array([3, -4])
# L2 范数（欧几里得范数）—— 默认
print(np.linalg.norm(v)) # 5.0 = sqrt(3^2 + (-4)^2)
# L1 范数（曼哈顿范数）
print(np.linalg.norm(v, ord=1)) # 7.0 = |3| + |-4|
# L∞ 范数（切比雪夫范数）
print(np.linalg.norm(v, ord=np.inf)) # 4.0 = max(|3|, |-4|)
# Lp 范数（p=3 为例）
print(np.linalg.norm(v, ord=3)) # (3^3 + 4^3)^(1/3) ≈ 4.4979
# 向量范数常用场景：L2 归一化
v_normed = v / np.linalg.norm(v)
print(v_normed) # [0.6, -0.8]，单位向量
7.2 矩阵范数
A = np.array([[1, 2], [3, -4]])
# Frobenius 范数（默认，逐元素平方和开方）
print(np.linalg.norm(A)) # 5.4772 ≈ sqrt(1+4+9+16)
# 矩阵 2-范数（谱范数）= 最大奇异值
print(np.linalg.norm(A, ord=2)) # 4.8541
# 矩阵 1-范数（最大列绝对值和）
print(np.linalg.norm(A, ord=1)) # 6.0 = max(4, 6)
# 矩阵 ∞-范数（最大行绝对值和）
print(np.linalg.norm(A, ord=np.inf)) # 7.0 = max(3, 7)
# Frobenius 范数与 SVD 关系
U, S, Vt = np.linalg.svd(A)
print(np.sqrt((S**2).sum())) # 5.4772（与 Frobenius 范数相等）
7.3 范数的实际应用
# L2 正则化（Ridge 回归）
from sklearn.linear_model import Ridge
# Ridge 的损失函数: ||y - Xw||^2 + α||w||^2
# L1 正则化（Lasso 回归）
from sklearn.linear_model import Lasso
# Lasso 的损失函数: ||y - Xw||^2 + α||w||_1
# 梯度裁剪
gradients = np.random.randn(100)
threshold = 5.0
grad_norm = np.linalg.norm(gradients)
if grad_norm > threshold:
gradients = gradients / grad_norm * threshold # 缩放到阈值范围
# 相对误差度量
y_true = np.array([1.0, 2.0, 3.0])
y_pred = np.array([1.1, 1.9, 3.2])
relative_error = np.linalg.norm(y_true - y_pred) / np.linalg.norm(y_true)
print(f"相对误差: {relative_error:.4f}") # 约 0.055
范数选择速查：L1 范数用于稀疏性约束和鲁棒性；L2 范数用于能量度量和欧几里得距离；无穷范数用于最大偏差度量；Frobenius 范数是矩阵最自然的"大小"度量。在正则化中，L1 产生稀疏解（特征选择），L2 产生均匀收缩（防止过拟合）。
八、核心要点总结
- 矩阵乘法：
@ 运算符是数据科学中推荐的标准写法；matmul 支持广播；dot 保持向后兼容但高维行为不同
- 逆与伪逆：方阵满秩用
inv；通用情况用 pinv；求解线性系统优先用 solve 或 lstsq
- SVD 是万能工具：数值最稳定、适用任意形状矩阵、直接支持截断近似，是数据降维的首选
- 特征值分解：对称矩阵务必用
eigh（快 2-4 倍），特征值升序排列；非对称用 eig
- 方程求解策略：方阵满秩用
solve（LU 分解）；超定/欠定用 lstsq（SVD 分解）
- 范数的工程意义：L1 正则化产生稀疏模型（特征选择）；L2 正则化控制模型复杂度（防止过拟合）；Frobenius 范数等价于奇异值平方和的开方
- 性能意识：NumPy 的 linalg 底层调用 BLAS/LAPACK，对大规模矩阵计算不要手动实现 Python 循环矩阵乘，应始终利用向量化运算
九、进一步思考
线性代数不仅是工具，更是一种思维方式。理解矩阵分解的本质比记住函数签名重要得多。例如，SVD 分解中左奇异向量 U 的列张成了列空间的一组标准正交基，右奇异向量 V 的列张成了行空间的标准正交基，奇异值 S 则编码了各个方向上的"能量"分布——这种几何直观在理解数据降维时不可或缺。
扩展学习方向：
- scipy.linalg 进阶：提供更多分解（LU、Schur、Hessenberg）、更灵活的矩阵函数（expm、logm、sqrtm、sinm 等矩阵函数）
- 稀疏线性代数（scipy.sparse.linalg）：处理百万级以上规模的稀疏矩阵，使用迭代求解器（CG、GMRES）和稀疏特征值求解（eigs、svds）
- GPU 加速（cupy、cublas）：通过 CuPy 无缝替换 NumPy，利用 GPU 并行计算大幅加速矩阵运算
- 自动微分框架中的线性代数：PyTorch/TensorFlow 的 linalg 模块支持反向传播，使矩阵分解操作可微分，这在现代深度学习中至关重要
实践建议：打开 Jupyter Notebook，用本文的代码示例逐个运行、修改、观察输出变化。尤其建议用真实数据（如 UCI 数据集或 sklearn 内置数据集）跑一遍 PCA 实现和线性回归的闭式解，这样才能真正建立"数学公式 → NumPy 代码 → 数据洞察"的完整链路。
"掌握了 NumPy 的线性代数模块，你就掌握了数据科学的数学语言。每一个函数背后，都是几百年来数学家们智慧的结晶。"