4. 线性代数基础选型

在 SciPy 生态里做线性代数，最先遇到的困惑往往不是“公式怎么写”，而是“到底该调哪个函数”。NumPy 提供了 numpy.linalg，SciPy 又提供了 scipy.linalg，两者接口高度相似，都能解方程、求逆、做最小二乘。选错的代价通常不是报错，而是精度下降、速度变慢，或者在稀疏问题上直接把内存撑爆。本节把选型拆成三个连续判断：先看库，再看矩阵规模与稀疏性，最后看任务本身需要的是解、逆还是最小二乘解。

一、两个库的真实差异

numpy.linalg 是 NumPy 自带的轻量实现，依赖底层 LAPACK/BLAS，覆盖面以常用分解和求解为主，优势是“无需额外依赖、随 NumPy 一起可用”。scipy.linalg 则是对 LAPACK 更完整的封装，函数更多、参数更细，并且提供了一批 NumPy 没有的能力：lu_factor/lu_solve 的分解复用、solve_triangular 的三角求解、lstsq 的多种驱动选择（gelsd、gelss、gelsy）、cho_factor/cho_solve 的 Cholesky 复用，以及 scipy.linalg.blas、lapack 的底层直调。对稠密中小矩阵，两者结果通常一致；差别体现在可控性、复用性和边界场景。

一个常被忽略的点是：scipy.linalg 默认会做更积极的输入检查与覆盖写控制（overwrite_a、check_finite），在循环内反复求解时，这些开关会明显影响性能。因此“用哪个库”不是信仰问题，而是看是否需要分解复用、是否需要精细控制。

二、先判断稀疏性，再谈函数

选型的第一道分水岭是稀疏性。若矩阵中非零元素占比很低（例如有限差分、图拉普拉斯、网络关联矩阵），用 numpy.linalg 或 scipy.linalg 的稠密例程会把大量零也存成浮点数，内存按 O(n²) 增长。此时应转向 scipy.sparse 体系：scipy.sparse.linalg.spsolve 做直接稀疏求解，scipy.sparse.linalg.cg、gmres 等做迭代求解，scipy.sparse.linalg.lsqr 处理稀疏最小二乘。判断阈值没有绝对标准，工程上常用“非零率低于 5%~10% 且规模上千”作为转向稀疏的参考信号，但更可靠的做法是实测内存与耗时。

稠密矩阵内部还要看规模。小矩阵（几十到几百阶）用哪个库差别不大，优先选调用更顺手的；中等规模（数百到数千阶）且需要多次求解同一系数矩阵时，分解复用带来的收益开始显现；大规模稠密矩阵则要关注是否值得改用迭代法或利用结构（对称、正定、带状、Toeplitz）。

三、solve、inv、lstsq 的取舍逻辑

solve 是默认答案。只要目标是“已知 A 和 b，求 x 使 Ax=b”，就应调用 solve，而不是先 inv(A) 再乘 b。原因很直接：求逆需要 O(n³) 的分解加回代，再乘 b 又增加一次矩阵乘法，而 solve 只做一次分解加回代，运算量和舍入误差都更小。数值线性代数的经验法则是“能不解逆就不解逆”。只有在确实需要逆矩阵本身参与后续运算（例如构造投影矩阵、计算协方差、做 Sherman-Morrison 更新）时，才用 inv。

inv 的合理场景要说得清。如果后续要对同一个 A 反复求解大量不同的 b，正确做法也不是先求逆，而是做一次分解并复用：scipy.linalg.lu_factor 得到 LU 因子，再用 lu_solve 对每个 b 回代；若 A 对称正定，用 cho_factor/cho_solve 更快更稳。这比 inv 再乘更精确，也比每次重新 solve 更省。

lstsq 处理非方阵与病态。当 A 是 m×n 且 m≠n，或者 A 接近奇异、方程无精确解时，应使用 lstsq 求最小二乘解。它返回解、残差平方和、秩和奇异值，后两者是判断结果可信度的关键：秩小于 n 说明解不唯一，奇异值跨度大说明病态。scipy.linalg.lstsq 还允许通过 lapack_driver 选择 gelsd（默认，基于 SVD，稳健）、gelss、gelsy（基于 QR，较快）。对病态问题，gelsd 配合 cond 截断小奇异值通常更可靠。

四、一段可直接运行的对照示例

下面用同一个 3×3 方程组对比三种写法，并给出一个超定最小二乘的例子。示例参数为教学构造，结果以实际运行输出为准。

import numpy as np
import scipy.linalg as sla

A = np.array([[3.0, 1.0, 0.0],
              [1.0, 4.0, 1.0],
              [0.0, 1.0, 5.0]])
b = np.array([9.0, 16.0, 21.0])

# 1) 推荐：直接求解
x_solve = sla.solve(A, b)

# 2) 不推荐但可对照：先求逆再乘
x_inv = sla.inv(A) @ b

# 3) 复用分解：同一 A 多次求解
lu, piv = sla.lu_factor(A)
x_lu = sla.lu_solve((lu, piv), b)

print(np.allclose(x_solve, x_lu))   # True
print(np.linalg.norm(x_solve - x_inv))  # 通常为 1e-15 量级

# 4) 超定最小二乘：4 个方程、2 个未知量
A_over = np.array([[1.0, 0.0],
                   [1.0, 1.0],
                   [1.0, 2.0],
                   [1.0, 3.0]])
y = np.array([1.0, 3.0, 5.0, 7.0])
sol, res, rank, sv = sla.lstsq(A_over, y)
print(sol, res, rank, sv)

运行后可以观察到：x_solve 与 x_lu 在浮点容差内一致，说明分解复用与直接求解等价；x_inv 的差异通常在 1e-15 量级，矩阵越大、条件数越差，这个差异越明显。最小二乘例子中，若数据恰好落在一条直线上，残差平方和接近 0，秩为 2，奇异值给出两个正数，说明解唯一且拟合良好。若把 y 改成带噪声的值，残差会变为正数，这就是判断“模型是否解释得了数据”的直接依据。

五、选型检查清单

• 矩阵非零率低、规模大：先考虑 scipy.sparse.linalg，不要直接上稠密例程。
• 方阵、非奇异、只求一个解：solve，不要 inv 再乘。
• 同一 A 反复求解：lu_factor/lu_solve，对称正定用 cho_factor/cho_solve。
• 非方阵、秩亏、病态：lstsq，并检查返回的 rank 与奇异值。
• 需要更细的 LAPACK 控制或更多分解类型：优先 scipy.linalg；只求随 NumPy 可用的基础能力：numpy.linalg 足够。