View Least Square Method (LSM) from Perspectives of Curve Fitting, Parameter Estimation, and Geometry Meaning of Solving Over-determined Equations
注：本文章仅仅关注线性系统，并未讨论更一般的非线性系统。
From perspective of curve fitting
假设有$N$个样本点：$D={(x_1,y_1), (x_2, y_2), \cdots, (x_N, y_N)}$，且 $x_i\in\mathbb{R}^p$，$y_i\in\mathbb{R}$，$i=1,2,\cdots, N$。设
\[\begin{equation} X=[x_1,x_2,\cdots, x_N]^T=\begin{bmatrix}x_1^T\\x_2^T\\\vdots\\x_N^T\end{bmatrix}_{N\times p}, \quad Y=\begin{bmatrix}y_1\\y_2\\\vdots\\y_N\end{bmatrix}_{N\times 1} \end{equation}\]假设有拟合曲线$f(w)=w^Tx$，则最小二乘估计（Least square estimation, LSE）的损失函数可以表示为：
\[L(w)=\sum_{i=1}^N|w^Tx_i-y_i|^2\]进一步展开可以得到：
\[\begin{align*}L(w)&=\sum_{i=1}^N|w^Tx_i-y_i|^2\\&=(w^Tx_1-y_1,w^Tx_2-y_2,\cdots,w^Tx_N-y_N) \begin{bmatrix}w^Tx_1-y_1\\w^Tx_2-y_2\\\vdots\\w^Tx_N-y_N\end{bmatrix}\\ &=(w^TX^T-Y^T)(Xw-Y)\\ &=w^TX^TXw-w^TX^TY-Y^TXw-Y^TY\\ &=w^TX^TXw-2w^TX^TY+Y^TY \end{align*}\]由LSM估计出的参数$\hat{w}$应当使该损失函数的值最小，因此有
\[\begin{aligned} \hat{w}&=\arg \min \limits_{w} L(w)\notag \\ &=\arg \min \limits_{w} \sum_{i=1}^N|w^Tx_i-y_i|^2\notag \\ &=\arg \min \limits_{w} \big[ w^TX^TXw-2w^TX^TY+Y^TY \big]\\ \end{aligned} \label{MLELoss}\] \[\begin{equation} \dfrac{\partial L(w)}{\partial w}=2X^TXw-2X^TY=0 \label{LS} \end{equation}\]因此参数$\hat{w}$的估计式为：
\[\hat{w}=(X^TX)^{-1}X^TY\label{1-2}\]其中，$X^{\dagger}=(X^TX)^{-1}X^T$称为$X$的伪逆（Pseudo Inverse）。实际上，矩阵$X^TX$也可能不可逆，这时我们需要直接求解$\eqref{LS}$所对应的线性方程组1。
From perspective of parameter estimation (Maximum Likelihood Estimation)
从参数估计角度讲，最小二乘法是解决这样一个问题：当我们已知噪声$\varepsilon$的分布和数据$x_i$和待估计的参数$w$，如何估计出参数$w$，使得$y_i$出现的概率最大。这正是极大似然估计法（Maximum Likelihood Estimation, MLE）的思想。
极大似然估计法（maximum likelihood estimation, MLE）2
设有总体分布$f(x;\theta_1, \cdots, \theta_k)$，$X_1, \cdots, X_n$为自这个总体中抽出的样本，则样本$(X_1, \cdots, X_n)$的分布，即其概率密度函数或概率函数）为：
\[L(x_1, \cdots, x_n;\theta_1, \cdots, \theta_k) = f(x_1;\theta_1, \cdots, \theta_k) f(x_2;\theta_1, \cdots, \theta_k) \cdots f(x_n;\theta_1, \cdots, \theta_k)\]对于函数$L(x_1, \cdots, x_n;\theta_1, \cdots, \theta_k)$：
- 若将函数$L$看作$x_1, \cdots, x_n$的函数时，$L$是一个概率密度或概率函数。如果$L(Y_1, \cdots, Y_n; \theta_1, \cdots, \theta_k) \gt L(X_1, \cdots, X_n;\theta_1, \cdots, \theta_k)$，则表明在观察时出现$(Y_1, \cdots, Y_n)$这个点的可能性比出现$X_1, \cdots, X_n$这个点的可能性大
- 反过来看，固定$X_1, \cdots, X_n$而把函数$L$看作$\theta_1, \cdots, \theta_k$的函数时，$L$被称为似然函数（Likelihood function），它的大小反映了在观察结果$(X_1, \cdots, X_n)$已知的条件下，$(\theta_1, \cdots, \theta_k)$的各种取值的“似然程度”。在这里，参数$\theta_1, \cdots, \theta_k$有一定的值（虽然未知），并非事件或者随机变量，无概率可言，故改用“似然”这个词。
根据MLE的定义，假设有噪声$\varepsilon\sim N(0, \sigma^2)$，并有拟合曲线$y=f(w)+\varepsilon=w^Tx+\varepsilon$，则
\[y\vert x_i,w \sim N(w^Tx,\sigma^2)=\dfrac1{\sqrt{2\pi}\sigma}\exp\{-\dfrac{(y-w^Tx)^2}{2\sigma^2}\}\]则$\log$似然函数为：
\[\begin{aligned} \mathcal{L}(w)&=\log P(Y|X,w)\\ &=\log\prod_{i=1}^NP(y_i|x_i,w)\\ &=\sum_{i=1}^N \log\dfrac1{\sqrt{2\pi}\sigma}-\dfrac1{2\sigma^2}(y_i-w^Tx_i)^2 \end{aligned}\]则根据极大似然估计法估计的参数$w$需要满足：
\[\begin{aligned} \hat{w}&=\arg \max \limits_{w} \mathcal{L}(w)\\ &=\arg \max \limits_{w}\sum_{i=1}^N-\dfrac1{2\sigma^2}(y_i-w^Tx_i)^2&\\ &=\arg \min \limits_{w}\sum_{i=1}^N \ (y_i-w^Tx_i)^2\\ \end{aligned}\]这与最小二乘估计法的公式$\eqref{MLELoss}$相吻合，该结果表明若假设噪声服从零均值的正态分布，MLE等价于LSM。
From perspective of geometry meaning of solving over-determined equations
对于线性方程组$A\boldsymbol{x}=\boldsymbol{b}$ ，如果方程个数$m$和未知变量个数$n$满足$m>n$，则称该线性方程组为超定方程组（也称为不一致方程组）。由于在实际工程中收集数据时常常存在误差，系数矩阵$A$和向量$\boldsymbol{b}$的值不够精确，导致超定方程组非常常见。
\[\begin{bmatrix}a&d\\b&e\\c&f\end{bmatrix} \begin{bmatrix}x\\y\end{bmatrix}=\begin{bmatrix}b_1\\b_2\\b_3\end{bmatrix}\]从几何上看，超定方程组试图寻找一个二维向量，使其经过给定线性变换后等于三维空间中给定的一个三维向量。我们知道，线性变换所对应矩阵的的列向量是经过变换后基向量的坐标，矩阵乘法就表示基向量的线性组合3，因此经过矩阵$A$所对应的线性变换后的空间中只具有两个基向量$\hat{i}=[a,b,c]$和$\hat{j}=[d,e,f]$，两个基向量只能张成一个二维空间（对应下图，基向量$\hat{i}$和$\hat{j}$张成绿色平面）。因此我们通常无法找到一个精确的二维向量，使其经过给定的线性变换后与给定的三维向量重合，即超定方程组没有精确解。
但是我们可以使用LSM计算一个近似解。从上图可以看到，我们可以在二维空间（紫色平面）中找到一个向量$\boldsymbol{\hat{x}}$，使其经过线性变换$A$后最接近空间中给定的向量$\boldsymbol{b}$ 。即向量$A\boldsymbol{\hat{x}}$具有以下特点：余项$\boldsymbol{b}-A\boldsymbol{\hat{x}}$和平面${A\boldsymbol{\hat{x}} \vert \boldsymbol{\hat{x}}\in \mathbb{R}^{n}}$垂直，即： \((\boldsymbol{b}-A\boldsymbol{\hat{x}})\bot \{A\boldsymbol{\hat{x}}\vert \boldsymbol{\hat{x}}\in \mathbb{R}^{n}\}\)
写成矩阵乘法有：
\[(A\boldsymbol{\hat{x}})^T(\boldsymbol{b}-A\boldsymbol{\hat{x}})=0\Rightarrow\boldsymbol{\hat{x}}^TA^T (\boldsymbol{b}-A\boldsymbol{\hat{x}})=0\label{3-1}\]式$\eqref{3-1}$表明$n$维向量$A^T (\boldsymbol{b}-A\boldsymbol{\hat{x}})$和$\mathbb{R}^{n}$中的其他的$n$维向量都垂直，也包含它自身（因为$\boldsymbol{\hat{x}}$是未知量，可任意取值），则只有一种可能的情况：
\[A^T (\boldsymbol{b}-A\boldsymbol{\hat{x}})=0\]即：
\[A^TA\boldsymbol{\hat{x}}=A^T\boldsymbol{b} \label{3-2}\]式$\eqref{3-2}$所定义的方程组有解，它的解$\boldsymbol{\hat{x}}$即是原方程组$A\boldsymbol{\hat{x}}=\boldsymbol{b}$的最小二乘解。式$\eqref{3-2}$也被称为原超定方程的法线方程。因此只要利用其他求解一致方程组的方法来求解该法线方程就得到解$\boldsymbol{\hat{x}}$即可。
\[\boldsymbol{\hat{x}}=(A^TA)^{-1}A^T\boldsymbol{b}\label{3-3}\]可以看出，公式$\eqref{3-3}$与由曲线拟合方法推出的的解析公式$\eqref{1-2}$是一致的。
若定义最小二乘解的余项为$\boldsymbol{r}=\boldsymbol{b}-A\boldsymbol{\hat{x}}$，并使用它的二范数的平方来衡量余项的大小：
\[\vert\vert\boldsymbol{r}\vert\vert^2_2=\vert\vert\boldsymbol{b}-A\boldsymbol{\hat{x}}\vert\vert_2^2=r_1^2+\cdots +r_m^2\label{3-4}\]式$\eqref{3-4}$就衡量了近似解$\boldsymbol{\hat{x}}$和精确解$\boldsymbol{\hat{x}}$之间的误差大小，最小化这个误差，就可以得到超定方程对应的最小二乘解。
\[\boldsymbol{\hat{x}}=\arg\min \limits_{\boldsymbol{\hat{x}}}||\boldsymbol{b}-A\boldsymbol{\hat{x}}||_2^2\label{3-5}\]公式$\eqref{3-5}$与公式$\eqref{MLELoss}$一致。
In closing
- 在噪声服从零均值的正态分布时，极大似然估计法等价于最小二乘法；
- 超定方程组(不一致方程组)没有精确解，但是可以使用最小二乘法求最接近原超定方程的最小二乘解。
References
- 3Blue1Brown线性代数系列视频.
- T. Sauer. 数值分析(原书第2版). 裴玉茹, 马赓宇, 译. 北京：机械工业出版社，2014.10.
- 【机器学习】【白板推导系列】P10.
-
陈希孺. 概率论与数理统计. 合肥: 中国科学技术大学出版社, 2009.2(2019.8重印). ˄