基本介绍

全部最小二乘解决形如如下的最小二乘问题

求解非零项转换为如下的特征值问题

可以最小化这个约数问题的就是最小的特征值对应的特征向量。它和的最末一个右特征向量是相同的。这是因为

通过选择最小的值，我们就可以达到约束这个问题的目的。
求解直线拟合问题的时候，为了所有的输入变量具有相同的噪声模型，需要先对变量进行处理，减去各自均值

然后通过最小化

来拟合2D直线方程。

最后在的。

BiSquare

给每个点施加一个权重，离群点的权重会变为零。
按matlab的文档中的写法

其中 为调谐常数， 为稳健标准偏差，由残差的中位绝对偏差 (MAD) 除以 得出。
第一次迭代的时候所有的权重都默认为1.
则在任意一次迭代中，直线的加权的最小全部最小二乘的形式变为

解的过程变为了

其中为对角矩阵，使得。
结果同样是取最小的值对应的来达到约束这个问题的目的。

实现步骤

• 1.通过加权最小二乘法拟合模型。对于第一次迭代，除非你指定权重，否则算法使用的权重等于 1。
  要先使用权重更新均值，然后是构造的求解矩阵:
  均值改为加权均值
  
  构造2D直线方程的加权最小完全平方公式的矩阵形式
• 2.计算调整后的残差并将其标准化。计算新的权重
  
  参考下图
  BiSquare

在实际操作过程中计算的过程为

• 3.如果拟合收敛，则退出迭代过程。否则，更新权重返回步骤1，执行下一次二平方权重拟合方法迭代。
  与使用稳健最小二乘拟合将异常值的影响降至最低相比，您可以将数据点标记为从拟合中剔除。

实现代码

int BiSquareLineFitting(
    const Pointf *vPts, const int n,
    double &xMean, double &yMean,
    float &a, float &b, float &c)
{
    double xsum{0.f}, ysum{0.f};

    // 初始化权重
    cv::Mat_<double> mWeight(n, 1);
    for (size_t i = 0; i < n; i++)
    {
        mWeight(i, 0) = 1.f;
    }

    cv::Mat_<double> mA(n, 2);

    // 迭代计算
    int maxIter = n * (n - 1) / 2;
    maxIter = min(maxIter, 50);

    cv::Mat_<double> mWA(n, 2);      ///< 实际计算时候，权重与A的乘积矩阵
    cv::Mat_<double> mR(n, 1);       ///< 偏差
    cv::Mat_<double> mRsort(n, 1);   ///< 偏差进行排序
    cv::Mat_<double> mMadsort(n, 1); ///< 偏差进行排序
    cv::Mat_<double> mWR(n, 1);      ///< 加权后的偏差
    cv::Mat w, u, vt;
    cv::Mat_<double> vk;
    double K = 4.685;
    double s{};
    double sumWR{}; // 加权后的偏差值
    for (size_t iter = 0; iter < maxIter; iter++)
    {
        // 这儿是不是可以根据权重先重新更新一下均值
        xsum = 0.0f;
        ysum = 0.0f;
        double w_iSum{0.f};
        for (size_t i = 0; i < n; i++)
        {
            xsum += vPts[i].x * mWeight(i, 0);
            ysum += vPts[i].y * mWeight(i, 0);
            w_iSum +=mWeight(i, 0);
        }
        // xMean = xsum / n;
        // yMean = ysum / n;

        xMean = xsum / w_iSum;
        yMean = ysum / w_iSum;

        for (size_t i = 0; i < n; i++)
        {
            mA(i, 0) = vPts[i].x - xMean;
            mA(i, 1) = vPts[i].y - yMean;
            // mA(i, 2) = 1.f;
        }

        // 更新计算权重后的矩阵
        for (size_t j = 0; j < n; j++)
        {
            mWA(j, 0) = mA(j, 0) * mWeight(j, 0);
            mWA(j, 1) = mA(j, 1) * mWeight(j, 0);
        }

        // 更新权重之后计算
        cv::SVD::compute(mWA, w, u, vt);
        vk = vt.t().col(1);
#ifdef _DEBUG
        std::cout << vk << std::endl;
#endif
        a = vk(0, 0);
        b = vk(1, 0);

        // 偏差矩阵
        for (size_t i = 0; i < n; i++)
        {
            // mR(i, 0) = mWeight(i, 0)*(a * mA(i, 0) + b * mA(i, 1));
            mR(i, 0) = (a * mA(i, 0) + b * mA(i, 1));
        }
        // 偏差中位数
        double median{};
        double wr_iSum{};
        w_iSum = 0.f;
        for (size_t i = 0; i < n; i++)
        {
            wr_iSum += mWeight(i, 0) * mR(i, 0);
            w_iSum += mWeight(i, 0);
        }
        median = wr_iSum / w_iSum;

        // mRsort = cv::abs(mRsort - median);
        mRsort = cv::abs(mR - median);
        cv::sort(mRsort, mMadsort, cv::SORT_EVERY_COLUMN + cv::SORT_ASCENDING);
        double mad = n % 2 == 0 ? (mMadsort(n / 2, 0) + mMadsort(n / 2 - 1, 0)) / 2 : mMadsort(n / 2, 0);
        s = mad / 0.6745;
        // 根据新计算出来的参数用bisquare更新权重
        for (size_t i = 0; i < n; i++)
        {
            double u_i = mR(i, 0) / (K * s);
            if (abs(u_i) < 1)
            {
                mWeight(i, 0) = pow(1 - pow(u_i, 2), 2);
            }
            else
            {
                mWeight(i, 0) = 0.f;
            }
            mWR(i, 0) = mWeight(i, 0) * pow(mR(i, 0), 2);
        }

        //
        sumWR = cv::sum(mWR)[0];
#ifdef _DEBUG
        // std::cout << "mWeight = " << mWeight << std::endl;
        std::cout << "sumWR = " << sumWR << std::endl;
#endif
    }

    a = vk(0, 0);
    b = vk(1, 0);
    c = -a * xMean - b * yMean;

    return 0;
}

参考

全部最小二乘
最小二乘法介绍（matlab文档翻译）