最小二乘法之加权最小二乘的应用

我在上文「对最小二乘法拟合直线的一些思考与改进」中介绍了通过距离的平方来拟合一条直线，这实际上也是一个PCA问题，但是我在结尾也说了，这种拟合方法的缺点是对离群噪点的干扰异常明显，稍不留神就把我们的目标给拉偏了，比如下面这张图片，在工业机器视觉检测中，为了测量某些产品的高精密尺寸，我们需要拟合出图中斜边边缘的准确直线方程，但是由于产品的缺陷，在边缘处往往有噪点的存在：

而如果我们不加以处理而直接拟合的话，最终的拟合结果是这样的：

绿线是我们根据上篇文章推导出来的公式进行拟合出来的结果，这显然不是我们所要的结果。由于噪点的干扰，我们所需的边缘直线方程被拉偏了。

这是由于在实际拟合过程中，离群较大的点它们都有很大的权重，而实际上每个点都应该有不同的权重，越偏离目标点它们的权重应该越小。因此我们迫切需要对每个点都定义一个权重，这就是今天我要介绍的加权最小二乘法。

在介绍这个算法之前，先回答一个问题，上篇文章中，有网友私信问我为什么那个参数方程要选取较小的特征值与特征向量：

这个问题我在上篇文中只是提到了一下，最小距离的总和刚好是数据集矩阵的特征值，但并没有说明具体的原因，这里给出一个更直观的证明方法：

如上图所示，我们定义一个直线的单位法向量  ,根据向量内积的几何意义，任意一点的坐标向量  与法向量  的内积刚好是这个点到直线的距离，但注意有方向正负，也即：

也即： 所以，最终 问题变成了，找到一个合适的单位法向量，让  

成立，注意，这里的  可根据我们上面推导出来的公式： 来进行计算，根据矩阵的相关计算，上式右边也即：

根据矩阵乘法结合律(矩阵乘法为什么满足结合律请参考我另外一篇文章「深入浅出线性代数的理解及应用」)，容易得知：  

为了方便，我们令 很明显，这里的运算实际上就是数据集的转置矩阵左乘其本身，计算结果是一个二阶方阵，也即： 注意，  是一个对称阵，因此它一定能够相似对角化，我们对其进行特征值分解：

其中  是  的正交阵，  是对角阵，且   ，这里  

因此：

我们令 由于  是正交阵，因此  的长度不变，也即 ，因此上式可化解为：  整理上式右边，得到：  由于约束条件中  ，因此上式中只有当  时才有最小值，最小值为  (较小的那个特征值)。也就是说，只有直线的法向量  的方向与数据集矩阵最小的特征值对应特征向量的方向一致的时候，距离平方之和才有最小值，最小值刚好为较小的那个特征值。这就是为什么要选取最小特征值的原因。

好了，我们再回到刚开始那个问题，我们究竟该对每个点如何定义权重值  呢？这是一个很有意思的问题，理论上，我们希望越偏离我们需要的目标值那么它的权值就应该越小，甚至为零，反之越大。按照这个理论，每个点的权重取值必须是基于点到直线的距离  (我们已经约束了拉格朗日乘子条件为  )来进行设定，因此权重一定是点到直线距离的函数。但是这里却有个问题，在没有拟合直线之前，我们并不知道每个点到直线的距离是多少，而为了拟合结果却必须知道每个点的距离，这看起来是一个矛盾的问题。解决这个问题的方法就是通过多次的迭代法来拟合直线。在第一次拟合过程中我们默认所有点的权重都为  ，然后默认它是正确的结果，再根据这个结果通过一个权重函数  来定义每个点的权重，然后将这些权重用于后续的迭代中。而权重函数的定义我们可以使用Huber权重函数：

参数  是我们定义的权重阈值，主要靠它来区分哪些是我们认为的离群噪点或者干扰点。在这个函数中，如果某点到直线的距离小于我们规定的值，那么它的权重为  ，否则它的权重会小于  ，而且会随着距离  的增加而减小。

但是我们究竟如何这个权重阈值  呢？一般情况下，我们可以通过这些数据集到直线距离的标准差来实现，而这个标准差可以通过下面这个公式来计算：这个式子可以用来计算离群值的标准偏差，右边的分子是数据点集合到直线距离的中值(注意是中值，而非均值)，分母的数字  究竟是如何来的，这里涉及到正态分布的知识，我有空会在下篇文章中详细介绍。标准差计算出来后，我们一般选取权值阈值是标准差的一个小的整数倍，例如.

那么我们该采用哪种方式进行拟合呢，上面说了，权重的取值跟点到直线的距离密切相关，因此我们依然采用点到直线距离的最小来拟合，这也是一个条件极值问题，同样需要适配拉格朗日乘子： 再分别对  求偏导，并让其等于  ，为了方便，我们直接套用上篇文章对最小二乘法拟合直线的一些思考与改进中推导出来的公式： 我们解出  的值，也就是我们需要的结果。 

下面给出实例代码(c++版)：

//函数名称：FitLineWeight                               
//返回类型：bool                         
//参数说明：vec:边缘点集合                                    
//功能：加权最小二乘拟合直线，拉格朗日乘数法
//作者：@刘亚曦
bool FitLineWeight(vector<POINT> vec, double &a, double &b, double &c)
{
  if (vec.size() < 2)
  {
    return false;
  }
  //首先进行第一次粗略计算
  FitLine(vec, a, b, c);

  //先默认结果是正确的，再计算权重
  for (size_t it = 0; it < times; it++)
  {
    if (vec.size() < 2)
    {
      return false;
    }
    double sum_x = 0.0f, sum_y = 0.0f;
    double avg_x = 0.0f, avg_y = 0.0f;
  
    //先计算集合点距离的中值
    vector<double> vecDistance;
    for (size_t i = 0; i < vec.size(); i++)
    {
      double d = abs(a * vec[i].x + b * vec[i].y + c);
      vecDistance.push_back(d);
    }
    //排序
    sort(vecDistance.begin(), vecDistance.end());

    //计算中值
    double median;
    if (vecDistance.size() % 2 == 0)
    {
      median = (vecDistance[vecDistance.size() / 2 - 1] + vecDistance[vecDistance.size() / 2]) / 2;
    }
    else
    {
      median = vecDistance[vecDistance.size() / 2];
    }
    //计算权重阈值分割点
    double t = 2 * (median / 0.6745);
  
    double sum_ww = 0.0f;
    vector<double> vecW;
    for (int i = 0; i < vec.size(); i++)
    {
      //计算本点权重
      double w;
      double d = abs(a * vec[i].x + b * vec[i].y + c);
      if (d <= t)
      {
        w = 1;
      }
      else
      {
        w = t / d;
      }

      vecW.push_back(w);
      sum_x += w * w * vec[i].x;
      sum_y += w * w * vec[i].y;
      sum_ww += w * w;
    }
    avg_x = sum_x / sum_ww;
    avg_y = sum_y / sum_ww;

    double L_xx = 0.0f;
    double L_yy = 0.0f;
    double L_xy = 0.0f;
    for (int i = 0; i < vec.size(); i++)
    {
      L_xx += pow(vecW[i], 2) * (vec[i].x - avg_x) * (vec[i].x - avg_x);
      L_xy += pow(vecW[i], 2) * (vec[i].x - avg_x) * (vec[i].y - avg_y);
      L_yy += pow(vecW[i], 2) * (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;
      }
    }
  }
}

按照上面的代码，分别调整代码里面的循环变量  ，从  开始，我们来试一下效果，还通过上面那张图片来说明：

不迭代结果(也就是没有权重分析，直接拟合，让变量  )，可以看到效果相当的差：

让  迭代一次的效果，比上张图片好一点，但还是有所偏离：

迭代两次又更好了：

迭代三次的结果：

可见，只经过了三次迭代，我们拟合的直线就收敛于我们想要的理想结果了。而事实也证明，基本上只要我们迭代三到五次，就能达到很好的效果。

说到这里，本文就算结束了，等我有空会再通过这种方法验证一下其他集合图形的拟合结果，例如圆或者椭圆，到时候再给大家回复。如果你想要本文的源代码，请扫描下面的图片，在我公众号内回复“加权最小二乘”六个字即可，我会将源代码都发送给你，方便你下去研究学习。