线性方程组的最小二乘解

1. 背景

对于线性方程组可以根据方程数量与未知数数量分为欠定方程、适定方程和超定方程，欠定方程是指方程数量少于未知数数量，适定方程是指方程数量等于未知数数量，超定方程是指方程数量多于未知数数量。对于存在存在唯一解的情况下，只需要简单的矩阵运算就可以得到解，如：
$Ax=b\ -> \ x=A^{-1}b$ ,但实际工程中几乎只会出现超定方程组，而且是一个无解的超定方程组。虽然通常情况下超定方程组是无解的，但为其求解一个近似解，也称作最小二乘解。

2. 线性方程组的最小二乘解-代数推导

在工程中通常有一堆测量数据，然后为这个系统找到一个数学模型：

$f(x_1,x_2···x_n)=y\\ a_1x_1+a_2x_2+···+a_nx_n=y$

将这一堆测量数据带入模型：

$a_{1}x_{11}+a_{2}x_{12}+···+a_{n}x_{1n}=y_1\\ a_{1}x_{21}+a_{2}x_{22}+···+a_{n}x_{2n}=y_2\\ ···\\ a_{1}x_{m1}+a_{2}x_{m2}+···+a_{n}x_{mn}=y_m$

上述线性方程组中，x并不是未知量，相反a才是未知量，所以上述是关于a的线性方程组，根据最小二乘的思想---a取何值时误差和最小：

$e(a_1,a_2···a_n)=\sum_{n=1}^m (a_{1}x_{i1}+a_{2}x_{i2}+···+a_{n}x_{in}-y_i)^2$

分别对a_1~a_n求偏导：

$\frac{\partial e(a_1,a_2···a_n)}{\partial a_1}=2\sum_{i=1}^m x_{i_1}(a_{1}x_{i1}+a_{2}x_{i2}+···+a_{n}x_{in}-y_i)=0\\ \frac{\partial e(a_1,a_2···a_n)}{\partial a_2}=2\sum_{i=1}^m x_{i_2}(a_{1}x_{i1}+a_{2}x_{i2}+···+a_{n}x_{in}-y_i)=0\\···\\ \frac{\partial e(a_1,a_2···a_n)}{\partial a_n}=2\sum_{i=1}^m x_{i_n}(a_{1}x_{i1}+a_{2}x_{i2}+···+a_{n}x_{in}-y_i)=0$

上式进一步化简得到：

$a_1\sum_{i=1}^m x_{i1}x_{i1}+a_2\sum_{i=1}^m x_{i1}x_{i2}+···+a_n\sum_{i=1}^m x_{i1}x_{in}=-\sum_{i=1}^m x_{i1}y_{i}\\ a_1\sum_{i=1}^m x_{i2}x_{i1}+a_2\sum_{i=1}^m x_{i2} x_{i1}+···+a_n\sum_{i=1}^m x_{i2}x_{in}=-\sum_{i=1}^m x_{i2}y_{i}\\ ···\\ a_1\sum_{i=1}^m x_{in}x_{i1}+a_2\sum_{i=1}^m x_{in} x_{i1}+···+a_n\sum_{i=1}^m x_{in}x_{in}=-\sum_{i=1}^m x_{in}y_{i}$

上式已变成了一个关于a为变量的适定方程组，写成矩阵形式：

$\begin{pmatrix} \sum_{i=1}^m x_{i1}x_{i1}& \sum_{i=1}^m x_{i1}x_{i2}&···&\sum_{i=1}^m x_{i1}x_{in}\\ \sum_{i=1}^m x_{i2}x_{i1}& \sum_{i=1}^m x_{i2}x_{i2}&···&\sum_{i=1}^m x_{i2}x_{in}\\ ·····\\ \sum_{i=1}^m x_{in}x_{i1}& \sum_{i=1}^m x_{in}x_{i2}&···&\sum_{i=1}^m x_{in}x_{in}\\ \end{pmatrix} \begin{pmatrix} a_1 \\a_2\\···\\a_n \end{pmatrix}= \begin{pmatrix} -\sum_{i=1}^m x_{i1}y_{i} \\-\sum_{i=1}^m x_{i2}y_{i}\\···\\-\sum_{i=1}^m x_{in}y_{i} \end{pmatrix}$

整理上述矩阵，如下：

$\begin{pmatrix} \sum_{i=1}^m x_{i1}x_{i1}& \sum_{i=1}^m x_{i1}x_{i2}&···&\sum_{i=1}^m x_{i1}x_{in}\\ \sum_{i=1}^m x_{i2}x_{i1}& \sum_{i=1}^m x_{i2}x_{i2}&···&\sum_{i=1}^m x_{i2}x_{in}\\ ·····\\ \sum_{i=1}^m x_{in}x_{i1}& \sum_{i=1}^m x_{in}x_{i2}&···&\sum_{i=1}^m x_{in}x_{in}\\ \end{pmatrix} \begin{pmatrix} a_1 \\a_2\\···\\a_n \end{pmatrix}= \begin{pmatrix} x_{11}&x_{21}&···&x_{m1}\\ x_{12}&x_{22}&···&x_{m2}\\···\\ x_{1n}&x_{2n}&···&x_{mn} \end{pmatrix} \begin{pmatrix} y_1\\y_2\\···\\y_m \end{pmatrix}$

上式可以使用矩阵形式：

$X^TX\mathbf{a}=X^Ty$

最终解得a:

$\mathbf{a} = (A^T A )^{-1}A^T \mathbf{b}$

3. 线性方程组的最小二乘解-矩阵微积分工具

对于超定线性方程组中：

$Ax=y$

根据最小二乘思想，寻找一个向量x使得Ax-y的长度最小，为了简化数学计算将长度的平方作为目标，即

$S(x)=\|Ax-y\|^2$

展开后：

$S(x)=(Ax-y)^T(Ax-y)\\ \ \\ S(x)=(x^TA^T-y^T)(Ax-y)\\ \ \\ S(x)=x^TA^TAx-x^TA^Ty-y^TAx+y^Ty$

上面的 $x^TA^Ty$ 与 $y^TAx$ 都是标量（点积满足交换律），且这两者相等，上式可再次化简：

$S(x)=x^TA^TAx-2x^TA^Ty+y^Ty$

为寻找上式最小值，利用矩阵微积分工具，对x进行求导，并令导数等于0：

$\frac{\partial S(x)}{\partial x}=2A^TAx-2A^Ty=0\\ \ \\ A^TAx=A^Ty$

当 $A^TA$ 满秩时，此时方程组的解x为：

$x=(A^TA)^{-1}A^Ty$

4. 线性方程组的最小二乘解-向量空间推导

对于超定线性方程组：

$a_{1}x_{11}+a_{2}x_{12}+···+a_{n}x_{1n}=y_1\\ a_{1}x_{21}+a_{2}x_{22}+···+a_{n}x_{2n}=y_2\\ ···\\ a_{1}x_{m1}+a_{2}x_{m2}+···+a_{n}x_{mn}=y_m$

将上式用矩阵乘法的列观点重写：

$a_1\begin{pmatrix}x_{11}\\x_{21}\\···\\x_{m1}\end{pmatrix}+ a_2\begin{pmatrix}x_{12}\\x_{22}\\···\\x_{m2}\end{pmatrix}+ \begin{pmatrix} ···\\···\\···\\···\\\end{pmatrix}+ a_n\begin{pmatrix} x_{1n}\\x_{2n}\\···\\x_{mn} \end{pmatrix}=\begin{pmatrix} y_1\\y_2\\···\\y_m \end{pmatrix}$

上述意义就变成了寻找一组系数a使得向量x的线性组合最最接近向量y，在空间中向量y与向量Xa最接近的位置即为y向量向平台X投影位置，为了更好理解，这里使用二维超定线性方程组来理解推导，假设我们有三个方程，两个未知数：

$\begin{aligned} a_{11}x_1 + a_{12}x_2 &= b_1, \\ a_{21}x_1 + a_{22}x_2 &= b_2, \\ a_{31}x_1 + a_{32}x_2 &= b_3. \end{aligned}$

用向量表示，我们有：

$\mathbf{a}_1 = \begin{pmatrix} a_{11} \\ a_{21} \\ a_{31} \end{pmatrix}, \quad \mathbf{a}_2 = \begin{pmatrix} a_{12} \\ a_{22} \\ a_{32} \end{pmatrix}, \quad \mathbf{b} = \begin{pmatrix} b_1 \\ b_2 \\ b_3 \end{pmatrix}.$

我们想找到 x1 和 x2，使得线性组合 $x_1\mathbf{a}_1 + x_2\mathbf{a}_2$ 尽量接近 b。向量 a1 和 a2 位于三维空间 R3 中，这两个向量的线性组合 $x_1\mathbf{a}_1 + x_2\mathbf{a}_2$ 构成一个平面，这个平面称为矩阵A 的列空间。我们需要找到 x1 和 x2，使得向量 b到这个平面的距离最小，换句话说，我们要将 b投影到这个平面上，投影点即为 $x_1\mathbf{a}_1 + x_2\mathbf{a}_2$ 。投影的条件是： $\mathbf{r} = \mathbf{b} - (x_1\mathbf{a}_1 + x_2\mathbf{a}_2)$ ，要使 r 最小化，我们要求它与平面正交，如下图所示：

这意味着：

$\mathbf{r} \cdot \mathbf{a}_1 = 0 \quad \text{and} \quad \mathbf{r} \cdot \mathbf{a}_2 = 0$

展开后：

$\mathbf{b} \cdot \mathbf{a}_1 - x_1 (\mathbf{a}_1 \cdot \mathbf{a}_1) - x_2 (\mathbf{a}_2 \cdot \mathbf{a}_1) = 0\\ \mathbf{b} \cdot \mathbf{a}_2 - x_1 (\mathbf{a}_1 \cdot \mathbf{a}_2) - x_2 (\mathbf{a}_2 \cdot \mathbf{a}_2) = 0$

上式可以写成矩阵形式：

$\begin{pmatrix} \mathbf{a}_1 \cdot \mathbf{a}_1 & \mathbf{a}_2 \cdot \mathbf{a}_1 \\ \mathbf{a}_1 \cdot \mathbf{a}_2 & \mathbf{a}_2 \cdot \mathbf{a}_2 \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} \mathbf{b} \cdot \mathbf{a}_1 \\ \mathbf{b} \cdot \mathbf{a}_2 \end{pmatrix}$

整理后：

$A^T A \mathbf{x} = A^T \mathbf{b}$

最终解得x:

$\mathbf{x} = (A^T A )^{-1}A^T \mathbf{b}$

5. 总结

上面从三个角度推到了线性方程组最小二乘解，但实际中一般会使用matlab这类计算工具计算线性方程组的解，一般有2种方法来求解这类超定线性方程组：左除（x = A \ b）和伪逆（x = pinv(A) * b），另外在一些关注计算效率以及稳定性情况下可能采用矩阵分解来求逆，如LU分解、QR分解、Cholesky 分解。