Table of Contents

数值计算简介

2023-03-29
8 min read
math

数值计算和数值线性代数息息相关，数值线性代数需要解决的问题主要有三类1：

  1. 解线性方程组：\(Ax=b\)
  2. 线性最小二乘问题：\(\|A x-b\|_2=\min \left\{\|A y-b\|_2: y \in \mathbf{R}^n\right\}\)​
  3. 矩阵特征值问题；奇异值问题

互联网上很多开源工具可以帮助实现数值计算，例如 Boost.Math（统计、数值、求根等）、Boost.Numeric.Odeint（ODE 求解等）。Python 提供了不同类型的矩阵计算库，例如 Numpy、Sympy 等，使用 Python 学习线性代数可以参考这里

误差

数值计算中出现的误差一般有如下几种：

  1. 建立数学模型时做的近似，例如忽略的一些非决定性因素以简化计算
  2. 经验或者计算造成的误差，例如万有引力常数的设置
  3. 工程计算过程造成的误差，例如下面介绍的浮点数误差
    1. 截断误差或者方法误差，例如数值微分常使用泰勒级数的前几级，必然造成一定的损失
    2. 舍入误差，计算机表示浮点数时只能表示有限位

当代绝大部分计算机表示浮点数时都遵循 IEEE754 标准，使用有限的比特位描述无限的浮点数，自然会存在精度损失。IEEE754 的介绍可以参考其他资料，本文只强调一些需要注意的事项

  1. 转换为 10 进制后，单精度浮点数有效小数位最多 7 位；双精度浮点数有效位最多 15 位
  2. 因为精度损失问题，IEEE 754 表示的浮点数，在实数坐标轴上的分布是不均匀的
  3. 为减少浮点数造成的误差，可以考虑如下技巧
    1. 避免中间计算结果出现上溢或下溢 / 避免“大数吃掉小数” / 避免符号相同的两相近数相减
    2. 注意简化步骤，减少运算次数

解方程

求解线性方程组的数值方法大体上可分为直接法和迭代法两大类。直接法是指在没有舍入误差的情况下经过有限次运算可求得方程组的 精确解的方法。因此，直接法又称为精确法。迭代法则是采取逐次逼近的方法，其极限才是方程组的精确解1

LU 分解

LU 分解后可以很方便的利用线性替换法求解方程，故 LU 分解十分重要。求解 LU 分解的方法有多种：

  1. Gauss 消去法（\(O(\frac 2 3 n^3)\)，\(A=LU\)） / 选主元三角分解（列主元消去法：\(PA=LU\)）
    1. 对于条件数较大的矩阵，简单的 Gauss 消去法无法获得正确的结果；可以考虑在消元的过程调整行列的位置，以选择合适的主元，为减少计算量，可以只调整列的位置，故主元三角分解多了一个矩阵 \(P\)
  2. 平方根法又叫 Cholesky 分解（\(O(\frac 1 3 n^3)\)），是求解对称正定线性方程组常用方法之一：\(A=LL^T\)​
    1. \(LDL^T\)​ 分解是 Cholesky 分解的变种，减少了 Cholesky 分解中的开方运算
  3. 分块 LU 分解（参考 BLAS 库）& 其他

因为 LU 分解时间复杂度为 \(O(n^3)\)，大尺寸矩阵使用 LU 方法时的耗时是无法接受的，故一些场景会使用迭代法

迭代法

迭代法是按照某种规则构造一个向量序列 \({x_k}\),使其极限向量 \(x_\infty\) 是方程组 \(Ax=b\) 的精确解1

迭代法的依据是不动点定理2 。以 Jacobi 迭代为例，解方程 \(Ax=b\)，迭代过程为：\(\mathbf{x}^{(k+1)}=D^{-1}\left(\mathbf{b}-(L+U) \mathbf{x}^{(k)}\right)\)

\[A=D+L+U \quad \text { where } \quad D=\left[\begin{array}{cccc} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & a_{n n} \end{array}\right] \text { and } L+U=\left[\begin{array}{cccc} 0 & a_{12} & \cdots & a_{1 n} \\ a_{21} & 0 & \cdots & a_{2 n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n 1} & a_{n 2} & \cdots & 0 \end{array}\right] \]

常用数值算法

最小二乘法

\[\|A x-b\|_2=\min \left\{\|A y-b\|_2: y \in \mathbf{R}^n\right\} \]

最小二乘（LEAST SQUARES REGRESSION）可用于拟合线性函数，也可以用于非线性方程的拟合。Python 中可用的拟合函数有：linalg.lstsq / optimize.curve_fit / numpy.polyval

理解最小二乘法求解可以从代数（解方程）、几何和分析三个角度进行解释与分析。先从代数角度理解最小二乘法，其基本思想是列满秩的矩阵 \(A\) ，\(A^TA\) 是方阵，且存在逆矩阵

\[\begin{aligned} \mathbf{X} \boldsymbol{\beta} & =\mathbf{y} \\ \left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{X} \boldsymbol{\beta} & =\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{y} \\ \boldsymbol{\beta} & =\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{y}\end{aligned} \]

因为数值精度原因，一般使用 \(QR\) 分解来求 \(\beta\)：

\[\begin{aligned} \mathbf{X} \boldsymbol{\beta} & =\mathbf{y} \\ \mathbf{Q R} \boldsymbol{\beta} & =\mathbf{y} \\ \mathbf{R} \boldsymbol{\beta} & =\mathbf{Q}^{\mathrm{T}} \mathbf{y} \\ \boldsymbol{\beta} & =\mathbf{R}^{-1} \mathbf{Q}^{\mathrm{T}} \mathbf{y} \end{aligned} \]

矩阵乘法可以从投影或者线性组合的角度进行思考。\(X\beta\) 意味着将 \(\beta\) 投影到 \(X\) 所在的空间，现实世界中的数据因为噪声的影响，\(X\beta\) 不可能与 \(y\) 完全拟合，即 \(y\) 并不处于 \(X\) 所在的空间。通过正交投影，可以尽可能保证 \(X\beta\) 和 \(y\) 的相似性，两者的误差为 \(\boldsymbol{\epsilon}=\mathbf{X} \boldsymbol{\beta}-\mathbf{y}\)。因为正交投影，可得 \(X^T\epsilon=0\)，分析过程如下（细节可参考这里）

\[\begin{aligned} \mathbf{X}^{\mathrm{T}} \boldsymbol{\epsilon} & =\mathbf{0} \\ \mathbf{X}^{\mathrm{T}}(\mathbf{y}-\mathbf{X} \boldsymbol{\beta}) & =\mathbf{0} \\ \mathbf{X}^{\mathrm{T}} \mathbf{y}-\mathbf{X}^{\mathrm{T}} \mathbf{X} \boldsymbol{\beta} & =\mathbf{0} \\ \mathbf{X}^{\mathrm{T}} \mathbf{X} \boldsymbol{\beta} & =\mathbf{X}^{\mathrm{T}} \mathbf{y} \\ \boldsymbol{\beta} & =\left(\mathbf{X}^{\mathrm{T}} \mathbf{X}\right)^{-1} \mathbf{X}^{\mathrm{T}} \mathbf{y} \end{aligned} \]

这里的 SQUARES REGRESSION 指的是 Squared errors，也就是欧几里得距离 \(\min _{\boldsymbol{\beta}}\|\mathbf{X} \boldsymbol{\beta}-\mathbf{y}\|^2\)，细节可以参考3第十一章

插值法

常见插值（INTERPOLATION）算法：Linear, Cubic Spline, Lagrange, Newton’s polynomial interpolation。一维与二维线性插值算法可以参考这里。如上图所示，一维线性插值计算方法比较直观：\(y=y_0+\left(x-x_0\right) \frac{y_1-y_0}{x_1-x_0}\)。二维线性插值算法需要三次插值，前两次插值用于求解过 P 点直线 \(P_0P_1\)​ 的一维方程，随后利用一维线性插值求解 P 点结果（双线性插值算法可以参考 wiki）

\[\begin{aligned} f\left(x, y_1\right) & =\frac{x_2-x}{x_2-x_1} f\left(Q_{11}\right)+\frac{x-x_1}{x_2-x_1} f\left(Q_{21}\right) \\ f\left(x, y_2\right) & =\frac{x_2-x}{x_2-x_1} f\left(Q_{12}\right)+\frac{x-x_1}{x_2-x_1} f\left(Q_{22}\right) \\ f(x, y) & =\frac{y_2-y}{y_2-y_1} f\left(x, y_1\right)+\frac{y-y_1}{y_2-y_1} f\left(x, y_2\right) \end{aligned} \]

泰勒级数

\[f(x)=\sum_{n=0}^{\infty} \frac{f^{(n)}(a)(x-a)^n}{n !} \]

求根

ROOT FINDING 使用数值算法求解函数的根（ROOT），常见的求根方法有二分法、牛顿法等。BISECTION METHOD 算法4，利用夹挤定理不断逼近根，使用二分法需要提供误差（tolerance）和迭代次数上限

牛顿-拉夫森

牛顿-拉夫森（NEWTON–RAPHSON）方法与梯度下降算法思想类似，不断迭代以接近局部最优解（\(x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\prime}\left(x_n\right)}\)），自然梯度下降有的缺点，牛顿法也有。初始化算法时，要给一个启发值（Guess）

from scipy.optimize import fsolve
f = lambda x: x**3-100*x**2-x+100
fsolve(f, [2, 80])

数值微分

\[\begin{aligned} f\left(x_{j+1}\right) & =\frac{f\left(x_j\right)\left(x_{j+1}-x_j\right)^0}{0 !}+\frac{f^{\prime}\left(x_j\right)\left(x_{j+1}-x_j\right)^1}{1 !}+\frac{f^{\prime \prime}\left(x_j\right)\left(x_{j+1}-x_j\right)^2}{2 !}+\cdots \\ f^{\prime}\left(x_j\right) & =\frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h}+\left(-\frac{f^{\prime \prime}\left(x_j\right) h}{2 !}-\frac{f^{\prime \prime \prime}\left(x_j\right) h^2}{3 !}-\cdots\right) \\ f^{\prime}\left(x_j\right) & =\frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h}+O(h) \\ & \approx \frac{f\left(x_{j+1}\right)-f\left(x_j\right)}{h} \end{aligned} \]

数值微分（NUMERICAL DIFFERENTIATION），用于使用数值方法求解数值函数的局部导数。数值微分的思想可以从不同的角度来理解：

  1. 图形法，数值微分时求解的是波形局部线性化后的切线斜率（forward / backward / central）
  2. 使用泰勒级数前几项导数进行近似

利用泰勒级数差分可以求解函数的多阶导数

\[\begin{aligned} & f\left(x_{j-1}\right)+f\left(x_{j+1}\right)=2 f\left(x_j\right)+h^2 f^{\prime \prime}\left(x_j\right)+\frac{h^4 f^{\prime \prime \prime \prime}\left(x_j\right)}{24}+\cdots \\ & f^{\prime \prime}\left(x_j\right) \approx \frac{f\left(x_{j+1}\right)-2 f\left(x_j\right)+f\left(x_{j-1}\right)}{h^2} \end{aligned} \]

对包含噪声的信号做数值微分，可能得不到有效的结果，因为噪声会让局部变化非常大

数值积分

数值积分使用的思想是黎曼积分，也就是常说的定积分

\[\int_a^b f(x) d x \approx \sum_{i=0}^{n-1} h f\left(x_i\right) \]

利用泰勒级数可以提升黎曼积分精度

\[\begin{aligned} & \int_{x_i}^{x_{i+1}} f(x) d x=\int_{x_i}^{x_{i+1}}\left(f\left(x_i\right)+f^{\prime}\left(x_i\right)\left(x-x_i\right)+\cdots\right) d x \\ & \int_{x_i}^{x_{i+1}} f(x) d x=h f\left(x_i\right)+\frac{h^2}{2} f^{\prime}\left(x_i\right)+O\left(h^3\right) \\ & \int_{x_i}^{x_{i+1}} f(x) d x=h f\left(x_i\right)+O\left(h^2\right) \end{aligned} \]

Trapezoid Rule（梯形法则）和 Simpson’s Rule（辛普森法则）可以用于提升数值积分精度

  1. 梯形法则是一种简单的数值积分方法，它将曲线下的区域划分为一系列梯形，然后计算这些梯形的面积之和来估计定积分的值
  2. 辛普森法则的基本思想是将曲线下的区域近似为一系列二次曲线（或者说是拟合的二次多项式），然后计算这些二次曲线的面积之和

常微分方程（ODEs）

科学技术中很多问题都可以用常微分方程（Ordinary Differential Equation(s)，ODE）的定解问题来描述，主要分为初值（Initial Value）问题和边值（Boundary Value）问题两大类。工程领域一些求不出解析解微分方程，可以用数值方法求出点序列解

初值问题

欧拉方法是常见的数值微分初值问题求解方法，欧拉方法是迭代算法。已知初始值和初始位置（即初始位置的导数）可以近似求解下一个点，随后使用新点求解下一个点。大部分数值微分方法都基于这个思想。求解高阶微分方程时需要将高阶微分方程降阶到一阶微分方程

普通的前向（Forward）微分或者其他微分方式，或多或少都会有一定的精度损失（截断误差，\(o(\Delta t)\)）。为了尽可能精确估算新点的位置，求解微分方程的时候可以使用一些方法与技巧，如梯度法（Trapezoidal Method，\(o(\Delta t^2)\)）、Second-Order Gear（\(o(\Delta t^2)\)）、Runge-Kutta 方法等。下面以一阶微分方程为例描述欧拉后向方法（BDF），微分方程可以整理成如下形式：

\[F\left(x,f(x),{\frac{d f(x)}{d x}},{\frac{d^{2}f(x)}{d x^{2}}},{\frac{d^{3}f(x)}{d x^{3}}},\ldots,{\frac{d^{n-1}f(x)}{d x^{n-1}}}\right)={\frac{d^{n}f(x)}{d x^{n}}} \]

\(n=1\) 时上式可以写为 \(F(x,f(x))=f'(x)\)，初值微分方程问题已知 \(x_0\)，则 \(f'(x_0)\) 已知。以状态方程 \(C\dot x(t)+Gx(t)=Bu(t)\) 为例（\(\dot x(t)=C^{-1}(Bu(t)-Gx(t))\)），系统零状态响应下 \(x(0)=0\)，则可得 \(C\dot x(0)=Bu(0)\)，进而可得 \(\dot x(0)=C^{-1}Bu(0)\)，因为 \(C,B,u\) 三个矩阵已知，故 \(\dot x_0\) 已知。设数值求解过程中的步长为 \(h\)，利用下面迭代公式可解得 \(f(x)\) 的数值解

\[f(x_{n+1})=f(x_n)+h\dot x_n=f(x_n)+hF(x_n,f(x_n)) \]

边值问题

暂略

傅立叶变换

FFT（略）

  1. 徐树方, 高立, and 张平文. 数值线性代数. 北京大学出版社, 2013. ↩︎

  2. 《不动点》 ↩︎

  3. Cohen, Mike X. Practical Linear Algebra for Data Science. " O’Reilly Media, Inc.", 2022. ↩︎

  4. Kong, Qingkai, Timmy Siauw, and Alexandre Bayen. Python programming and numerical methods: A guide for engineers and scientists. Academic Press, 2020. ↩︎

#数值分析 #math #矩阵 #线代