7.2 linalg矩阵运算与求解

7.2 linalg矩阵运算与求解

本节摘要：@ 运算符执行矩阵乘，内部按缓存友好的分块算法执行，是深度学习与统计建模的基石。np.linalg 提供 solve、inv、det、eig、svd 等工具：解方程组应优先 solve 而非先求逆再乘；SVD 是降噪、降维、压缩的万能分解。本节贯穿小案例并给出大规模场景的换库建议。

矩阵乘法与逐元素乘的区别

先分清两个乘法：

import numpy as np

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

print(A * B)    # 逐元素：[[ 5 12] [21 32]]
print(A @ B)    # 矩阵乘：[[19 22] [43 50]]
print(np.dot(A, B))          # @ 的函数形式，结果同上
print(np.matmul(A, B).shape) # matmul 与 @ 等价

# 一维数组在 @ 里被当作向量处理
v = np.array([1, 1])
print(A @ v)    # [3 7]，矩阵乘向量
print(v @ v)    # 2，向量点积（返回标量）
​

矩阵乘为什么值得单独优化？看维度账本：两个 n 阶方阵相乘要做 n 的三次方次乘加，但只读 2 份输入写 1 份输出——计算量远大于数据量，瓶颈在 CPU 缓存与寄存器调度。NumPy 底层调用 BLAS 库，按分块策略把大矩阵切成缓存装得下的小块计算，这也是它对连续性敏感的原因之一：

import numpy as np
import time

n = 1200
A = np.random.rand(n, n)
B = np.random.rand(n, n)

t0 = time.perf_counter()
C1 = A @ B
print("连续矩阵乘:", round(time.perf_counter() - t0, 3), "秒")   # 约 0.1 至 0.3 秒

# 非连续视图参与运算，BLAS 通常会先内部拷贝，变慢
At = np.ascontiguousarray(A.T)
t0 = time.perf_counter()
C2 = At.T @ B
print("转置视图参与:", round(time.perf_counter() - t0, 3), "秒")  # 略慢或相近
print(np.allclose(C1, C2))   # False，注意这里 A.T 就是 A 自身转置的近似演示
​

结论：矩阵乘前不必手工转置优化（BLAS 自己处理），但保持输入连续是好习惯。

solve 与 inv：解方程的正确姿势

线性方程组 Ax = b，教科书教"逆矩阵左乘"，工程上请直接 solve：

import numpy as np

A = np.array([[3.0, 1.0],
              [1.0, 2.0]])
b = np.array([9.0, 8.0])

x = np.linalg.solve(A, b)
print(x)               # [2. 3.]，验算 3*2+1*3=9，1*2+2*3=8

# 对比：先求逆再乘（更慢且数值误差更大）
A_inv = np.linalg.inv(A)
x2 = A_inv @ b
print(x2)              # [2. 3.]

# 顺手的诊断工具
print("行列式:", np.linalg.det(A).round(3))     # 5.0，非零说明可逆
print("条件数:", np.linalg.cond(A).round(2))    # 约 2.6，越小说明数值越稳
​

solve 更快（不显式构造逆）、更稳（带部分选主元的 LU 分解）。条件数是信任度的体温计：条件数到了 10 的十次方量级，解基本不可信，应考虑重新建模而非硬算。

特征分解与 SVD 实战

import numpy as np

M = np.array([[4.0, 1.0],
              [1.0, 3.0]])
vals, vecs = np.linalg.eig(M)
print("特征值:", vals.round(3))    # [4.618 2.382]
# 验证：M @ v = lambda * v
print(np.allclose(M @ vecs[:, 0], vals[0] * vecs[:, 0]))   # True
​

SVD 是更通用的分解，任何形状矩阵都能做，也是降噪压缩的主力：

import numpy as np

rng = np.random.default_rng(4)
data = rng.rand(20, 8)                      # 20 个样本 8 个特征
U, S, Vt = np.linalg.svd(data, full_matrices=False)
print(S.round(3)[:4])   # 奇异值从大到小排列

# 用前 k 个分量重建（秩 k 近似）
k = 3
approx = U[:, :k] * S[:k] @ Vt[:k]
print("重建误差:", np.linalg.norm(data - approx).round(3))
print("原数据能量:", np.linalg.norm(data).round(3))
# 输出示例：重建误差: 2.19，原数据能量: 4.55
​

线代工具速查与规模边界

案例：用最小二乘拟合直线

完整过程。背景：传感器采到 20 个带噪声的 (x, y) 点，拟合 y = kx + b 并评估。

操作：

import numpy as np

rng = np.random.default_rng(8)
x = np.linspace(0, 10, 20)
y = 2.5 * x + 1.0 + rng.normal(0, 1.5, 20)   # 真值斜率2.5 截距1

# 组装正规方程的最小二乘解：A = [x, 1] 列向量堆叠
A = np.column_stack([x, np.ones_like(x)])
k, b = np.linalg.lstsq(A, y, rcond=None)[0]
print("斜率估计:", k.round(3), "截距估计:", b.round(3))
# 输出示例：斜率估计: 2.538 截距估计: 0.812

residual = y - (k * x + b)
print("残差标准差:", residual.std().round(3))   # 约 1.4，接近噪声水平1.5
​

结果与解读：估计值 2.538 与 0.812 逼近真值 2.5 与 1.0，残差标准差回到注入噪声的水平，拟合有效。lstsq 内部就是 SVD，对病态设计矩阵也比正规方程稳。变式：拟合二次曲线只需在 column_stack 里加一列 x 平方，多项式回归的通用模板一模一样；批量拟合多组数据时把 A 一次算好，对每个 y 调 solve 即可。

⚠️ 常见坑：矩阵乘维度不匹配时报的 LinAlgError 或 ValueError 常在深夜出现，先打三个 shape 再读报错——十有八九是最后两个维度没对上（(n,m) @ (m,k) 的 m 必须一致）。

本节要点回顾

• 两种乘法：* 逐元素同形状，@ 矩阵乘末两维对齐；一维数组按向量参与
• solve 优于 inv：更快更稳，inv 只在真的需要逆矩阵本身时用
• 诊断工具：det 判可逆、cond 判数值信任度，超阈值先怀疑建模
• SVD 万能：任意形状可分解，低秩近似实现降噪与压缩，lstsq 是它的工程封装
• 规模边界：稠密千阶内放心用，稀疏换 scipy，GPU 规模换张量框架

最后一章收束全册：合并与分割、文件读写与内存映射，以及一张可以贴在工位的性能清单。