对最小二乘法拟合直线的一些思考与改进

最小二乘法是我在工作中做图像处理最常用的一种拟合算法，同时也是在机器学习中线性回归的一种重要算法。我们先来回顾一下传统的最小二乘拟合直线模型算法：

如上图所示，为了从   五个点中拟合出最佳的直线方程，按照传统的做法，假设这些点符合线性关系   ,要拟合出最佳的直线方程，也就是让每个点在   处与   的偏差都很小，也即让总偏差的和最小，总结一句话就是如何找到合适的   与   ，让

最小，也即：

 但是这个式子带着绝对值，有正有负不太方便，所以我们一般让它们的平方最小，也即：

这就是最小二乘名字的由来（额，再鄙视一下这个蹩脚的翻译，为什么翻译成「最小二乘」，我觉得翻译成「最小平方」才能表达其本意）。

为了求出合适的   与   ，我们定义函数

根据多元函数求极值方法，我们把   与   看成未知量，并分别对他们求偏导，让它们等于0，得到：

化简，整理：

也即：

我们分别求出   与   的值即可得到我们想要的最「理想」的直线方程。

我们进行一次测试：下张图片是汽车电子的某个工业零件，为了测量左边的下半圆孔到右边斜边的距离，我们首先要拟合出那条斜边的直线方程，我们使用上面刚推导出来的斜截式公式进行拟合，代码如下(C++描述，下面雷同)：

by @yaxi_liu
void 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_liu
void 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;
}
}
}

我们来看下拟合的结果：

上图左侧绿线是我们拟合出来的结果，发现效果非常不错，而如果我们用斜截式来拟合的话，由于斜率无穷大，根本拟合不出来结果。

但是这里还有个问题，由于我们采用的是距离的平方最小进行拟合的，因此这个算法对离群的一些噪点是非常敏感的，尤其是对离群较大的点集合来说，它们的权重非常大，稍不注意就有可能把我们所需要的结果给拉偏了，例如下图，为了拟合出斜边方程，但是此斜边周围有一些常见的噪点，如果我们对这些噪点不加以处理，那么结果也就不太理想：

上图由于噪点的干扰，我们所需要的的直线明显被拉偏了。因此我们迫切需要对这些噪点进行处理，具体的做法就是对我们获取的这些点做权重分析，对离群较大的点给他们较小的权重，但这里涉及到权重大小的取值问题，具体做法我会在下篇文章中细致分析。

如果你想要本文的源代码，请关扫描下方二维码，注我的公众号，回复“最小二乘”四个字即可领取。