最小二乘法求解直线拟合问题

本文瀏覽次數

发布日期

2022年1月4日

摘要
本文以直线拟合问题为例,介绍了最小二乘法的计算方法。文章基于SVD分解和拉格朗日乘数法,给出了直线拟合问题的闭式解。

1801年1月1日,意大利天文学家Giuseppe Piazzi发现了谷神星。在连续观测了几十天后,谷神星湮没在太阳的眩光之中。为了再次找到它,24岁的高斯凭借现在被称为最小二乘法的方法,预测了谷神星的轨道。1801年12月31日,天文学家在高斯预测的轨道附近再度找回了谷神星。这就是最小二乘法在历史上首次登场的故事。

最小二乘(least squares)的含义是最小化误差的平方和(该术语受到日文影响。其中的“二乘”即为“平方”)。给定一系列样本点,以最小二乘为目标,求模型参数的优化方法叫做最小二乘法。其中有线性解的方法叫做线性最小二乘法

下面的小程序给出了一个例子。你可以通过鼠标点击在方框中添加二维点。如果方框中有两个以上的点,程序就会应用线性最小二乘法,找到使误差平方和最小的直线。这个例子中,样本点是我们手动添加的二维坐标,模型参数是直线的系数。下文将会详细介绍程序的计算原理。

浏览器不支持!!

直线拟合问题

上面展示的是一个求解直线拟合问题的程序。我们先形式化地描述这个问题,然后给出数学求解步骤。

问题:已知二维平面上的有一系列点 \((x_k, y_k)\)\(k=1,2,3,\cdots\)要求用最小二乘法拟合直线。

解: 设直线方程为\(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)=0 \]
方程满足\(A\v{x}=0\)的形式。其中\(A\)\(k\times 3\)的矩阵,\(x=(a, b, c)^T\).

所谓最小二乘法,即是求解\(x\)使得 \[||Ax||^2=(Ax)^T Ax=x^T A^T A x\] 的值最小。

这个问题显然有一个解\(x=0\),但这个解没有意义,因为它不能构成直线的系数,人们把这个解叫做“平凡解(trivial solution)”。

为了排除平凡解,引入一个约束条件\[x^Tx=1\]

问题变为约束条件下求最小值的问题,可以运用拉格朗日乘数法。令 \[L(x)=x^TA^TAx-\lambda(x^Tx-1)\]

\(L(x)\)最小时,\[{\partial L(x)\over \partial x} = 0,\] 因此 \[A^TAx - \lambda x = 0\]

因此\(\lambda\)\(A^TA\)的特征值,\(x\)\(A^T A\)对应\(\lambda\)的一个单位特征向量,记\(x=\boldsymbol{e_\lambda}\)。代入得\(L(\boldsymbol{e_\lambda})=\lambda \boldsymbol{e_\lambda}^T \boldsymbol{e_\lambda}=\lambda\)\(\lambda\)最小时,\(L(x)\)取得最小值。换言之,\(A^TA\)的最小特征值对应的单位特征向量即是所求的解。

求矩阵的特征特征向量有几种方法。注意这里的\(A^T A\)是实对称矩阵,可以用Jacobi eigenvalue algorithm求特征值。另一种常见的方法则是使用SVD分解。

用SVD分解求解最小二乘问题

根据SVD分解, \[A=UDV^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] \] 注意到\(A^TA=(UDV^T)^T(UDV^T)=VDU^TUDV^T=V\Lambda V^T\) 其中,\(\Lambda\)的第i行的非0元素\(\lambda_i=\sigma_i^2\) 两边右乘\(V\)得到 \[(A^TA)V=V\Lambda\] 即是说,\(V\)的每一个列向量\(\boldsymbol{v_j}\)都是\(A^TA\)的特征向量,对应的特征值是\(\lambda_j\).

总结

为了解决文中所述的带平凡解的最小二乘问题,本文先引入约束条件\(x^Tx=1\),然后运用拉格朗日乘数法来求解最小值。实际上,还有其它种形式的约束条件,例如,可以约定\(b=1\),但这样一来, 直线方程从一般式退化成斜截式,就无法表示平行于y轴的直线了。尽管有缺陷,这种形式的约束条件也还是很常用。

若是采用斜截式,经过求导等数学推导步骤,一样可以求出直线的参数。但不论采用何种约束条件,这两种方法都是以最小化误差的平方和为目的的,都可以叫做最小二乘法。

最小二乘法不仅可以用于直线的拟合,也可以推广到圆、椭圆、二项式、平面……等各类模型的拟合。

参考资料

By @執迷 in
Tags : #最小二乘法,