数值方法

数值方法（Numerical Methods）研究如何用有限精度的计算来近似求解数学问题。在实际应用中，许多数学问题（如求解高维积分、偏微分方程、大规模线性系统）没有封闭形式的解析解，必须依赖数值近似。数值方法关注的核心问题是：精度（近似解与精确解的差距）、效率（计算复杂度）和稳定性（误差是否会在计算过程中放大）。

在 AI 领域，数值方法是深度学习优化器、Neural ODE、物理信息神经网络（PINN）等技术的数学基础。本笔记涵盖浮点数表示、插值、数值积分、线性方程组求解、非线性方程求根、ODE 数值解以及数值稳定性分析。

---

1. 浮点数表示与误差

1.1 IEEE 754 浮点数

计算机使用有限位数表示实数，IEEE 754 标准定义了浮点数格式：

\[ x = (-1)^s \times 1.f \times 2^{e - \text{bias}} \]

格式	符号位	指数位	尾数位	bias	精度（十进制）
float32	1	8	23	127	~7 位
float64	1	11	52	1023	~16 位
float16	1	5	10	15	~3 位
bfloat16	1	8	7	127	~2 位

深度学习中的数值精度

• 训练通常使用 混合精度（float16 计算 + float32 累加），可加速 2-3 倍。
• bfloat16 的指数范围与 float32 相同，溢出风险低，被 TPU 和新一代 GPU 广泛采用。
• 量化推理（INT8/INT4）进一步降低精度以换取速度和内存。

1.2 误差类型

误差类型	来源	示例
截断误差	用有限项近似无穷级数	Taylor 展开取前 \(n\) 项
舍入误差	浮点数精度有限	\(0.1 + 0.2 \neq 0.3\)（在 float64 中）
传播误差	误差在计算过程中累积	长序列浮点数求和

绝对误差：\(|x - \hat{x}|\)，相对误差：\(\frac{|x - \hat{x}|}{|x|}\)

---

2. 插值

插值的目标是构造一个函数 \(p(x)\)，使得 \(p(x_i) = y_i\) 对所有已知数据点 \((x_i, y_i)\) 成立。

2.1 拉格朗日插值

给定 \(n+1\) 个数据点，拉格朗日插值多项式为：

\[ L(x) = \sum_{i=0}^{n} y_i \prod_{j \neq i} \frac{x - x_j}{x_i - x_j} \]

• 唯一性：\(n+1\) 个点确定唯一的 \(n\) 次多项式
• 缺点：增加数据点时需要重新计算所有基函数

2.2 牛顿插值

使用差商（divided differences）递推构造：

\[ N(x) = f[x_0] + f[x_0, x_1](x - x_0) + f[x_0, x_1, x_2](x - x_0)(x - x_1) + \cdots \]

其中 \(f[x_i, x_{i+1}] = \frac{f[x_{i+1}] - f[x_i]}{x_{i+1} - x_i}\)

优点：新增数据点时只需计算新的差商项。

2.3 样条插值

高次多项式插值可能产生 Runge 现象（端点处剧烈振荡）。三次样条在每个子区间使用三次多项式，并要求：

• 插值条件：\(S(x_i) = y_i\)
• 连续性：\(S, S', S''\) 在节点处连续
• 边界条件：自然样条（\(S''(x_0) = S''(x_n) = 0\)）或固定斜率

---

3. 数值积分

3.1 梯形法则

将积分区间 \([a, b]\) 分为 \(n\) 等份，用梯形面积近似：

\[ \int_a^b f(x)\,dx \approx \frac{h}{2}\left[f(a) + 2\sum_{i=1}^{n-1}f(x_i) + f(b)\right], \quad h = \frac{b-a}{n} \]

误差为 \(O(h^2)\)。

3.2 辛普森法则

在每两个子区间上使用二次多项式近似：

\[ \int_a^b f(x)\,dx \approx \frac{h}{3}\left[f(a) + 4\sum_{\text{odd } i} f(x_i) + 2\sum_{\text{even } i} f(x_i) + f(b)\right] \]

误差为 \(O(h^4)\)，精度远高于梯形法则。

3.3 高斯求积

选择最优节点（Gauss 点），使得 \(n\) 个点可以精确积分 \(2n-1\) 次多项式：

\[ \int_{-1}^{1} f(x)\,dx \approx \sum_{i=1}^{n} w_i f(x_i) \]

点数	精确积分的最高次数	典型节点
1	1	\(x_1 = 0\)
2	3	\(x_{1,2} = \pm 1/\sqrt{3}\)
3	5	\(x_1 = 0, x_{2,3} = \pm\sqrt{3/5}\)

---

4. 线性方程组

4.1 高斯消元

将增广矩阵 \([A | b]\) 通过行变换化为上三角形，再回代求解：

• 时间复杂度：\(O(n^3)\)
• 部分主元选取（Partial Pivoting）：每步选列中绝对值最大的元素作为主元，减小舍入误差

4.2 LU 分解

将矩阵分解为 \(A = LU\)（下三角 \(\times\) 上三角），等价于将高斯消元的过程"存储"下来：

\[ Ax = b \;\Rightarrow\; Ly = b \;\text{（前向替换）},\; Ux = y \;\text{（回代）} \]

优势：对同一个 \(A\)、不同的 \(b\)，只需一次分解，多次替换。

4.3 迭代法

对于大型稀疏系统，直接法（\(O(n^3)\)）开销过大，迭代法更为实用。

Jacobi 迭代：

\[ x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j \neq i} a_{ij} x_j^{(k)} \right) \]

Gauss-Seidel 迭代：使用最新计算出的分量值，通常收敛更快。

方法	并行性	收敛速度	存储需求
Jacobi	高（可并行）	较慢	需存两组向量
Gauss-Seidel	低（顺序依赖）	较快	原地更新
SOR	低	最快（\(\omega\) 选取合适时）	原地更新

收敛条件：迭代矩阵的谱半径 \(\rho(M) < 1\)。对角优势矩阵保证 Jacobi 和 Gauss-Seidel 收敛。

---

5. 非线性方程求根

求解 \(f(x) = 0\)。

5.1 二分法

• 前提：\(f\) 在 \([a, b]\) 上连续，\(f(a) \cdot f(b) < 0\)
• 每步将区间减半，收敛速度：线性，每步减少 1 个二进制位
• 优点：简单可靠；缺点：收敛慢

5.2 牛顿法

\[ x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \]

• 收敛速度：二次（在单根附近），即误差每步平方缩小
• 需要计算导数 \(f'(x)\)
• 可能不收敛（初始值选取不当或 \(f'(x_n) \approx 0\) 时）

牛顿法在深度学习中的变体

• 一阶优化器（SGD, Adam）可视为"不完整"的牛顿法：只使用梯度，不使用 Hessian。
• 二阶优化器（L-BFGS, K-FAC）近似使用 Hessian 信息，收敛更快但内存开销大。

5.3 割线法

用差商代替导数：

\[ x_{n+1} = x_n - f(x_n) \cdot \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} \]

收敛阶约为 \(\phi \approx 1.618\)（黄金比例），介于线性和二次之间。

---

6. ODE 数值解

求解初值问题 \(\frac{dy}{dt} = f(t, y)\)，\(y(t_0) = y_0\)。

6.1 欧拉法（Euler Method）

\[ y_{n+1} = y_n + h \cdot f(t_n, y_n) \]

• 一阶方法，局部截断误差 \(O(h^2)\)，全局误差 \(O(h)\)
• 简单但精度低，不适合刚性问题

6.2 Runge-Kutta 法（RK4）

经典四阶 Runge-Kutta 方法：

\[ \begin{aligned} k_1 &= f(t_n, y_n) \\ k_2 &= f(t_n + h/2, y_n + h k_1/2) \\ k_3 &= f(t_n + h/2, y_n + h k_2/2) \\ k_4 &= f(t_n + h, y_n + h k_3) \\ y_{n+1} &= y_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4) \end{aligned} \]

局部截断误差 \(O(h^5)\)，全局误差 \(O(h^4)\)。

6.3 刚性问题与隐式方法

刚性问题：ODE 中同时存在快变和慢变分量，显式方法需要极小步长才能保持稳定。

• 隐式欧拉法：\(y_{n+1} = y_n + h \cdot f(t_{n+1}, y_{n+1})\)，A-稳定，适合刚性问题
• BDF（向后差分公式）方法：高阶隐式多步法，广泛用于刚性 ODE 求解器

方法	阶数	稳定性	适用场景
显式欧拉	1	条件稳定	教学演示
RK4	4	条件稳定	非刚性问题
隐式欧拉	1	A-稳定	刚性问题
BDF	1-6	A(\(\alpha\))-稳定	刚性问题

---

7. 数值稳定性

7.1 条件数

矩阵 \(A\) 的条件数衡量输入扰动对输出的影响：

\[ \kappa(A) = \|A\| \cdot \|A^{-1}\| \]

• \(\kappa(A) \approx 1\)：良态问题
• \(\kappa(A) \gg 1\)：病态问题，微小输入变化可能导致巨大输出变化

7.2 向后误差分析

核心思想：数值计算得到的结果 \(\hat{x}\) 不是原问题的精确解，但它是某个"扰动问题"的精确解。

\[ (A + \Delta A)\hat{x} = b + \Delta b \]

如果 \(\|\Delta A\|\) 和 \(\|\Delta b\|\) 足够小（与机器精度同阶），则算法是向后稳定的。

---

8. 在 AI 中的应用

8.1 Neural ODE

Neural ODE 将残差网络的离散层推广为连续的微分方程：

\[ \frac{dh}{dt} = f_\theta(h(t), t) \]

• 使用 ODE 求解器（如自适应 RK45）进行前向传播
• 通过伴随方法（adjoint method）计算梯度，内存复杂度 \(O(1)\)
• 应用：连续正则化流、时间序列建模

8.2 物理信息神经网络（PINN）

PINN 将物理方程作为正则化项嵌入神经网络训练：

\[ \mathcal{L} = \mathcal{L}_{\text{data}} + \lambda \cdot \mathcal{L}_{\text{PDE}} \]

其中 \(\mathcal{L}_{\text{PDE}}\) 利用自动微分计算 PDE 残差。

应用场景：流体力学模拟、材料科学、气象预测等。

8.3 数值线性代数与深度学习

• SVD（奇异值分解）：用于低秩近似（LoRA）、数据降维（PCA）
• 特征值分解：Hessian 谱分析，理解损失曲面
• Krylov 子空间方法：大规模稀疏矩阵的高效求解

---

参考资料

• Burden & Faires, Numerical Analysis
• Trefethen & Bau, Numerical Linear Algebra
• Higham, Accuracy and Stability of Numerical Algorithms
• Chen et al., "Neural Ordinary Differential Equations" (NeurIPS 2018)
• Raissi et al., "Physics-informed neural networks" (JCP 2019)

---

评论