最小二乘问题详解1：线性最小二乘

1. 引言

最小二乘可以说是现代科学与工程的“隐形骨架”，几乎无处不在。比如：

1. 测绘与空间信息科学：摄影测量平差、GNSS/RTK定位、控制网平差。
2. 机器人与自动驾驶：视觉SLAM/LiDAR SLAM、传感器融合、手眼标定、运动学/动力学参数辨识。
3. 计算机视觉与计算机图形学：图像配准、三维重建、光照估计、纹理映射。
4. 机器学习与数据科学：线性回归、岭回归、主成分分析、支持向量机、神经网络训练。

笔者之前对最小二乘问题也只是一知半解，这里就详细学习总结一下。

2. 最小二乘

2.1 定义

最小二乘是一种从有误差的数据中寻找最佳拟合模型的数学方法，它的核心思想是让模型的预测值与实际观测值之间的“误差平方和”最小。

比如经典的最小二乘拟合直线的问题：给定一组有噪声的数据点，需要拟合一条直线$y=kx+b$，那么不可能所有点都正好在一条直线上，合理的方案是找到最佳的斜率$k$和截距$b$，使得所有点到这条直线的竖直距离的平方和最小。

最小二乘的数学表达为：

$$\underset{\theta }{\min }\sum _{i=1}^m{r}_i(\theta {)}^2=\underset{\theta }{\min }\Vert \mathbf{r}(\theta ){\Vert }^2$$

其中：

• $r_i(\theta )$ 是第 $i$ 个观测的残差（residual）：
  $r_i=y_i-f(x_i;\theta )$
• $\mathbf{r}(\theta )$ 是所有残差组成的向量
• $\theta$ 是待估计的参数向量

虽然定义出来了，但是另一个问题是——为什么最小二乘用“平方和”而不是“绝对值和”、“四次方和”或其他方式？这背后其实有深刻的数学原理：

• 从统计学的角度上来讲，最小二乘就是在误差服从高斯分布时的最大似然估计。
• 从几何的角度上来讲，平方和是欧氏距离的平方，是最自然的距离度量。

不过要说清楚这两点有点麻烦，我们可以先暂时通过高数知识来简单的理解。函数$f(r)=r^2$是一个凸函数，所谓凸函数，直观来说就是任意两点之间的线段始终在函数图像之上，只有一个“谷底”，这个“谷底”就是全局最小值。这意味着任何局部最小值就是全局最小值，在求解优化问题的时候，可以通过梯度下降等算法收敛到全局最优。

2.2 线性

最小二乘问题可以分为线性最小二乘和非线性最小二乘来讨论。首先，我们先来讨论一个比较本质的问题，什么叫做线性？在《初等线性代数》中，线性指的是可加性和齐次性，例如一个变换$T$能满足如下两个条件：

1. $T(x+y)=T(x)+T(y)$
2. $T(\alpha x)=\alpha T(x)$

突然地引入数学上的定义确实有点难以理解，不过我们只需要明白，线性是一种非常优良的性质。比如说，满足线性的函数/变换显然是连续的、可导的以及光滑的，这意味着这个函数/变换不仅结构简单，也易于预测和控制。科学家和工程师都喜欢假设问题的模型是线性的开始研究，即使真实世界的问题模型大多数是非线性的，也会通过数学方法将非线性问题转换成线性问题。因此，要研究最小二乘，首先需要理解线性最小二乘。

3. 线性最小二乘

3.1 定义

需要明确指出的是，问题模型的线性还是非线性，是相对于待定参数$\theta$而言的，而不是已知参数$x$。线性最小二乘的问题模型可以写成如下形式：

$$f(x;\theta )=A\theta$$

那么，线性最小二乘的数学表达为：

$$\underset{\theta }{\min }\Vert A\theta -b{\Vert }^2\text{\hspace{1em}(2)}$$

其中：

• $A$：设计矩阵（$m\times n$，$m$ 是数据点数，$n$ 是参数数）
• $\theta$：未知参数向量（$n\times 1$）
• $b$：观测向量（$m\times 1$）

3.2 具体化

数学上的概念比较抽象，这里还是结合前面最小二乘拟合直线的例子来理解。给定一组有噪声的数据点：

$$(x_1,y_1),(x_2,y_2),\dots ,(x_m,y_m)$$

我们希望拟合一条直线：

$$y=kx+b_0$$

其中 $k$ 是斜率，$b_0$ 是截距。很显然，对于待定参数$k$和$b_0$来说，这个问题模型是线性的，需要使用线性最小二乘来估计参数。

将数据点带入这个问题模型，可得方程组：

$$\{\begin{array}{l}{y}_1=k{x}_1+b_0\\ y_2=k{x}_2+b_0\\ ⋮\\ y_m=k{x}_m+b_0\end{array}$$

将方程组写成矩阵形式：

$$\left[\begin{array}{cc}{x}_1& 1\\ x_2& 1\\ ⋮& ⋮\\ x_m& 1\end{array}\right]\left[\begin{array}{c}k\\ b_0\end{array}\right]=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]\text{\hspace{1em}(1)}$$

令：

• 设计矩阵：$A=\left[\begin{array}{cc}{x}_1& 1\\ x_2& 1\\ ⋮& ⋮\\ x_m& 1\end{array}\right]$（$m\times 2$）
• 参数向量：$\theta =\left[\begin{array}{c}k\\ b_0\end{array}\right]$（$2\times 1$）
• 观测向量：$b=\left[\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right]$（$m\times 1$）
• 问题模型函数：$f(x;\theta )=kx+b_1=\left[\begin{array}{cc}x& 1\end{array}\right]\left[\begin{array}{c}k\\ b_1\end{array}\right]$

那么问题模型的残差就是：

$$\mathbf{r}(\theta )=A\theta -b$$

线性最小二乘问题可归纳为：

$$\underset{\theta =[k,b_1{]}^T}{\min }\Vert A\theta -b{\Vert }^2$$

3.3 求解

先不谈如何求解最小二乘公式(2)的问题，先说说如何解决方程组(1)，毕竟如果能正确求解方程组(1)，那么这个问题就解决了。很显然，方程组(1)就是《初等线性代数》中的线性方程组，根据《初等线性代数》中的知识，这种方程个数$m$比未知数多的线性方程组$n$是没有解的。但是，归结到具体的显式问题中来说，这个方程组应该要有解：假设所有的数据点$(x_i,y_i)$都没有噪声，那么选取任意$n$组数据即可计算出唯一解。但是真实世界的数据是有噪声的，不能这么做。

回忆《初等线性代数》中的知识，求解线性方程组$A\theta =b$最容易理解就是矩阵求逆法，但是这个方程组$m$要远大于$n$，明显是没办法求解逆矩阵的。但是我们可以改造这个方程组，在两边都乘以相同的矩阵$A^T$:

$$A^{T}A\theta =A^{T}b$$

这个方程就是正规方程。$A^{T}A$是方阵，在满秩的情况下可以求逆矩阵，其解为：

$${\theta }^{\ast }=(A^{T}A{)}^{-1}{A}^{T}b\text{\hspace{1em}(3)}$$

这个解其实就是最小二乘公式(2)的解，即最小二乘解。

3.4 原理

为什么说上文的式(3)恰好就是式(2)的最小二乘解呢？为什么我们会知道在两边都乘以相同的矩阵$A^T$呢？这里就来推导一下。

3.4.1 代数推导

之前已经提到过，最小二乘是“误差平方和”，是一个凸函数，可以求它的极小值。令

$$J(\theta )=\Vert A\theta -b{\Vert }^2$$

根据《高等数学》中的知识，要求函数的极小值，需要对 $\theta$ 求导，并令梯度（导数）为 0：

$$\frac{\partial J}{\partial \theta }=2(A\theta -b)\frac{\partial A\theta -b}{\partial \theta }=0$$

根据矩阵微积分的知识，$f(\theta )=a^T\theta$的导数是$a$，因此：

$$2A^T(A\theta -b)=0$$

调换位置，也就得到了正规方程：

$$A^{T}A\theta =A^{T}b$$

3.4.2 几何推导

在回答这个问题之前，我们必须要对《线性代数》中的矩阵有更深刻的认识：矩阵的列向量张成了一个​列空间（Column Space）​​，由该矩阵所有列向量的线性组合所构成。而矩阵与向量相乘的结果，正是这些列向量以向量中对应分量为系数的线性组合。例如，设矩阵 $A=\left[\begin{array}{cccc}{\mathbf{a}}_1& {\mathbf{a}}_2& \cdots & {\mathbf{a}}_n\end{array}\right]$，其中 ${\mathbf{a}}_i$ 是列向量。对于任意向量 $\mathbf{x}=\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]$，有：

$$A\mathbf{x}=x_1{\mathbf{a}}_1+x_2{\mathbf{a}}_2+\cdots +x_n{\mathbf{a}}_n$$

这个结果 $A\mathbf{x}$ 显然是矩阵 $A$ 的列向量的一个线性组合，因此它属于列空间。所以，矩阵乘以一个向量的结果，是其列向量的一个线性组合，且这个结果落在矩阵的列空间中。

那么，对于线性最小二乘问题$A\theta =b$中来说，观测向量$b$会落到设计矩阵$A$的列空间中吗？由于噪声的存在，肯定是不行的，只能尽量寻找一个$\theta$，使得$A\theta$尽量靠近$b$。那么什么样的$\theta$才能满足尽可能接近的要求呢？答案很简单，就是做正交投影。形象的解释就是，一个向量$b$投影平面$A$的影子$A{\theta }^{\ast }$才是最接近$b$的，并且最接近的投影方式是正交投影，而这个$\theta$就是最小二乘解${\theta }^{\ast }$。

所谓正交投影，指的是一个点向一个平面（或直线）作垂线，垂足就是投影点；也就是说，$b-A\theta$应该垂直于$A$的列空间。这也意味着，$b-A\theta$与$A$的每一个列向量都正交，那么就有

$$A^T(b-A\theta )=0$$

调换位置，同样得到正规方程：

$$A^{T}A\theta =A^{T}b$$

以上推论也说明了一个原理：在欧几里得空间中，点到子空间的最短欧式距离，是通过正交投影实现的，最小二乘利用的就是这个原理。