3.1 线性代数运算与分解

3.1 线性代数运算与分解

本节摘要：scipy.linalg 是 SciPy 的稠密线性代数模块，构建在 LAPACK 之上，提供 LU、Cholesky、QR、SVD 四种分解，带假设检查的方程组求解、最小二乘、条件数诊断与 expm、sqrtm 矩阵函数。它与 numpy.linalg 功能重叠但覆盖面更全，本节用可运行示例讲清两者差异、分解选型原则，以及病态矩阵里误差是怎么被条件数放大的。同时给出 lu_factor 与 lu_solve 的分解复用写法、特征值与范数的诊断用法。读完本节，你能独立完成从矩阵结构判断到求解方案落地的完整选型，并解释条件数如何决定误差上限。

读前必看

阅读完本节，你应当能够：

1. 列出 scipy.linalg 相对 numpy.linalg 的主要增量，并解释这些增量解决什么问题；
2. 根据矩阵特性（是否方阵、是否对称正定、是否满秩）选出正确的分解；
3. 用 lu_factor 加 lu_solve 复用分解结果，避免反复求解时的重复计算；
4. 用条件数判断病态程度，并说清误差是怎么被放大的；
5. 用 lstsq 处理超定与秩亏问题，读懂残差、秩与奇异值输出；
6. 用 expm 与 sqrtm 计算矩阵函数，说清与逐元素运算的区别。

一、问题与直觉

先讲一个真实场景。你在做电路瞬态仿真，每个时间步都要解同一个系数矩阵的方程组，变化的只有右端项。起初你顺手用 numpy.linalg.solve，一切正常。等节点数涨到上万，程序明显变慢——你才意识到，每次调用 solve 都把 LU 分解从头做一遍，而分解占了求解九成以上的计算量。你翻遍 numpy.linalg 的文档，发现它没有"分解一次、多次求解"的两步式接口。

这就是 scipy.linalg 存在的意义。它不是 numpy.linalg 的替代品，而是它的完整版：同样基于 LAPACK，但把底层能力几乎全部暴露出来。我们按老规矩，先跑一段最小代码，再逐个拆解。

import numpy as np
from scipy import linalg

A = np.array([[2.0, 1.0, 1.0],
              [1.0, 3.0, 2.0],
              [1.0, 0.0, 0.0]])
b = np.array([4.0, 6.0, 1.0])

x = linalg.solve(A, b)
print(x)                       # 解
print(np.allclose(A @ x, b))   # 验证：A 乘 x 是否等于 b
​

跑通之后你会发现，求解本身和 numpy.linalg.solve 没有区别。区别藏在接口的纵深里：能不能拆两步、能不能声明结构、能不能拿到中间量。这三件事决定了同样一个求解问题，在不同规模、不同结构的矩阵面前，代码和性能会走上完全不同的路。

二、核心原理

2.1 与 numpy.linalg 的分工

两个模块都调用 LAPACK，但暴露的接口深度完全不同。numpy.linalg 是"够用"的子集，追求接口稳定；scipy.linalg 是"全量"封装，把分解、求解、诊断拆得更细，还附带 numpy 里根本没有的矩阵函数。

一句话总结：凡是你需要"分解结果拿来复用"或"针对矩阵结构选择算法"的场景，scipy.linalg 都有对应接口；numpy.linalg 则适合快速调用、不求纵深。

2.2 四种分解：选型先于调用

分解的本质是把一个难算的问题拆成几个好算的问题。选错分解不会报错，只会慢或者不稳，所以选型要先于调用。判断顺序从矩阵形状开始：

LU 分解：适用于一般方阵，带部分主元，是 solve 的默认路径。scipy.linalg.lu 返回置换矩阵 P、下三角 L、上三角 U，满足 A 等于 P 乘 L 乘 U。工程上更常用两步式：lu_factor 做一次分解，lu_solve 反复求解。

P, L, U = linalg.lu(A)
print(np.allclose(P @ L @ U, A))   # 验证分解
lu, piv = linalg.lu_factor(A)      # 分解一次
x = linalg.lu_solve((lu, piv), b)  # 换右端项时只做回代
​

Cholesky 分解：只接受对称正定矩阵，把 A 拆成 L 乘 L 的转置。因为利用了矩阵结构，计算量约为 LU 的一半。用错在非正定矩阵上会直接报错，这反而是好事——它替你做了结构校验。

A_spd = np.array([[4.0, 12.0, -16.0],
                  [12.0, 37.0, -43.0],
                  [-16.0, -43.0, 98.0]])
Lc = linalg.cholesky(A_spd, lower=True)
print(np.allclose(Lc @ Lc.T, A_spd))
cL, low = linalg.cho_factor(A_spd)   # 分解一次
x2 = linalg.cho_solve((cL, low), b)  # 反复求解
​

QR 分解：把 A 拆成正交矩阵 Q 与上三角 R，是非方阵最小二乘的标准路线，数值稳定性优于直接解正规方程。pivoting 参数在秩亏时给出更稳的结果。

SVD 分解：把任意矩阵拆成 U、奇异值对角阵、V 的共轭转置。它不要求方阵、不要求满秩，是最通用的分解，也是条件数、伪逆的诊断基础。代价是四种分解里最贵的，能用 LU 或 QR 解决时不必上 SVD。

U, s, Vh = linalg.svd(A_rect)      # A 等于 U 乘对角 s 乘 Vh
print(linalg.svdvals(A_rect))      # 只取奇异值，不白算 U 和 Vh
​

2.3 solve 的假设检查：把已知结构告诉求解器

scipy.linalg.solve 的 assume_a 参数可以声明矩阵结构：gen 是默认的一般矩阵，sym 是对称，her 是共轭对称，pos 是对称正定。声明之后求解器换用对应专用路径——pos 走 Cholesky，sym 走 LDL 分解。对一万阶对称正定矩阵，这个参数能带来接近一倍的加速，而且数值上更稳。

x_pos = linalg.solve(A_spd, b, assume_a='pos')   # 对称正定路径
x_sym = linalg.solve(A_sym, b, assume_a='sym')   # 对称路径
​

另一个参数 check_finite 默认检查输入是否含非有限值。数据可靠、追求极限性能时可以关掉，省一次全矩阵扫描，但代价是遇到 NaN 时静默出错。我的建议：调试期别关，上线后确认数据源干净再关。

声明结构还有个副作用：它替你做了一半的输入校验。assume_a 声明与矩阵实际结构不符时，求解器多半会报错而不是给出错误解——比如对非正定矩阵声明 pos，Cholesky 路径会在某个对角元上崩掉。这比默默算出垃圾结果好得多，也算是一种廉价的防御式编程。

2.4 条件数与病态矩阵：残差会骗人

求解误差不只看算法，更看矩阵本身。条件数等于最大奇异值除以最小奇异值，衡量矩阵对误差的放大倍数：解的相对误差大约等于条件数乘上数据误差与机器精度之和。条件数接近 1 是良态，超过 1e12 基本不可信。

经典的病态例子是希尔伯特矩阵，元素是行号加列号减一的倒数：

H = linalg.hilbert(10)
b_h = H @ np.ones(10)
x_h = linalg.solve(H, b_h)
print(np.max(np.abs(x_h - 1.0)))   # 真实解是全是 1 的向量
print(linalg.cond(H))              # 十阶希尔伯特的条件数
​

结果会让你意外：A 乘 x_h 与 b_h 的残差很小，但解和真解差出好几个数量级。这就是病态矩阵的典型特征——残差小不代表解准。反过来说，看到"残差正常但结果离谱"的报告，第一反应就该是查条件数。处理手段按顺序是：先检查建模是否合理（网格过密、单位悬殊都会放大病态），再考虑换基或预处理，最后才轮到换求解器。

2.5 矩阵函数：expm 与 sqrtm

矩阵函数把标量函数推广到矩阵上，是常微分方程数值解的地基。比如线性系统状态对时间的导数等于系数矩阵乘状态，解析解就是矩阵指数乘初始状态。这里的指数必须用 linalg.expm，它内部做缩放平方加 Padé 逼近；直接写 np.exp 是逐元素指数，完全是另一回事，算出来的矩阵通常没有任何物理意义。

u_t = linalg.expm(A_drift * dt) @ u0   # 常微分方程的解
S = linalg.sqrtm(A_spd)                # 矩阵平方根，协方差变换常用
​

sqrtm 求矩阵平方根，对对称正定矩阵有明确意义；logm 是其逆运算；funm 可以套任意标量函数。需要提醒的是，矩阵函数计算成本不低，别对巨型矩阵直接用，先确认规模。

2.6 特征值、行列式与范数

这一节把剩下的常用操作收个尾。特征值用 eig，返回的特征向量按列排布；如果矩阵对称或 Hermite，改用 eigh，它利用结构更快也更稳，特征值排序也有保证。行列式 det 和范数 norm 都是诊断工具：det 接近零提示奇异，但尺度敏感的矩阵里 det 会骗人——一个良态但数值很大的矩阵 det 也可以很大，所以判断奇异性优先看条件数而不是行列式。norm 支持多种范数，默认的 Frobenius 范数等于所有元素平方和的平方根，2 范数等于最大奇异值，两者在误差分析里都常用。这些函数在 numpy.linalg 里也有对应版本，结果一致，按前文的建议统一用 scipy.linalg 即可。

三、工程实践要点

3.1 分解选型速查

3.2 常见问题

问：什么时候用 numpy.linalg 就够了？
答：只做一次性求解、特征值、SVD，且不需要分解复用与假设参数时，两边结果一致，用哪个都行。我倾向在科学计算脚本里统一用 scipy.linalg，少记一套接口。

问：cond 算出来很大，接下来怎么办？
答：先检查建模是否合理，比如网格是否过密、单位是否悬殊；再考虑换基或预处理；最后才轮到换求解器。条件数是症状，不是病因。

问：方程组无解或者解不唯一时，solve 会怎样？
答：方阵且奇异时 solve 会报奇异警告或者给出充满 NaN 的解；方程组不相容或未知数多于方程时，应该改用 lstsq，它返回最小二乘意义下的解、残差、秩和奇异值，让你看清问题到底出在哪一步。

⚠️ 常见坑：用逆矩阵求解，即 x 等于 inv 乘 b。求逆本身要花一次分解的钱，数值稳定性还更差；正确做法是直接 solve，或分解后回代。
💡 关键直觉：求解前先 cond 一下，一行代码，能省掉大量排查"结果为什么不对"的时间。条件数是线性代数问题里最便宜的体检。

要点速记

• 要点一：scipy.linalg 是 numpy.linalg 的超集，核心增量是分解复用接口、assume_a 假设参数与矩阵函数。
• 要点二：LU 是一般方阵求解的默认路线，lu_factor 加 lu_solve 把分解缓存，换右端项只做回代。
• 要点三：Cholesky 只收对称正定矩阵，计算量约为 LU 的一半，用错矩阵会直接报错。
• 要点四：QR 与 SVD 是非方阵最小二乘的主武器，SVD 额外给出秩与奇异值诊断。
• 要点五：条件数决定误差放大倍数，希尔伯特矩阵这类病态矩阵里，残差小不代表解准。
• 要点六：expm 与 sqrtm 是矩阵函数，与逐元素 np.exp 完全不同，常微分方程求解依赖它们。
• 要点七：工程顺序是 cond 诊断在前、选分解在中、验证残差在后。

下一节我们把规模推高一个数量级：当稠密矩阵装不进内存时，scipy.sparse 的三种存储格式与迭代求解器登场。

图：矩阵分解选择示意