CE7453 Numerical Algorithms 期末高分超详细攻略（第五章：Least Squares）

本章内容：最小二乘法、正规方程、QR分解、非线性最小二乘、GPS应用
适用对象：零基础/考前冲刺/快速查漏补缺
内容特色：详细原理、公式推导、算法流程、例题全解、英文关键词

1. 基本概念与背景

• 最小二乘法（Least Squares）：用于拟合数据、解超定方程组（方程数多于未知数），使误差平方和最小。
• 实际应用：数据拟合、回归分析、信号处理、GPS定位、机器学习等。

2. 线性最小二乘法

2.1 问题描述

• 给定 $m$ 个方程 $Ax\approx b$，$A$ 为 $m\times n$ 矩阵（$m>n$），通常无精确解。
• 目标：找到 $x$ 使 $||Ax-b|{|}^2$ 最小。

2.2 正规方程（Normal Equation）

• 推导： $$\underset{x}{\min }||Ax-b|{|}^2=(Ax-b{)}^T(Ax-b)$$ 对 $x$ 求导，令导数为0，得 $$A^{T}Ax=A^{T}b$$
• 解法：解 $n$ 阶方程组 $A^{T}Ax=A^{T}b$

2.3 QR分解法（QR Decomposition）

• 原理：将 $A$ 分解为 $A=QR$，$Q$ 为正交矩阵，$R$ 为上三角矩阵。
• 步骤：
  1. $A=QR$
  2. $Rx=Q^{T}b$，回代求解 $x$
• 优点：数值稳定性高，适合大规模问题

2.4 典型例题

例题1：用最小二乘法拟合直线 $y=ax+b$，已知数据点 $(1,2),(2,3),(3,5)$

详细解答：

1. 问题建模：目标是找到参数 $a$ 和 $b$，使得直线 $y=ax+b$ 尽可能接近给定的数据点。每个数据点对应一个方程：

   • 对于 $(1,2)$：$a\cdot 1+b=2$
   • 对于 $(2,3)$：$a\cdot 2+b=3$
   • 对于 $(3,5)$：$a\cdot 3+b=5$

2. 矩阵形式：将上述方程组写成矩阵形式 $Ax=b$： $$\left[\begin{array}{cc}1& 1\\ 2& 1\\ 3& 1\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}2\\ 3\\ 5\end{array}\right]$$ 这里 $A$ 是 $3\times 2$ 矩阵，$x=[a,b{]}^T$ 是待求参数向量，$b=[2,3,5{]}^T$ 是观测值向量。

3. 正规方程：由于 $A$ 不是方阵，无法直接求逆，我们使用正规方程 $A^{T}Ax=A^{T}b$：

   • 计算 $A^T$： $$A^T=\left[\begin{array}{ccc}1& 2& 3\\ 1& 1& 1\end{array}\right]$$
   • 计算 $A^{T}A$： $$A^{T}A=\left[\begin{array}{ccc}1& 2& 3\\ 1& 1& 1\end{array}\right]\left[\begin{array}{cc}1& 1\\ 2& 1\\ 3& 1\end{array}\right]=\left[\begin{array}{cc}1\cdot 1+2\cdot 2+3\cdot 3& 1\cdot 1+2\cdot 1+3\cdot 1\\ 1\cdot 1+1\cdot 2+1\cdot 3& 1\cdot 1+1\cdot 1+1\cdot 1\end{array}\right]=\left[\begin{array}{cc}14& 6\\ 6& 3\end{array}\right]$$
   • 计算 $A^{T}b$： $$A^{T}b=\left[\begin{array}{ccc}1& 2& 3\\ 1& 1& 1\end{array}\right]\left[\begin{array}{c}2\\ 3\\ 5\end{array}\right]=\left[\begin{array}{c}1\cdot 2+2\cdot 3+3\cdot 5\\ 1\cdot 2+1\cdot 3+1\cdot 5\end{array}\right]=\left[\begin{array}{c}23\\ 10\end{array}\right]$$

4. 求解线性方程组：解 $A^{T}Ax=A^{T}b$，即： $$\left[\begin{array}{cc}14& 6\\ 6& 3\end{array}\right]\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}23\\ 10\end{array}\right]$$ 使用高斯消元法或矩阵求逆：

   • 矩阵行列式：$14\cdot 3-6\cdot 6=42-36=6$
   • 求逆矩阵： $$(A^{T}A{)}^{-1}=\frac{1}{6}\left[\begin{array}{cc}3& -6\\ -6& 14\end{array}\right]=\left[\begin{array}{cc}0.5& -1\\ -1& \frac{14}{6}\end{array}\right]=\left[\begin{array}{cc}0.5& -1\\ -1& \mathrm{2.333...}\end{array}\right]$$
   • 计算 $x=(A^{T}A{)}^{-1}{A}^{T}b$： $$\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{cc}0.5& -1\\ -1& \mathrm{2.333...}\end{array}\right]\left[\begin{array}{c}23\\ 10\end{array}\right]=\left[\begin{array}{c}0.5\cdot 23+(-1)\cdot 10\\ (-1)\cdot 23+\mathrm{2.333...}\cdot 10\end{array}\right]=\left[\begin{array}{c}11.5-10\\ -23+\mathrm{23.333...}\end{array}\right]=\left[\begin{array}{c}1.5\\ \mathrm{0.333...}\end{array}\right]$$ 因此，$a=1.5$，$b=\frac{1}{3}\approx 0.333$。

5. 结果：拟合直线为 $y=1.5x+0.333$。可以通过计算残差 $||Ax-b|{|}^2$ 验证拟合效果。

3. 非线性最小二乘（Nonlinear Least Squares）

3.1 问题描述

• 拟合模型 $y=f(x,\theta )$，$\theta$ 为参数，$f$ 非线性
• 目标：最小化 $S(\theta )={\sum }_{i=1}^m(y_i-f(x_i,\theta ){)}^2$

3.2 Gauss-Newton 方法

• 思想：用泰勒展开线性化，迭代求解
• 步骤：
  1. 初始猜测 ${\theta }_0$
  2. 计算残差 $r_i=y_i-f(x_i,\theta )$
  3. 计算雅可比矩阵 $J_{ij}=\frac{\partial f(x_i,\theta )}{\partial {\theta }_j}$
  4. 解正规方程 $J^{T}J\Delta \theta =J^{T}r$
  5. 更新 $\theta \leftarrow \theta +\Delta \theta$
  6. 迭代至收敛

3.3 典型例题

例题2：用最小二乘法拟合二次曲线 $y=a{x}^2+bx+c$，数据点 $(1,2),(2,3),(3,7),(4,8)$

详细解答：

1. 问题建模：虽然二次曲线是非线性函数，但对于参数 $a,b,c$ 来说，模型是线性的，因此可以用线性最小二乘法直接求解。每个数据点对应一个方程：

   • 对于 $(1,2)$：$a\cdot 1^2+b\cdot 1+c=2$
   • 对于 $(2,3)$：$a\cdot 2^2+b\cdot 2+c=3$
   • 对于 $(3,7)$：$a\cdot 3^2+b\cdot 3+c=7$
   • 对于 $(4,8)$：$a\cdot 4^2+b\cdot 4+c=8$

2. 矩阵形式：将方程组写成矩阵形式 $Ax=b$： $$\left[\begin{array}{ccc}1& 1& 1\\ 4& 2& 1\\ 9& 3& 1\\ 16& 4& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{c}2\\ 3\\ 7\\ 8\end{array}\right]$$

3. 正规方程：计算 $A^{T}A$ 和 $A^{T}b$：

   • $A^T=\left[\begin{array}{cccc}1& 4& 9& 16\\ 1& 2& 3& 4\\ 1& 1& 1& 1\end{array}\right]$
   • $A^{T}A=\left[\begin{array}{ccc}1+16+81+256& 1+8+27+64& 1+4+9+16\\ 1+8+27+64& 1+4+9+16& 1+2+3+4\\ 1+4+9+16& 1+2+3+4& 1+1+1+1\end{array}\right]=\left[\begin{array}{ccc}354& 100& 30\\ 100& 30& 10\\ 30& 10& 4\end{array}\right]$
   • $A^{T}b=\left[\begin{array}{c}1\cdot 2+4\cdot 3+9\cdot 7+16\cdot 8\\ 1\cdot 2+2\cdot 3+3\cdot 7+4\cdot 8\\ 1\cdot 2+1\cdot 3+1\cdot 7+1\cdot 8\end{array}\right]=\left[\begin{array}{c}205\\ 61\\ 20\end{array}\right]$

4. 求解线性方程组：解 $A^{T}Ax=A^{T}b$，即： $$\left[\begin{array}{ccc}354& 100& 30\\ 100& 30& 10\\ 30& 10& 4\end{array}\right]\left[\begin{array}{c}a\\ b\\ c\end{array}\right]=\left[\begin{array}{c}205\\ 61\\ 20\end{array}\right]$$ 使用数值方法或矩阵求逆，得到：

   • $a\approx 0.000$（几乎为0）
   • $b\approx 2.200$
   • $c\approx -0.500$

5. 结果：拟合曲线为 $y=2.2x-0.5$。这表明数据点更适合用直线而非二次曲线进行拟合。可以通过计算残差验证拟合效果。

例题3：用非线性最小二乘法拟合指数模型 $y=a{e}^{bx}$，数据点 $(0,1),(1,2),(2,5)$

详细解答：

1. 问题建模：目标是找到参数 $a$ 和 $b$，使得 $y=a{e}^{bx}$ 拟合数据点。模型为非线性，无法直接用线性最小二乘法求解，因此使用Gauss-Newton迭代方法。

2. 目标函数：最小化残差平方和： $$S(a,b)=\sum _{i=1}^3(y_i-a{e}^{b{x}_i}{)}^2$$

3. 初始猜测：选择初始值 $a=1,b=0.5$（可以根据数据趋势粗略估计）。

4. 迭代步骤：

   • 计算残差：对于当前参数 $a,b$，计算每个数据点的残差 $r_i=y_i-a{e}^{b{x}_i}$。
     • $x=0$：$r_1=1-1\cdot e^{0.5\cdot 0}=1-1=0$
     • $x=1$：$r_2=2-1\cdot e^{0.5\cdot 1}\approx 2-1.6487=0.3513$
     • $x=2$：$r_3=5-1\cdot e^{0.5\cdot 2}\approx 5-2.7183=2.2817$
   • 计算雅可比矩阵：$J_{ij}=\frac{\partial f(x_i,\theta )}{\partial {\theta }_j}$，其中 $f(x,a,b)=a{e}^{bx}$，${\theta }_1=a$，${\theta }_2=b$。
     • $\frac{\partial f}{\partial a}=e^{bx}$
     • $\frac{\partial f}{\partial b}=ax{e}^{bx}$
     • 对于 $x=0$：$J_{11}=e^{0.5\cdot 0}=1$，$J_{12}=0\cdot 1\cdot e^{0.5\cdot 0}=0$
     • 对于 $x=1$：$J_{21}=e^{0.5\cdot 1}\approx 1.6487$，$J_{22}=1\cdot 1\cdot e^{0.5\cdot 1}\approx 1.6487$
     • 对于 $x=2$：$J_{31}=e^{0.5\cdot 2}\approx 2.7183$，$J_{32}=2\cdot 1\cdot e^{0.5\cdot 2}\approx 5.4366$
     • 雅可比矩阵： $$J=\left[\begin{array}{cc}1& 0\\ 1.6487& 1.6487\\ 2.7183& 5.4366\end{array}\right]$$
   • 解正规方程：计算 $J^{T}J$ 和 $J^{T}r$，解 $J^{T}J\Delta \theta =J^{T}r$，得到参数更新量 $\Delta a,\Delta b$。
   • 更新参数：$a\leftarrow a+\Delta a$，$b\leftarrow b+\Delta b$。
   • 重复迭代：重复上述步骤，直到参数收敛或残差足够小。

5. 结果：经过多次迭代，参数收敛到 $a\approx 0.884$，$b\approx 0.864$，拟合模型为 $y\approx 0.884e^{0.864x}$，接近真实指数增长趋势。

4. GPS定位中的最小二乘

• 原理：通过测量与多颗卫星的距离，建立超定方程组，利用最小二乘法估算接收机位置。
• 模型：$||x-s_i||=d_i$，$x$为未知位置，$s_i$为卫星位置，$d_i$为测距
• 步骤：线性化后用Gauss-Newton迭代

4.1 典型例题

例题4：假设接收机与三颗卫星的距离测量如下，估算接收机位置 $(x,y)$：

• 卫星1：位置 $(0,0)$，距离 $d_1=5$
• 卫星2：位置 $(8,0)$，距离 $d_2=5$
• 卫星3：位置 $(4,6)$，距离 $d_3=5$

详细解答：

1. 问题建模：接收机位置为 $(x,y)$，与各卫星的距离方程为：

   • $\sqrt{(x-0{)}^2+(y-0{)}^2}=5⟹x^2+y^2=25$
   • $\sqrt{(x-8{)}^2+(y-0{)}^2}=5⟹(x-8{)}^2+y^2=25$
   • $\sqrt{(x-4{)}^2+(y-6{)}^2}=5⟹(x-4{)}^2+(y-6{)}^2=25$

2. 非线性方程组：上述方程是非线性的，直接求解较复杂，因此使用Gauss-Newton方法进行迭代求解。

3. 初始猜测：根据几何图形，接收机可能位于三颗卫星的中心附近，初始猜测为 $(x,y)=(4,2)$。

4. 线性化与迭代：

   • 定义残差函数：
     • $r_1=\sqrt{x^2+y^2}-5$
     • $r_2=\sqrt{(x-8{)}^2+y^2}-5$
     • $r_3=\sqrt{(x-4{)}^2+(y-6{)}^2}-5$
   • 计算雅可比矩阵（偏导数）：
     • $\frac{\partial r_1}{\partial x}=\frac{x}{\sqrt{x^2+y^2}}$，$\frac{\partial r_1}{\partial y}=\frac{y}{\sqrt{x^2+y^2}}$
     • $\frac{\partial r_2}{\partial x}=\frac{x-8}{\sqrt{(x-8{)}^2+y^2}}$，$\frac{\partial r_2}{\partial y}=\frac{y}{\sqrt{(x-8{)}^2+y^2}}$
     • $\frac{\partial r_3}{\partial x}=\frac{x-4}{\sqrt{(x-4{)}^2+(y-6{)}^2}}$，$\frac{\partial r_3}{\partial y}=\frac{y-6}{\sqrt{(x-4{)}^2+(y-6{)}^2}}$
   • 在初始点 $(4,2)$ 处计算残差和雅可比矩阵，解正规方程更新参数。
   • 经过多次迭代，位置收敛到 $(x,y)\approx (4,1.5)$。

5. 结果：接收机位置为 $(4,1.5)$，可以通过代入原方程验证各卫星到该位置的距离是否接近5。

5. 常见考点与易错点

• 正规方程推导与矩阵乘法
• QR分解步骤与正交化方法（Gram-Schmidt）
• 非线性最小二乘的迭代过程
• 残差（residual）、条件数（condition number）理解
• 英文术语拼写与公式记忆

6. 英文关键词与表达

• least squares, overdetermined system, normal equation, QR decomposition, orthogonalization, Gram-Schmidt, residual, regression, nonlinear least squares, Gauss-Newton, GPS positioning, condition number

如需更详细例题推导或某一算法的代码实现，可随时补充！