4.3 最小二乘近似

passxgx

于 2024-07-18 22:52:15 发布

阅读量788

点赞数 24

分类专栏： # 第4章正交文章标签：线性代数

本文链接：https://blog.csdn.net/passxgx/article/details/140363245

版权

第4章正交专栏收录该内容

4 篇文章 0 订阅

订阅专栏

一、最小二乘解

$A\boldsymbol x=\boldsymbol b$ 经常无解，一般是因为方程太多了。矩阵 $A$ 的行比列要多，即方程要多余未知数（ $m > n$ ）。 $n$ 个列只能张成 $m$ 空间的一小部分，除非测量的值都非常完美，否则 $\boldsymbol b$ 将会落在 $A$ 的列空间之外。此时使用消元法将会得到不可能存在的方程，因为测量值包含噪声，我们要找到最优解。
重复：我们不可能能一直把误差 $\boldsymbol e=\boldsymbol b-A\boldsymbol x$ 降为零，当 $\boldsymbol e$ 为零时， $\boldsymbol x$ 是 $A\boldsymbol x=\boldsymbol b$ 的确切解；否则我们要使得 $\boldsymbol e$ 的长度尽可能的小，此时 $\boldsymbol {\hat x}$ 就是最小二乘解（least squares solution）。本节的目的就是计算并使用 $\boldsymbol {\hat x}$ ，因为有许多实际问题需要答案。
以前强调 $\boldsymbol p$ （投影），现在强调 $\boldsymbol {\hat x}$ （最小二乘解），它们由 $\boldsymbol p=A\boldsymbol {\hat x}$ 联系起来。基本的方程任然是 $A^TA\boldsymbol {\hat x}=A^T\boldsymbol b$ ，这里有一个非正式的方法得到这个 “正态方程”：

$当\,A\boldsymbol x=\boldsymbol b\,没有解时，两边同时左乘\,A^T\,然后求解\,\colorbox{white}{$A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$}$

【例1】最小二乘的一个重要应用就是对于 $m$ 个点拟合成一条直线。从三个点开始：求出最接近点 $(0, 6), (1, 0)$ 和 $(2, 0)$ 的直线。
没有任何一条直线可以同时过着三个点，假设直线方程为 $b = C + D t$ ，我们要求出两个数字 $C$ 和 $D$ ，它满足三个方程： $n = 2$ 且 $m = 3$ 。这三个方程在 $t = 0, 1, 2$ 时要适配给定的值 $b = 6, 0, 0$ ： $t=0\kern 10pt如果\,\colorbox{lightgray}{$C+D\cdot 0=6$}，则第一个点在直线\,b=C+Dt\,上\\t=1\kern 10pt如果\,\colorbox{lightgray}{$C+D\cdot1=0$}，则第二个点在直线\,b=C+Dt上\\t=2\kern 10pt如果\,\colorbox{lightgray}{$C+D\cdot 2=0$}，则第三个点在直线\,b=C+Dt\,上$ 这个 $3\times2$ 的系统是没有解的： $\boldsymbol b=(6,0,0)$ 不是列 $(1, 1, 1)$ 和 $(0, 1, 2)$ 的线性组合。从这些方程中得到 $A，\boldsymbol x$ 和 $\boldsymbol b$ ： $A=\begin{bmatrix}1&0\\1&1\\1&2\end{bmatrix}\kern 8pt\boldsymbol x=\begin{bmatrix}C\\D\end{bmatrix}\kern 8pt\boldsymbol b=\begin{bmatrix}6\\0\\0\end{bmatrix}\kern 10ptA\boldsymbol x=\boldsymbol b\,无解$ 我们计算出 $\boldsymbol {\hat x}=(5,-3)$ ，这两个数字是最好的 $C$ 和 $D$ ，所以 $5 - 3 t$ 是对于这三个点最优的直线。我们要让投影和最小二乘联系起来，就要解释为什么 $A^TA\boldsymbol {\hat x}=A^T\boldsymbol b$ 。
在实际问题中，可能很轻易就有 $m = 100$ 个点而不是 $m = 3$ ，它们不会精确的匹配到任何直线 $C + D t$ 。本例中的数字 $6, 0, 0$ 夸大了误差，误差 $\boldsymbol e_1，\boldsymbol e_2$ 和 $\boldsymbol e_3$ 如 Figure 4.6所示。

二、最小化误差

我们如何让误差 $\boldsymbol e=\boldsymbol b-A\boldsymbol x$ 尽可能小呢？这个问题很重要，答案也很美妙。最好的 $\boldsymbol x$ （称为 $\boldsymbol{\hat x}$ ）可以通过几何学（误差 $\boldsymbol e$ 与 $A$ 列空间的夹角是 $90°$ ）；也可以由代数： $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ；还可以通过微积分得到相同的 $\boldsymbol {\hat x}$ ：误差 $||A\boldsymbol x-\boldsymbol b||^2$ 的导数在 $\boldsymbol{\hat x}$ 处为零。
由几何学： 每个 $A\boldsymbol x$ 都落在由列 $(1, 1, 1)$ 和 $(0, 1, 2)$ 所张成的平面内。在这个平面中，我们要找到最靠近 $\boldsymbol b$ 的点，这个最近的点就是投影 $\boldsymbol p$ 。
$A\boldsymbol {\hat x}$ 最好的选择就是 $\boldsymbol p$ ，最小的误差是 $\boldsymbol e=\boldsymbol b-\boldsymbol p$ ，它垂直于 $A$ 的列。这三个点的高度为 $p_1,p_2,p_3)$ 时，它们是在一条直线上，因为 $\boldsymbol p$ 在 $A$ 的列空间。拟合成一条直线时， $\boldsymbol{\hat x}$ 最好的选择就是 $(C, D)$ 。
由代数： 每个向量 $\boldsymbol b$ 都可以分为两部分，在列空间中的部分是 $\boldsymbol p$ ，还有垂直的部分是 $\boldsymbol e$ 。我们无法求解方程 $A\boldsymbol x=\boldsymbol b$ ，但是我们可以求解方程 $A\boldsymbol{\hat x}=\boldsymbol p$ （移除 $\boldsymbol e$ 然后求解 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ）： $A\boldsymbol x=\boldsymbol b=\boldsymbol p+\boldsymbol e\,无解\kern 10ptA\boldsymbol{\hat x}=\boldsymbol p\,有解\kern 10pt\boldsymbol{\hat x}\,就是\,(A^TA)^{-1}A^T\boldsymbol b\kern 15pt(4.3.1)$ $A\boldsymbol{\hat x}=\boldsymbol p$ 的解得到最小可能误差 $\boldsymbol e$ ： $\pmb{任意\,\boldsymbol x\,的长度平方}\kern 12pt||A\boldsymbol x-\boldsymbol b||^2=||A\boldsymbol x-\boldsymbol p||^2+||\boldsymbol e||^2\kern 14pt(4.3.2)$ 这就是直角三角形中的勾股定理 $c^2=a^2+b^2$ ，向量 $A\boldsymbol x-\boldsymbol p$ 在列空间中，它与在左零空间中的 $\boldsymbol e$ 垂直。我们通过选择 $\boldsymbol x=\boldsymbol{\hat x}$ 将 $A\boldsymbol x-\boldsymbol p$ 减小至零，留下了我们不能再减小的最小可能误差 $\boldsymbol e=(e_1,e_2,e_3)$ 。
注意 “最小” 的意思是 $A\boldsymbol x-\boldsymbol b$ 长度的平方最小化： $\pmb{最小二乘解\,\boldsymbol{\hat x}\,使得\,E=||A\boldsymbol x-\boldsymbol b||^2\,尽可能的小。}$ Figure 4.6a 展示了最近的直线，它的误差距离是 $e_1,e_2,e=1,-2,1$ 。这些都是竖直的距离，最小二乘直线最小化 $E=e_1^2+e_2^2+e_3^2$ 。
Figure 4.6b 显示了在三维空间（ $\boldsymbol b,\boldsymbol p,\boldsymbol e$ 的空间）中的同样问题。向量 $\boldsymbol b$ 不在 $A$ 的列空间中，这也就是为什么 $A\boldsymbol x=\boldsymbol b$ 的原因，没有任何直线可以同时穿过这三点。最小可能的误差就是垂直向量 $\boldsymbol e=\boldsymbol b-A\boldsymbol{\hat x}$ ，误差向量 $(1, - 2, 1)$ 在这三个方程中，它就是与这条最优直线的垂直距离。这两张图的背后就是基本方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。在这里插入图片描述
注意误差 $1, - 2, 1$ 相加等于零。这是因为：误差 $\boldsymbol e=(e_1,e_2,e_3)$ 垂直于 $A$ 的第一列 $(1, 1, 1)$ ，它们的点积就是 $e_1+e_2+e_3=0$ 。
由微积分： 大部分函数的最小值都是用微积分解决！图形降到最低点并且在每一个方向的导数都为零。这里要使得误差函数 $E$ 最小化， $E$ 就是平方和 $e_1^2+e_2^2+e_3^2$ （每个方程中误差的平方）： $E=||A\boldsymbol x-\boldsymbol b||^2=(C+D\cdot0-6)^2+(C+D\cdot1)^2+(C+D\cdot2)^2\kern 18pt(4.3.3)$ 未知数是 $C$ 和 $D$ ，两个未知数有两个偏导数，它们在极小值处都为零。称为偏导数是因为求 $\displaystyle\frac{\partial E}{\partial C}$ 时将 $D$ 当做常数，求 $\displaystyle\frac{\partial E}{\partial D}$ 时将 $C$ 当做常数： $\frac{\partial E}{\partial C}=2(C+D\cdot0-6)+2(C+D\cdot1)+2(C+D\cdot2)=0\\\frac{\partial E}{\partial D}=2(C+D\cdot0-6)\cdot(0)+2(C+D\cdot1)\cdot(1)+2(C+D\cdot2)\cdot(2)=0$ 因为求导的链式法则， $\displaystyle\frac{\partial E}{\partial D}$ 含有额外的因子 $0, 1, 2$ （最后的 $C+2D)^2$ 的导数是 $2$ 乘 $C + 2 D$ 再乘 $2$ ）。这些因子在 $\displaystyle\frac{\partial E}{\partial C}$ 中就是 $1, 1, 1$ 。
毫无意外的 $||A\boldsymbol x-\boldsymbol b||$ 导数的因子 $1, 1, 1$ 和 $0, 1, 2$ 就是 $A$ 的列，现在消去每一项的 $2$ 然后将 $C$ 和 $D$ 的同类项进行合并： $\begin{matrix}C\,的导数为零：3C+3D=6\\D\,的导数为零：3C+5D=0\end{matrix}\kern 15pt\pmb{矩阵}\begin{bmatrix}3&3\\3&5\end{bmatrix}就是\,A^TA\kern 16pt(4.3.4)$ 这些方程与 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 完全相同。 最优的 $C$ 和 $D$ 就是 $\boldsymbol{\hat x}$ 的分量。从微积分得到的方程与从线性代数得到的 “正规方程” 完全一致。这是最小二乘的关键方程：

$当\,A^TA\boldsymbol{\hat x}=A^T\boldsymbol b\,时，||A\boldsymbol x-\boldsymbol b||^2\,的偏导数为零。$

解是 $C = 5 ， D = - 3$ ，因此 $b = 5 - 3 t$ 就是最优的直线 —— 它与这三点最近。当 $t = 0, 1, 2$ 时，这条直线经过 $\boldsymbol p=5,2,-1$ ，它不通过 $\boldsymbol b=6,0,0$ 。误差是 $1, - 2, 1$ ，就是误差向量 $\boldsymbol e$ ！

三、最小二乘的大图

Figure 4.3 是关键的一张图，它显示了 $4$ 个基本子空间和矩阵的真正作用，左边的向量 $\boldsymbol x$ 到右边的 $\boldsymbol b=A\boldsymbol x$ ，这张图将 $\boldsymbol x$ 分成 $\boldsymbol x_r+\boldsymbol x_n$ ， $A\boldsymbol x=\boldsymbol b$ 可能存在很多解。
现在的情况与上述正好相反， $A\boldsymbol x=\boldsymbol b$ 无解，这里我们不再分解 $\boldsymbol x$ ，而是将 $\boldsymbol b$ 分成 $\boldsymbol b=\boldsymbol p+\boldsymbol e$ 。Figure 4.7 显示了最小二乘的大图，我们不解 $A\boldsymbol x=\boldsymbol b$ ,而是求解 $A\boldsymbol{\hat x}=\boldsymbol p$ ，误差 $\boldsymbol e=\boldsymbol b-\boldsymbol p$ 无法避免。在这里插入图片描述
注意这里的零空间 $\pmb N(A)$ 很小 —— 仅有一个点，这是因为 $A$ 的列线性无关， $A\boldsymbol x=\boldsymbol 0$ 的唯一解就是 $\boldsymbol x=\boldsymbol 0$ ，此时 $A^TA$ 是可逆的。方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 就决定了最优向量 $\boldsymbol{\hat x}$ ，对于误差有 $A^T\boldsymbol e=\boldsymbol 0$ 。

四、拟合一条直线

最小二乘最常见的应用就是拟合一条直线，从 $m > 2$ 个点开始，希望找到一条最接近这些点的直线。在时间点 $t_1,t_2,\cdots,t_m$ 时，这 $m$ 个点的高度是 $b_1,b_2,\cdots,b_m$ 。最优直线 $C + D t$ 距这些点的垂直距离是 $e_1,e_2,\cdots,e_m$ 。没有完全经过这些点的直线，我们要使得这些误差的平方 $E=e_1^2+e_2^2+\cdots+e_m^2$ 最小。
第一个例子是 Figure 4.6 中的三点问题，现在我们将允许 $m$ 个点（ $m$ 可以很多）， $\boldsymbol{\hat x}$ 的两个分量仍然是 $C$ 和 $D$ 。
如果我们要求 $A\boldsymbol x=\boldsymbol b$ 的确切解时，它需要一条完全通过这 $m$ 个点的直线，通常这是不现实的。两个未知数 $C$ 和 $D$ 决定了一条直线，所以 $A$ 仅有 $n = 2$ 列。为了拟合这 $m$ 个点，我们需要求解 $m$ 个方程（但是只有两个未知数！） $A\boldsymbol x=\boldsymbol b\kern 5pt是\kern 5pt\begin{matrix}C+Dt_1=b_1\\C+Dt_2=b_2\\\vdots\\C+Dt_m=b_m\end{matrix}\kern 10pt其中\,A=\begin{bmatrix}1&t_1\\1&t_2\\\vdots&\vdots\\1&t_m\end{bmatrix}\kern 15pt(4.3.5)$ $A$ 的列空间很细，这使得大部分的 $\boldsymbol b$ 都落在列空间外。如果 $\boldsymbol b$ 恰好落在列空间中，那么这些点就很正好落在一条直线上。这种情况下 $\boldsymbol b=\boldsymbol p$ ， $A\boldsymbol x=\boldsymbol b$ 可解，误差是 $\boldsymbol e=(0,0,\cdots,0)$ 。 ${\color{blue}\begin{array}{l}最接近的直线\,C+Dt\,的高度是\,p_1,p_2,\cdots,p_m，误差是\,e_1,e_2,\cdots,e_m.\\求解\,A^TA\boldsymbol{\hat x}=A^T\boldsymbol b，得到\,\boldsymbol{\hat x}=(C,D)，误差是\,e_i=b_i-C-Dt_i.\end{array}}$ 直线拟合非常重要，我们给定方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ， $A$ 的两列是线性无关的（除非所有的点 $t_i$ 都一样），我们转到最小二乘并求解 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。 $\pmb{点积矩阵}\kern 6ptA^TA=\begin{bmatrix}1&1&\cdots&1\\t_1&t_2&\cdots&t_m\end{bmatrix}\begin{bmatrix}1&t_1\\1&t_2\\\vdots&\vdots\\1&t_m\end{bmatrix}=\begin{bmatrix}m&\sum t_i\\\sum t_i&\sum t^2_i\end{bmatrix}\kern 13pt(4.3.6)$ 正规方程的右侧是 $2\times1$ 的向量 $A^T\boldsymbol b$ ： $A^T\boldsymbol b=\begin{bmatrix}1&1&\cdots&1\\t_1&t_2&\cdots&t_m\end{bmatrix}\begin{bmatrix}b_1\\b_2\\\vdots\\b_m\end{bmatrix}=\begin{bmatrix}\sum b_i\kern 8pt\\\sum t_ib_i\end{bmatrix}\kern 15pt(4.3.7)$ 在特定问题中，这些数字都是给定的，最优的 $\boldsymbol{\hat x}=(C,D)$ 是 $(A^TA)^{-1}A^T\boldsymbol b$ .

当 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 时，直线 $C + D t$ 使得 $e_1^2+e_2^2+\cdots+e_m^2=||A\boldsymbol x-\boldsymbol b||^2$ 最小： $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b\kern 20pt\begin{bmatrix}m&\sum t_i\\\sum t_i&\sum t^2_i\end{bmatrix}\begin{bmatrix}C\\D\end{bmatrix}=\begin{bmatrix}\sum b_i\\\sum t_ib_i\end{bmatrix}\kern 18pt(4.3.8)$

这 $m$ 个点与这条直线的垂直误差是 $\boldsymbol e=\boldsymbol b-\boldsymbol p$ 的分量，这个误差向量 $\boldsymbol b-A\boldsymbol{\hat x}$ 与 $A$ 的列是垂直的（几何学），误差向量在 $A^T$ 的零空间中（线性代数）。最优的 $\boldsymbol{\hat x}=(C,D)$ 使得总误差 $E$ 最小，即平方和最小（微积分）： $E(\boldsymbol x)=||A\boldsymbol x-\boldsymbol b||^2=(C+Dt_1-b_1)^2+(C+Dt_2-b_2)+\cdots(C+Dt_m-b_m)^2$ 微积分中令 $\displaystyle\frac{\partial E}{\partial C}$ 和 $\displaystyle\frac{\partial E}{\partial D}$ 等于零，就可以得到 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。
其它的最小二乘问题未知数可能多于两个，拟合成最优的抛物线有 $n = 3$ 个系数 $C, D, E$ 。一般情况下，我们用 $n$ 个参数拟合 $m$ 个数据点，矩阵 $A$ 有 $n$ 列，且 $n < m$ 。 $||A\boldsymbol x-\boldsymbol b||^2$ 的导数可以得到 $n$ 个方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。平方的导数是线性的。这就是为什么最小二乘法会流行。

【例2】当测量的时间点 $t_i$ 的和为零时， $A$ 有正交列。假设在时间点 $t = - 2, 0, 2$ 时 $b = 1, 2, 4$ ，这些时间点相加为零。 $A$ 的列点积为零： $(1, 1, 1)$ 与 $(- 2, 0, 2)$ 正交。 $\begin{array}{l}C+D(-2)=1\\C+D(0)=2\\C+D(2)=4\end{array}或\kern 5ptA\boldsymbol x=\begin{bmatrix}1&-2\\1&0\\1&2\end{bmatrix}\begin{bmatrix}C\\D\end{bmatrix}=\begin{bmatrix}1\\2\\4\end{bmatrix}$ 当 $A$ 的列正交时， $A^TA$ 是一个对角矩阵： $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b\kern 5pt是\kern 5pt\begin{bmatrix}3&\pmb0\\\pmb0&8\end{bmatrix}\begin{bmatrix}C\\D\end{bmatrix}=\begin{bmatrix}7\\6\end{bmatrix}\kern 15pt(4.3.9)$ 重点：因为 $A^TA$ 是对角矩阵，我们可以分开求出 $C=\displaystyle\frac{7}{3}$ ， $D=\displaystyle\frac{6}{8}$ 。 $A^TA$ 中的零是 $A$ 垂直列的点积。这个对角矩阵 $m = 3$ ， $t_1^2+t_2^2+t_3^2=8$ ，它与单位矩阵基本上有差不多的性质。
正交矩阵非常有用，我们可以通过减去时间的平均 $\bar t=(t_1+t_2+\cdots+t_m)/m$ 来进行时间的移位。如果原始时间点是 $1, 3, 5$ ，则它们的平均值是 $\bar t=3$ ，平移后的时间点 $T=t-\bar t=t-3$ ，它们的和是零！ $\begin{array}{l}T_1=1-3=-2\\T_2=3-3=0\\T_3=5-3=2\end{array}\kern 12ptA_{new}=\begin{bmatrix}1&T_1\\1&T_2\\1&T_3\end{bmatrix}\kern 12ptA^T_{new}A_{new}=\begin{bmatrix}3&\pmb0\\\pmb0&8\end{bmatrix}$ 此时 $C$ 和 $D$ 可以由简单的方程（4.3.9）直接得到，最优的直线 $C + D T$ 就是 $C+D(t-\bar t)=C+D(t-3)$ 。

五、A 有相关列： $\boldsymbol{\hat x}$ 是什么？

从一开始，我们就假设 $A$ 的列是无关的，那么 $A^TA$ 可逆， $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 可以得到 $A\boldsymbol x=\boldsymbol b$ 的最小二乘解。
如果 $A$ 有相关列，最优的 $\boldsymbol{\hat x}$ 是什么呢？下面是一个特殊的例子：在这里插入图片描述测量值 $b_1=3$ 和 $b_2=1$ 都是在相同的时间点 $T$ ！直线 $C + D t$ 不可能同时经过这两个点。此时我们将 $\boldsymbol b=(3,1)$ 投影到 $A$ 的列空间得到 $\boldsymbol p=(2,2)$ ，这样就把方程 $A\boldsymbol x=\boldsymbol b$ 变成了 $A\boldsymbol{\hat x}=\boldsymbol p$ ，将一个无解的方程变成了一个又无穷多解的方程，这个问题的原因就是 $A$ 有相关列且 $(1, - 1)$ 在它的零空间中。
我们应该选择哪个解 $\boldsymbol{\hat x}$ 呢？图中所有的虚线在 $T = 1$ 处都有两个误差 $1$ 和 $- 1$ ，误差 $\boldsymbol e=\boldsymbol b-\boldsymbol p=(1,-1)$ 需要尽可能小，但是这里我们无法知道哪条是最优的直线。
在直觉上会选择高度为 $2$ 的水平线，如果最优直线的方程是 $b = C + D t$ ，则可以选择 $\hat x_1=C=2$ ， $\hat x_2=D=0$ ，如果最优方程是 $b = c t + d$ （交换了一下 $C$ 和 $D$ 的位置），则水平线有 $\hat x_1=c=0$ ， $\hat x_2=d=2$ ，即不同的 $\boldsymbol{\hat x}$ 得到了相同的直线！这样同样无法选择最优的 $\boldsymbol{\hat x}.$
$A$ 的 “伪逆矩阵（pseudoinverse）” 会选出 $A\boldsymbol{\hat x}=\boldsymbol p$ 的最短解。这里的最短解是 $\boldsymbol x^+=(1,1)$ ，它是 $A$ 行空间中的特解， $\boldsymbol x^+$ 的长度是 $\sqrt2$ （上述两个 $\boldsymbol{\hat x}=(2,0)$ 和 $(0, 2)$ 的长度都为 $2$ ）。
当 $A$ 的列无关时，零空间仅包含零向量，伪逆矩阵就是我们正常的左逆矩阵 $L=(A^TA)^{-1}A^T$ 。
注意：MATLAB 的奇异矩阵得到的是 $\textrm{Inf}$ 或 $\textrm{NaN}$ （Not a Number）或 $10^{16}$ （不好的数），每种情况都会出现警告！ $\textrm{Inf}$ 和 $\textrm{NaN}$ 和 $10^{16}$ 可能是来自 $\boldsymbol {0x}=\boldsymbol b$ 和 $\boldsymbol{0x}=\boldsymbol 0$ 和 $\boldsymbol{10^{-16}x}={\boldsymbol 1}$ 。
这三个例子是三大困难：奇异且无解，奇异有无穷多解，非常非常接近奇异。

六、抛物线拟合

如果我们扔出去一颗球，要用一条直线去拟合球的轨迹是疯狂的行为。一条抛物线 $b=C+Dt+Et^2$ 将允许这个球先向上再向下（ $b$ 是在时间 $t$ 时的高度）。实际轨迹并不是一条完美的抛物线，但是投影的全部定理都是由近似开始的。
当伽利略从比萨斜塔扔下一块石头开始，它会加速，距离包含一个二次项 $\displaystyle\frac{1}{2}gt^2$ （伽利略的观点与石头的重量无关）。如果没有 $t^2$ 我们将无法将卫星送入轨道。尽管这里出现了一个非线性函数 $t^2$ ，未知数 $C, D, E$ 仍然是线性的！最优抛物线的拟合仍然是线性代数的一个问题。
问题： 在时间 $t_1,t_2,\cdots,t_m$ 时，高度是 $b_1,b_2,\cdots,b_m$ ，用一条抛物线拟合这些点。
解：当 $m > 3$ 个点时，这 $m$ 个方程一般是无解的： $\begin{matrix}C+Dt_1+Et_1^2=b_1\\C+Dt_2+Et_2^2=b_2\\\cdots\\C+Dt_m+Et_m^2=b_m\end{matrix}\kern 5pt\begin{matrix}是\,A\boldsymbol x=\boldsymbol b，其中\,\\A\,是\,m\times3\,的矩阵\end{matrix}\kern 5ptA=\begin{bmatrix}1&t_1&t_1^2\\1&t_2&t_2^2\\\vdots&\vdots&\vdots\\1&t_m&t_m^2\end{bmatrix}\kern 15pt(4.3.10)$ 最小二乘： 最接近的抛物线 $C+Dt+Et^2$ 是 $\boldsymbol{\hat x}=(C,D,E)$ 满足三个正规方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。
下面将这个变成投影问题： $A$ 的列空间维度是 $3$ ， $\boldsymbol b$ 的投影是 $\boldsymbol p=A\boldsymbol{\hat x}$ ，其中 $\boldsymbol p$ 是系数 $C, D, E$ 对三列进行组合。第一个数据点的误差是 $e_1=b_1-C-Dt_1-Et_1^2$ ，误差的总平方和是 $e_1^2+e_2^2+\cdots+e_m^2$ 。如果用微积分来最小化误差，就是计算 $E$ （误差总的平方和，与系数 $E$ 不同）对于 $C, D, E$ 的偏导数，当 $\boldsymbol{\hat x}=(C,D,E)$ 是 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 这个 $3\times3$ 的方程系统的解时，这三个偏导数等于零。
傅里叶级数也是最小二乘的应用 —— 这里是近似函数而不是近似向量，函数的近似从误差的平方和 $e_1^2+e_2^2+\cdots+e_m^2$ 变成了误差平方的积分。

【例3】一个抛物线 $b=C+Dt+Et^2$ 当 $t = 0, 1, 2$ 时，经过三个高度 $b = 6, 0, 0$ ，关于 $C, D, E$ 的三个方程是 $\begin{matrix}C+D\cdot0+E\cdot0^2=6\\C+D\cdot1+E\cdot1^2=0\\C+D\cdot2+E\cdot2^2=0\end{matrix}\kern 35pt(4.3.11)$ 这个就是 $A\boldsymbol x=\boldsymbol b$ ，它是有解的，这三个数据点得到三个方程和一个方阵，解是 $\boldsymbol x=(C,D,E)=(6,-9,3)$ ，如 Figure 4.8a 所示，这条抛物线 $b=6-9t+3t^2$ 通过这三个点。
上述对于投影来说意味着什么？矩阵有三列，它张成整个 $R3 \pmb{\textrm R}^3$ 空间，投影矩阵就是单位矩阵， $\boldsymbol b$ 的投影就是 $\boldsymbol b$ ，误差是零。这里不需要 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ，因为 $A\boldsymbol x=\boldsymbol b$ 有解。当然两边也可以左乘 $A^T$ ，只是我们没什么理由这样做。
Figure 4.8 显示了在 $t_4$ 时刻的第四点是 $b_4$ ，如果它正好落在抛物线上，则新的 $A\boldsymbol x=\boldsymbol b$ （ $4$ 个方程）任然有解，如果第四点不在抛物线上，我们就需要 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ，那么最小二乘的抛物线就不一定让四个点都落在它上面。
最小二乘会平衡四个误差得到三个关于 $C, D, E$ 的方程。在这里插入图片描述

七、主要内容总结

最小二乘解 $\boldsymbol{\hat x}$ 使得 $||A\boldsymbol x-\boldsymbol b||^2=\boldsymbol x^TA^TA\boldsymbol x-2\boldsymbol x^T\boldsymbol A^T\boldsymbol b+\boldsymbol b^T\boldsymbol b$ 最小，这个就是 $E$ ，它是 $m$ 个方程误差的平方和（ $m > n$ ）。
最优的 $\boldsymbol{\hat x}$ 由正规方程 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 得来， $E$ 是最小值。
通过直线 $b = C + D t$ 拟合 $m$ 个点，正规方程可以得到 $C$ 和 $D$ 。
最优直线的高度是 $\boldsymbol p=(p_1,p_2,\cdots,p_m)$ ，和数据点的竖直距离就是误差 $\boldsymbol e=(e_1,e_2,\cdots,e_m)$ ，关键方程是 $A^T\boldsymbol e=\boldsymbol 0$ 。
如果有 $m$ 个数据点（ $m > n$ ），我们可以得到 $n$ 个方程，这 $n$ 个方程 $A\boldsymbol x=\boldsymbol b$ 通常是无解的。但是 $n$ 个 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 方程可以得到最小二乘解 —— 最小 MSE（均方值：mean square error）的组合。

八、例题

【例4】在时刻 $t=1,2,\cdots,9$ 时， $9$ 个测量值 $b_1$ 到 $b_9$ 全为零。第 $10$ 个测量值 $b_{10}=40$ ，它是一个离群值。找到一条最优的水平直线 $y = C$ 来拟合这 $10$ 个点 $(1,0),(2,0),\cdots,(9,0),(10,40)$ ，使用下面三种情况的误差 $E$ ：
（1）最小二乘误差 $E_2=e_1^2+e_2^2+\cdots+e^2_{10}$ （ $C$ 的正规方式是线性的）
（2）最小极大误差（least maximum error） $E_\infty=|e_{\textrm{max}}|$
（3）最小误差和 $E_1=|e_1|+|e_2|+\cdots+|e_{10}|$
解：（1）最小二乘法对 $0,0,\cdots,0,40$ 这 $10$ 个点拟合出来的水平直线是 $C = 4$ ： $A\,是\,10\times1\,的全\,1\,矩阵\kern 10ptA^TA=10\kern 10ptA^T\boldsymbol b=\sum^{10}_{i=1}b_i=40\kern 10pt所以\,10C=40$ （2）最小极大误差是 $C = 20$ ，是 $0$ 和 $40$ 的一半。
（3）最小误差和是 $C = 0$ 。误差和是 $9∣ C ∣ + ∣40 - C ∣$ ，可得当 $C = 0$ 时最小。
最小和来自中位测量（median measurement）， $0,0,\cdots,0,40$ 的中位数是 $0$ 。统计学中，最小二乘解受离群值影响很大，例如此例中的 $b_{10}=40$ ，有时使用最小误差和会比较好。但是这样的话，方程就是一个非线性的了。
下面使用最小二乘法用一条直线 $C + D t$ 来拟合这 $10$ 个点， $(1, 0)$ 到 $(10, 40)$ ： $A^TA=\begin{bmatrix}10&\sum t_i\\\sum t_i&\sum t^2_i\end{bmatrix}=\begin{bmatrix}10&55\\55&385\end{bmatrix}\kern 20ptA^T\boldsymbol b=\begin{bmatrix}\sum b_i\\\sum t_ib_i\end{bmatrix}=\begin{bmatrix}40\\400\end{bmatrix}$ 这些就是由方程（4.3.8）而来，然后求解 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ ，得到 $C = - 8 ， D = 24/11$ 。
如果将 $\boldsymbol b=(0,0,\cdots,40)$ 乘上 $3$ 再加上 $30$ 得到 $\boldsymbol b_{\textrm{new}}=(30,30,\cdots,150)$ ，那么 $C$ 和 $D$ 会发生什么样的变化呢？线性将允许我们重新设定 $\boldsymbol b$ 的比例， $3$ 乘 $\boldsymbol b$ 可以得到 $C$ 和 $D$ 都乘 $3$ ，对所有的 $b_i$ 都加上 $30$ ，相当于对 $C$ 加上 $30$ ，得到新的 $C_{\textrm {new}}=3C+30=6$ ， $D_{\textrm{new}}=3D=72/11$ 。

【例5】在时刻 $t = - 2, - 1, 0, 1, 2$ 时， $\boldsymbol b=(0,0,1,0,0)$ ，找到一条离这些点最近（最小二乘误差）的抛物线 $C+Dt+Et^2$ 。首先写出 $5$ 个方程 $A\boldsymbol x=\boldsymbol b$ ，它有三个未知数 $\boldsymbol x=(C,D,E)$ ，这些方程分别过这 $5$ 个点。但是它们没有解，因为并不存在这样的抛物线，我们要求解 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 。
这里可以预测 $D = 0$ ，最优的抛物线应该是关于 $t = 0$ 对称。在 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 中，方程 $2$ 中的 $D$ 与方程 $1$ 和 $3$ 是解耦的。
解： $A\boldsymbol x=\boldsymbol b$ 的 $5$ 个方程有一个矩形的范德蒙德（Vandermonde）矩阵 $A$ ： $\begin{matrix}C+D(-2)+E(-2)^2=0\\C+D(-1)+E(-1)^2=0\\C+D\kern 8pt(0)+E\kern 8pt(0)^2=1\\C+D\kern 8pt(1)+E\kern 8pt(1)^2=0\\C+D\kern 8pt(2)+E\kern 8pt(2)^2=0\end{matrix}\kern 10ptA=\begin{bmatrix}1&-2&4\\1&-1&1\\1&\kern 7pt0&0\\1&\kern 7pt1&1\\1&\kern 7pt2&4\end{bmatrix}\kern 10ptA^TA=\begin{bmatrix}5&\pmb0&10\\\pmb0&10&\pmb0\\10&\pmb0&34\end{bmatrix}$ $A^TA$ 中的零表示 $A$ 的列 $2$ 和它的列 $1$ 和列 $3$ 正交，在 $A$ 中可以直接看出来（时间 $- 2, - 1, 0, 1, 2$ 是对称的）。抛物线 $C+Dt+Et^2$ 中最优的 $C, D, E$ 由 $A^TA\boldsymbol{\hat x}=A^T\boldsymbol b$ 解得，这里 $D$ 和 $C$ 与 $E$ 是解耦的： $\begin{bmatrix}5&0&10\\0&10&0\\10&0&34\end{bmatrix}\begin{bmatrix}C\\D\\E\end{bmatrix}=\begin{bmatrix}1\\0\\0\end{bmatrix}\kern 5pt解得\kern 5pt\begin{array}{l}C=\displaystyle\frac{17}{35}\\D=0，与预测一致\\E=-\displaystyle\frac{1}{7}\end{array}$