Badge

本文目标:让没有任何数值分析基础的读者,从零开始系统掌握数值分析的核心概念——浮点数表示、误差分析、数值微分与积分、常微分方程数值解法(欧拉法/RK4)、线性方程组求解、以及TensorRT量化的数学原理。全文采用”概念→定义→性质→历史→应用”教学法,所有内容连接到飞行动力学仿真、光线追踪和深度学习推理。


一、什么是数值分析——数学公式如何在计算机上”算得准”

1.1 数值分析回答的核心问题

理论上完美的数学公式,在计算机上计算时会遇到三个问题:

问题 例子 后果
浮点误差 1+1016=11 + 10^{-16} = 1(单精度) 累计误差导致发散
病态问题 条件数大的矩阵 Ax=bAx=b 很小的输入噪声导致很大的输出误差
截断误差 欧拉方法近似微分方程 仿真结果与真实值有偏差

数值分析就是研究如何在计算机上高效、稳定、精确地计算数学问题

1.2 发展简史

历史

  • 1671年:牛顿发明牛顿法求解方程,是最早的数值算法之一
  • 1805年:高斯发明最小二乘法
  • 1901年:龙格(Runge)和库塔(Kutta)发明RK4方法
  • 1947年:冯·诺依曼研究数值误差传播,标志数值分析成为独立学科
  • 1960年代:威尔金森(Wilkinson)出版《代数过程中的舍入误差》,奠定理论基础
  • 2010年代:深度学习兴起,数值精度与模型压缩成为热点

二、浮点数——计算机如何表示实数

2.1 浮点数的二进制表示

计算机用浮点数近似表示实数。IEEE 754标准定义了浮点数格式:

value=(1)s×(1+f)×2ebias \text{value} = (-1)^s \times (1 + f) \times 2^{e - \text{bias}}

其中:

  • ss:符号位(0正1负)
  • ee:指数部分
  • ff:尾数部分(小数)

2.2 常用浮点精度

格式 总位数 指数位 尾数位 有效数字 范围
FP32(单精度) 32 8 23 约7位 ±3.4×1038\pm 3.4 \times 10^{38}
FP64(双精度) 64 11 52 约16位 ±1.8×10308\pm 1.8 \times 10^{308}
FP16(半精度) 16 5 10 约3位 ±6.5×104\pm 6.5 \times 10^{4}
BF16(脑浮点) 16 8 7 约3位 ±3.4×1038\pm 3.4 \times 10^{38}
INT8(整数量化) 8 0位 128127-128 \sim 127

2.3 机器精度

机器精度 ϵm\epsilon_m:1与大于1的最小可表示浮点数之间的差值。

精度 ϵm\epsilon_m
FP32 1.19×1071.19 \times 10^{-7}
FP16 9.77×1049.77 \times 10^{-4}
BF16 7.81×1037.81 \times 10^{-3}

重要规则

  • 如果 aba \ll ba+bba + b \approx b(小数字被”吃掉”)
  • 两个接近的数相减会灾难性抵消(相对误差急剧放大)
  • 浮点数加法不满足结合律:(a+b)+ca+(b+c)(a+b)+c \neq a+(b+c)

2.4 TensorRT量化——从FP32到INT8

为什么量化? 将FP32权重压缩到INT8:模型大小缩小4倍,推理速度提升2-4倍。

量化公式

xINT8=round(xFP32/s)+z x_{\text{INT8}} = \text{round}(x_{\text{FP32}} / s) + z

其中 ss 是缩放因子,zz 是零点偏移。

校准方法——如何确定最佳缩放因子:

方法 原理 精度损失
MinMax 记录激活值的最小/最大值 对异常值敏感
Entropy 最小化量化前后的KL散度 最常用,精度损失最小
Percentile 去除异常值(如99.9%) 适合有离群点的分布

在YOLO中的应用:YOLOv8使用TensorRT INT8量化后,mAP下降通常不超过1%,但推理速度提升约2.7倍。


三、误差分析

3.1 误差的分类

类型 来源 特点
舍入误差 浮点数精度有限 每一步微小,千百万步后可能显著
截断误差 用有限步近似无限过程 步长越小,误差越小
测量误差 传感器精度有限 与外因有关,无法消除
模型误差 数学模型不完美 理论建模阶段引入

3.2 条件数——问题本身的”坏”程度

定义:条件数衡量一个问题的解对输入扰动的敏感程度。

对于线性方程组 Ax=bA\mathbf{x} = \mathbf{b},条件数为:

κ(A)=AA1=σmaxσmin \kappa(A) = \|A\| \cdot \|A^{-1}\| = \frac{\sigma_{\max}}{\sigma_{\min}}

解释

  • κ(A)1\kappa(A) \approx 1:良态问题,输入的小扰动只引起解的小变化
  • κ(A)1\kappa(A) \gg 1:病态问题,输入的小扰动可能导致解的巨大偏差
  • κ(A)\kappa(A) \to \infty:矩阵接近奇异,求逆结果不可信

在飞行动力学中的应用:控制分配矩阵 MM 的条件数决定了四旋翼能否精确实现期望力矩。如果某个方向的条件数极大,说明该方向难以控制。

历史:条件数的概念由英国数学家图灵(Alan Turing)在1948年引入。

3.3 数值稳定性

定义:一个算法如果对输入的小误差不敏感,就是数值稳定的。

前向误差 vs 后向误差

  • 前向误差y^y\|\hat{y} - y\|(计算解与真实解的差异)
  • 后向误差min{δ:y^=f(x+δ)}\min\{\delta : \hat{y} = f(x + \delta)\}(需要多大的输入扰动才能得到计算解)

后向稳定算法:后向误差很小(O(ϵm)O(\epsilon_m))的算法。

案例:求解 Ax=bAx = b 时,高斯消去(带部分选主元)是后向稳定的,但Cramer法则(用行列式)极不稳定。


四、数值微分与数值积分

4.1 数值微分——当导数公式未知时

前向差分

f(x)f(x+h)f(x)h f'(x) \approx \frac{f(x+h) - f(x)}{h}

中心差分(更精确):

f(x)f(x+h)f(xh)2h f'(x) \approx \frac{f(x+h) - f(x-h)}{2h}

误差分析:中心差分的截断误差为 O(h2)O(h^2),前向差分为 O(h)O(h)

选择步长 hh 的权衡

  • hh 太大 → 截断误差大
  • hh 太小 → 舍入误差大(f(x+h)f(x+h)f(x)f(x) 几乎相等,相减导致灾难性抵消)

最优 hh 通常在 ϵmx\sqrt{\epsilon_m} \cdot |x| 附近。

在深度学习中的应用:反向传播就是数值微分的精确替代——它用微积分中的链式法则计算精确梯度,而不是用数值近似。

4.2 数值积分——当原函数未知时

梯形法则

abf(x)dxh2[f(a)+2i=1n1f(xi)+f(b)] \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]

Simpson法则(更精确):

abf(x)dxh3[f(a)+4i=1n/2f(x2i)+2i=1n/21f(x2i+1)+f(b)] \int_a^b f(x)dx \approx \frac{h}{3}\left[f(a) + 4\sum_{i=1}^{n/2}f(x_{2i}) + 2\sum_{i=1}^{n/2-1}f(x_{2i+1}) + f(b)\right]

误差:Simpson法则的截断误差为 O(h4)O(h^4),梯形的为 O(h2)O(h^2)


五、常微分方程数值解法——仿真计算的基石

5.1 欧拉法——最简单的数值解法

问题:求解 dydt=f(t,y)\frac{dy}{dt} = f(t, y),已知 y(t0)=y0y(t_0) = y_0

欧拉法

yn+1=yn+hf(tn,yn) y_{n+1} = y_n + h f(t_n, y_n)

其中 hh 是时间步长。

截断误差O(h)O(h)(一阶精度)

例子:假设加速度恒为 9.8m/s29.8\text{m/s}^2Δt=0.1s\Delta t = 0.1\text{s},初速度 v0=0v_0=0

时间 速度
0 0 0
1 0.1 0+0.1×9.8=0.980 + 0.1 \times 9.8 = 0.98
2 0.2 0.98+0.1×9.8=1.960.98 + 0.1 \times 9.8 = 1.96
3 0.3 1.96+0.1×9.8=2.941.96 + 0.1 \times 9.8 = 2.94

精确值:v=at=9.8×0.3=2.94v = at = 9.8 \times 0.3 = 2.94——完全精确!但对于更复杂的微分方程,欧拉法会有累积误差。

5.2 四阶龙格-库塔法(RK4)

RK4是工程中最常用的ODE求解器,精度远高于欧拉法。

算法

k1=f(tn,yn)k2=f(tn+h/2,yn+hk1/2)k3=f(tn+h/2,yn+hk2/2)k4=f(tn+h,yn+hk3)yn+1=yn+h6(k1+2k2+2k3+k4) \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(h4)O(h^4)(四阶精度)

为什么RK4优于欧拉法? RK4在每一步内取了4个点的导数,加权平均后作为最终增量。这相当于用更高阶的多项式拟合真实解。

在飞行动力学中的应用:PX4仿真和Gazebo中求解六自由度运动方程时,使用的就是RK4或更先进的变步长求解器。

5.3 步长选择

方法 精度 每步计算量 适用场景
欧拉法 O(h)O(h) 1次 快速原型、定性分析
RK4 O(h4)O(h^4) 4次 工程通用标准
变步长RK 自适应 约6次 刚性问题、高精度需求

六、线性方程组的数值求解

6.1 高斯消去法——最直接的方法

步骤:将增广矩阵 [Ab][A|\mathbf{b}] 通过行变换变为上三角矩阵,然后回代求解。

计算量23n3\frac{2}{3}n^3 次浮点运算。

需要部分选主元(Partial Pivoting):在每一步选择当前列绝对值最大的元素作为主元,避免除以接近0的数。

6.2 LU分解——更高效的方法

A=LU A = LU

其中 LL 是下三角矩阵,UU 是上三角矩阵。

优点:一旦分解完成,对不同的右端项 b\mathbf{b} 求解只需 O(n2)O(n^2) 次运算(回代),而不是重新做 O(n3)O(n^3) 的高斯消去。

6.3 Cholesky分解——对称正定矩阵的专用方法

A=LLT A = LL^T

条件AA 必须是对称正定矩阵。

优点:计算量约为LU分解的一半,数值稳定性更好。

在无人机中的应用:卡尔曼滤波中的协方差矩阵更新就使用了Cholesky分解。


七、TensorRT量化——数值分析在AI中的应用

7.1 量化流程

1
FP32模型 → 校准数据集 → 确定缩放因子 → INT8模型 → 推理

7.2 量化对精度的影响

以YOLOv8m在RTX 4090上的测试为例:

精度 延迟 mAP@0.5:0.95 加速比
FP32 5.2ms 45.4% 1.0×
FP16 2.8ms 45.3% 1.9×
INT8 1.9ms 44.5% 2.7×

关键观察:INT8量化损失约0.9% mAP,但速度提升2.7倍——在边缘部署(Jetson Orin)中是一个值得的权衡。

7.3 量化误差的来源

  • 裁剪误差:超出INT8表示范围的激活值被裁剪
  • 舍入误差:FP32到INT8的精度损失
  • 校准误差:校准集不能完全代表推理时的分布

八、核心公式速查卡

公式 含义 应用
ϵm=1.19×107\epsilon_m = 1.19 \times 10^{-7} (FP32) 机器精度 浮点计算误差边界
κ(A)=AA1=σmax/σmin\kappa(A) = \|A\|\|A^{-1}\| = \sigma_{\max}/\sigma_{\min} 条件数 矩阵病态程度、控制分配稳定性
f(x)f(x+h)f(xh)2hf'(x) \approx \frac{f(x+h)-f(x-h)}{2h} 中心差分 数值微分
yn+1=yn+h6(k1+2k2+2k3+k4)y_{n+1} = y_n + \frac{h}{6}(k_1+2k_2+2k_3+k_4) RK4 飞行动力学仿真
A=LUA = LU LU分解 线性方程组求解
xINT8=round(xFP32/s)+zx_{\text{INT8}} = \text{round}(x_{\text{FP32}}/s) + z 量化 TensorRT推理加速

参考文献与推荐资源

  1. Trefethen, L. N., & Bau, D. (1997). Numerical Linear Algebra. SIAM. — 数值线性代数经典
  2. Burden, R. L., & Faires, J. D. (2010). Numerical Analysis. 9th ed., Brooks/Cole. — 最全面的数值分析教材
  3. Press, W. H., et al. (2007). Numerical Recipes. 3rd ed., Cambridge University Press. — 数值算法实用指南
  4. Higham, N. J. (2002). Accuracy and Stability of Numerical Algorithms. 2nd ed., SIAM. — 数值稳定性权威著作
  5. Wilkinson, J. H. (1963). Rounding Errors in Algebraic Processes. — 数值误差分析奠基之作
  6. Stevens, B. L., et al. (2015). Aircraft Control and Simulation. 3rd ed., Wiley. — RK4在飞行动力学中的应用