\[ \def \v#1{{\boldsymbol{#1}}} \def \m#1{{\mathbf{#1}}} \def \argmax#1{\underset{#1}{\operatorname{argmax}}} \]
1801年1月1日，意大利天文学家 Giuseppe Piazzi 发现了谷神星。在连续观测了几十天后，谷神星湮没在太阳的眩光之中。为了再次找到它，24岁的高斯凭借现在被称为最小二乘法的方法，预测了谷神星的轨道。1801年12月31日，天文学家在高斯预测的轨道附近再度找回了谷神星。这就是最小二乘法在历史上首次登场的故事。
最小二乘（least squares）的含义是最小化误差的平方和（该术语受到日文影响。其中的”二乘”即为”平方”）。给定一系列样本点，以最小二乘为目标，求模型参数的优化方法叫做最小二乘法。其中有线性解的方法叫做线性最小二乘法。
下面的小程序给出了一个例子。你可以通过鼠标点击在方框中添加二维点。如果方框中有两个以上的点，程序就会应用线性最小二乘法，找到使误差平方和最小的直线。这个例子中，样本点是我们手动添加的二维坐标，模型参数是直线的系数。下文将会详细介绍程序的计算原理。
最小二乘法求解直线拟合问题
上面展示的是一个求解直线拟合问题的程序。我们先形式化地描述这个问题，然后给出数学求解步骤。
问题： 已知二维平面上的有一系列点 \((x_k, y_k)\)，\(k=1,2,3,\cdots\)，要求用最小二乘法拟合直线。
针对这个问题，根据设置的目标函数不同，解法也不同。下面将给出几种不同的目标函数，说明其特点，然后给出对应的解法。
斜截式
在这个例子中，我们设直线方程为 \(y=ax+b\)，其中 \(a\) 为斜率，\(b\) 为截距。这种类型的直线方程被称为斜截式。
根据方程可列出如下式子：
\[ \left(\begin{matrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \end{matrix} \right) \left(\begin{matrix} a \\ b \\ \end{matrix}\right) = \left(\begin{matrix} y_1\\ y_2\\ \vdots \end{matrix}\right) \]
该方程满足 \(\mathbf A \v x = \v b\) 的形式，其中 \(\mathbf A\in \mathbb R^{k\times 2}\)，\(\v x=(a,b)^T\), \(\v b \in \mathbb R^{k\times 1}\)。
解：
我们的目标为最小化
\[ \begin{aligned} L(\v x) &= ||\mathbf A\v x - \v b||^2 \\ &= \v x^T \mathbf A^T \mathbf A \v x - 2 \v x ^T \mathbf A^T \v b +\v b^T \v b \end{aligned} \]
令其导数为 \(0\)
\[ \frac{\partial (L(\v x))}{\partial x} = 0 \]
得
\[ \begin{aligned} &\mathbf A^T \mathbf A \v x - \mathbf A^T \v b = 0 \\ \Rightarrow &\mathbf A^T \mathbf A \v x = \mathbf A^T \v b\\ \Rightarrow &\v x= (\mathbf A^T \mathbf A)^{-1} \mathbf A^T \v b \end{aligned} \]
这样就得到了该问题的最优解。
在这个解法中，矩阵 \((\mathbf A^T \mathbf A)^{-1} \mathbf A^T\) 是 \(\mathbf A\) 的伪逆矩阵。这是因为 \((\mathbf A^T \mathbf A)^{-1} \mathbf A^T \mathbf A= \mathbf I\)。
\((\mathbf A^T \mathbf A)^{-1} \mathbf A^T\) 出现在 \(\mathbf A\) 的左边，被称为 \(\mathbf A\) 的左逆。
斜截式的局限性在于其不能表示垂直于 \(x\) 轴的直线。下面的方法能够解决这个问题。
一般式
设直线方程为 \(ax+by+c=0\)，则有
\[ \left(\begin{matrix} x_1 & y_1 & 1\\ x_2 & y_2 & 1\\ \vdots & \vdots & \vdots\end{matrix}\right) \left(\begin{matrix} a\\b\\c \end{matrix}\right)=\v 0 \]
方程满足 \(\mathbf A\v{x}=\v 0\) 的形式。其中 \(\mathbf A\) 是 \(k\times 3\) 的矩阵，\(\v x=(a, b, c)^T\)。
解：
与开头的方法不同，此方法采用一般式来表示直线，因此能够正常表示垂直于 \(x\) 轴的直线。
根据给定的目标，求解 \(\v x\) 使得
\[ ||\mathbf A\v x||^2=(\mathbf A \v x)^T\mathbf A\v x=\v x^T \mathbf A^T \mathbf A \v x \]
的值最小。
这个问题显然有一个解 \(\v x=\v 0\)，但这个解没有意义，因为它不能构成直线的系数，人们把这个解叫做”平凡解（trivial solution）“。
为了排除平凡解，引入一个约束条件：
\[ \v x^T \v x=1 \]
问题变为约束条件下求最小值的问题，可以运用拉格朗日乘数法。令
\[ L(\v x)=\v x^T\mathbf A^T\mathbf A\v x-\lambda(\v x^T\v x-1) \]
\(L(\v x)\) 最小时，
\[ {\partial L(\v x)\over \partial \v x} = \v 0, \]
因此
\[ \mathbf A^T\mathbf A\v x - \lambda \v x = \v 0 \]
因此 \(\lambda\) 是 \(\mathbf A^T\mathbf A\) 的特征值，\(\v x\) 是 \(\mathbf A^T \mathbf A\) 对应 \(\lambda\) 的一个单位特征向量，记 \(\v x=\boldsymbol{e_\lambda}\)。代入得 \(L(\boldsymbol{e_\lambda})=\lambda \boldsymbol{e_\lambda}^T \boldsymbol{e_\lambda}=\lambda\)。\(\lambda\) 最小时，\(L(\v x)\) 取得最小值。换言之，\(\mathbf A^T\mathbf A\) 的最小特征值对应的单位特征向量即是所求的解。
求矩阵的特征向量有几种方法。注意这里的 \(\mathbf A^T\mathbf A\) 是实对称矩阵，可以用 Jacobi eigenvalue algorithm 求特征值。另一种常见的方法则是使用 SVD 分解。
用 SVD 分解求解最小二乘问题
根据 SVD 分解，
\[ \m A=\m U\m D\m V^T= \left[\begin{matrix} \boldsymbol{u_0} | \cdots | \boldsymbol{u_{p-1}} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{\sigma_0} & & \\ & \ddots & \\ & & \boldsymbol{\sigma_{p-1}} \end{matrix}\right] \left[\begin{matrix} \boldsymbol{v_0^T}\\ \hline \vdots\\ \hline \boldsymbol{v_{p-1}^T} \end{matrix}\right] \]
注意到
\[ \m A^T\m A=(\m U\m D\m V^T)^T(\m U\m D\m V^T)=\m V\m D\m U^T\m U\m D\m V^T=\m V\m \Lambda \m V^T \]
其中，\(\m \Lambda\) 的第 \(i\) 行的非 0 元素 \(\lambda_i=\sigma_i^2\)。两边右乘 \(\m V\) 得到
\[ (\m A^T\m A)\m V=\m V\m \Lambda \]
即是说，\(\m V\) 的每一个列向量 \(\boldsymbol{v_j}\) 都是 \(\m A^T\m A\) 的特征向量，对应的特征值是 \(\lambda_j\)。
局限性
前文中我们分斜截式和一般式两种情况讨论了直线拟合问题，并指出一般式具有能够拟合任意方向直线的优点。但还应当注意到两种方法都存在有共同的局限性：其并非以最小化点到直线的垂直距离为目标，即最小化
\[ \sum_i {|ax_i + by_i + c|\over \sqrt{a^2 + b^2}} \]
“正交回归”方法可以解决这一问题。但正交回归已经超出了本文要讨论的范畴，这里就不再继续展开。
总结
为了解决文中所述的带平凡解的最小二乘问题，本文先使用斜截式建立直线方程并给出解，然后指出了斜截式的局限性。
而后本文尝试了使用一般式表示直线，并引入约束条件 \(\v x^T \v x=1\)，然后运用拉格朗日乘数法来求解最小值。
最后，文章讨论了此类方法的局限性。
虽然本文围绕著直线拟合问题展开，但最小二乘法不仅可以用于直线的拟合，也可以推广到圆、椭圆、二项式、平面……等各类模型，是一个非常实用的数学工具。
