对最小二乘法拟合直线的一些思考与改进
最小二乘法是我在工作中做图像处理最常用的一种拟合算法，同时也是在机器学习中线性回归的一种重要算法。我们先来回顾一下传统的最小二乘拟合直线模型算法：
如上图所示，为了从 五个点中拟合出最佳的直线方程，按照传统的做法，假设这些点符合线性关系 ,要拟合出最佳的直线方程，也就是让每个点在 处与 的偏差都很小，也即让总偏差的和最小，总结一句话就是如何找到合适的 与 ，让
最小，也即：
但是这个式子带着绝对值，有正有负不太方便，所以我们一般让它们的平方最小，也即：
这就是最小二乘名字的由来（额，再鄙视一下这个蹩脚的翻译，为什么翻译成「最小二乘」，我觉得翻译成「最小平方」才能表达其本意）。
为了求出合适的 与 ，我们定义函数
根据多元函数求极值方法，我们把 与 看成未知量，并分别对他们求偏导，让它们等于0，得到：
化简，整理：
也即：
我们分别求出 与 的值即可得到我们想要的最「理想」的直线方程。
我们进行一次测试：下张图片是汽车电子的某个工业零件，为了测量左边的下半圆孔到右边斜边的距离，我们首先要拟合出那条斜边的直线方程，我们使用上面刚推导出来的斜截式公式进行拟合，代码如下(C++描述，下面雷同)：
by @yaxi_liuvoid FitLine(vector<POINT> vec, double &k, double &b){ double sum_xx = 0.0f; double sum_x = 0.0f; double sum_y = 0.0f; double sum_xy = 0.0f; int N = vec.size(); for (int i = 0; i < N; i++) { sum_xx += vec[i].x*vec[i].x; sum_x += vec[i].x; sum_y += vec[i].y; sum_xy += vec[i].x * vec[i].y; } k=(sum_xy*N-sum_x*sum_y) (sum_xx*N-sum_x*sum_x); b =(sum_xx*sum_y-sum_xy*sum_x) (sum_xx*N-sum_x*sum_x);}拟合结果如下：
绿线是我们根据拟合出来的直线方程结果，可见拟合出来的结果还算不错。
但是，在这里不知道大家有没有产生过这样的疑问：
第一，为什么我们想要让纵向的「距离」最小，而不是垂直的「距离」最小，后者不是更直观吗？
第二，当我们需要拟合的直线是垂直线或者接近垂直线(在图像处理中我经常遇到这种情况)，这种方法可行吗？
对于第一个问题，有一种的解释是，在实际工程应用中， 一般是没有误差的固定值，而只有 才有误差值，所以纵向的距离更有意义一些。
而第二个问题，答案是显而易见的，由于这个回归模型的直线方程 是本身不能表达垂直方程的，因此对于垂直线无能为力。所以我们不禁又想到了第一个问题：如果使用经典模型 ，让每个点到这条直线的垂直距离为最小，这样拟合出来的效果怎么样呢？
对于解决这个问题，与上面一样，基本框架思路是一样的，只是我们需要把问题简化成求距离的最小，也即我们需要找到合适的a,b,c使得
成立.
那么我们首先需要定义直线方程 ，但是这实际上已经过参数化了，我们需要对a,b加以限制，至少不能让他们同时为0，所以我们添加限制条件： 之所以这样限制，是为了在我们求距离的分母直接变为1，在计算的时候非常方面。因此，上面的问题等价于在上面的条件限制下，如何让下面的距离平方和最小：
很明显，这是一个条件极值问题，为了解决这个问题，我们需要对上述的限制条件作为格朗日乘子加以约束，再定义一个函数
我们按照惯例，对 分别求偏导，并让他们等于0，得到：
我们将上方程组的第三个计算的结果：
其中 分别为 的均值，将其代入方程组的1与2中，得到：
为了方便，我们令 ，得到：
上式化简为：
注意到左边的矩阵是一个对称阵，因此它一定有对应的特征值与特征向量，但是这个矩阵明显有两个不同的特征值，我们到底该取哪一个呢？注意观察，这里有一个很有意思的巧合，我们为了让
最小，我们再对上述方程组左乘特征向量的转置矩阵 ,得到的结果刚好是我们需要的最小值，结果也恰好是这个特征值，因此我们选择比较小的那个特征向量：
我们按照这个新公式再次进行测试，下图片还是上面某个汽车电子的零部件，我们为了求出椭圆孔左侧近乎垂直的直线方程，实现代码如下(C++描述）：
by @yaxi_liuvoid FitLine(vector<POINT> vec, double &a, double &b, double &c){ if (vec.size() < 2) { return; } double sum_x=0.0f, sum_y=0.0f; double avg_x = 0.0f, avg_y = 0.0f; for (int i = 0; i < vec.size(); i++) { sum_x += vec[i].x; sum_y += vec[i].y; } avg_x = sum_x / vec.size(); avg_y = sum_y / vec.size(); double L_xx = 0.0f; double L_yy = 0.0f; double L_xy = 0.0f; int N = vec.size(); for (int i = 0; i <N; i++) { L_xx += (vec[i].x - avg_x) * (vec[i].x - avg_x); L_xy += (vec[i].x - avg_x) * (vec[i].y - avg_y); L_yy += (vec[i].y - avg_y) * (vec[i].y - avg_y); } double lamd1, lamd2; double lamd; double m, n; m = L_xx + L_yy; n = L_xx*L_yy-L_xy*L_xy; lamd1=(m+sqrt(m*m-4*n))/2; lamd2=(m-sqrt(m*m-4*n))/2; lamd=lamd1<lamd2?lamd1:lamd2; double d=sqrt((L_xx-lamd)* (L_xx-lamd)+L_xy*L_xy); if (abs(d) < 1e-6) { a = 1; b = 0; c = -a * avg_x - b * avg_y; } else { if (lamd >= L_xx) { a = L_xy / d; b = (lamd - L_xx) / d; c = -a * avg_x - b * avg_y; } else { a = -L_xy / d; b = (L_xx - lamd) / d; c = -a * avg_x - b * avg_y; } }}我们来看下拟合的结果：
上图左侧绿线是我们拟合出来的结果，发现效果非常不错，而如果我们用斜截式来拟合的话，由于斜率无穷大，根本拟合不出来结果。
但是这里还有个问题，由于我们采用的是距离的平方最小进行拟合的，因此这个算法对离群的一些噪点是非常敏感的，尤其是对离群较大的点集合来说，它们的权重非常大，稍不注意就有可能把我们所需要的结果给拉偏了，例如下图，为了拟合出斜边方程，但是此斜边周围有一些常见的噪点，如果我们对这些噪点不加以处理，那么结果也就不太理想：
上图由于噪点的干扰，我们所需要的的直线明显被拉偏了。因此我们迫切需要对这些噪点进行处理，具体的做法就是对我们获取的这些点做权重分析，对离群较大的点给他们较小的权重，但这里涉及到权重大小的取值问题，具体做法我会在下篇文章中细致分析。
如果你想要本文的源代码，请关扫描下方二维码，注我的公众号，回复“最小二乘”四个字即可领取。