多项式拟合-最小二乘法

代数角度

最小二乘法使用多项式拟合样本数据，其基本代数原理以下。

假设使用的是m阶多项式拟合，那么在的估计值为：

其误差表示为：

那么对于n个样本点，使用平方和误差量化，该误差就是数据的优化目标，自变量为各阶的权值w，有：

求E(w)极小值，需对其求偏导，为：

有:

权重系数与i无关，可以写成：

同样的，此为一个约束条件，对于m阶的多项式，存在m+1项的偏导联合约束，所以有：

写成矩阵表达式:

C++实现

如下： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
#include <opencv2/opencv.hpp>
#include <iostream>
#include <math.h>

using namespace std;
using namespace cv;

//样本点集、阶数，输出拟合点
std::vector<cv::Point> PolyFit(const std::vector<cv::Point>& points, int m){
    cv::Mat A = cv::Mat::zeros(m+1, m+1, CV_32FC1);
    cv::Mat B = cv::Mat::zeros(m+1,1, CV_32FC1);

    int n = points.size();

    for(int i=0; i<m+1; i++){
        for(int j=0; j<m+1; j++){
            for(int k=0; k<n; k++){
                A.at<float>(i, j) += std::pow(points[k].x, i+j);
            }
        }
    }

    for(int i=0; i<m+1; i++){
        for(int k=0; k<n; k++)
            B.at<float>(i,0) += points[k].y * std::pow(points[k].x, i);
    }


    cv::Mat W = cv::Mat::zeros(m+1, 1, CV_32FC1); //解出权重矩阵
    if(!cv::solve(A,B,W, cv::DECOMP_LU)){
        cout << "Failed To Solve" << endl;
        return {};
    }

    //计算拟合点
    vector<cv::Point> result(n, cv::Point(0,0));
    for(int i=0; i<n; i++){
        result[i].x = points[i].x;
        for(int j=0; j<m+1; j++){
            result[i].y += W.at<float>(j,0) * std::pow(points[i].x, j);
        }
    }
    return result;
}

int main(int argc, char **argv) {
    const int height = 600;
    const int width = 1000;
    cv::Mat canvas = cv::Mat(height, width, CV_8UC3, Scalar(255,255,255));

    vector<cv::Point> points;
    for(int i=0; i<width; i+=10){ 
        int x = i;
        double noise = theRNG().uniform(-30.0, 30.0);
        int y = std::round((double)(height - ((double)height/width) * x) + noise);
        points.emplace_back(cv::Point(x,y));
        cv::circle(canvas, cv::Point(x,y), 5, Scalar(255,0,0), -1);
    }

    vector<cv::Point>fitPoints = PolyFit(points, 3); //三阶拟合
    cv::polylines(canvas, fitPoints, false, Scalar(0,0,255), 4);

    namedWindow("canvas", cv::WINDOW_NORMAL);
    cv::imshow("canvas",canvas); 

    cout << "Done" << endl;
    waitKey(0);
    return 0;
}

效果： 最小二乘法拟合

矩阵角度看待最小二乘求解

除了代数推导，从矩阵维度描述的最小二乘法有更加简便的求解方式，主要是正规方程求解以及QR分解/SVD分解等。

正规矩阵求解

假设有m阶拟合，其矩阵可以表示成：

其估计误差可以表示成一个列向量：

于是其平方误差可以表示为：

注意，此处的E已经是一个标量而非矩阵了。

对W求导：

  1. 来自前述规律，A为对称矩阵情况下

  2. E是标量，每一项和式也应该是标量，因此其和其转置都应该相等：

所以其中两项求导均是

于是：

至此就得到了最小二乘法的正规方程：

构造该方程，求解W即可，比较简单。

C++实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
///正规方程方法
//输入样本点集、阶数，输出拟合点
std::vector<cv::Point> PolyFit_Matrix(const std::vector<cv::Point>& points, int m){
    int n = points.size();
    cv::Mat A = cv::Mat::zeros(n, m+1, CV_32FC1);
    cv::Mat W = cv::Mat::zeros(m+1, 1, CV_32FC1); //待求权重矩阵
    cv::Mat B = cv::Mat::zeros(n, 1, CV_32FC1);


    for(int i=0; i<n; i++){
        for(int j=0; j<m+1; j++){
            A.at<float>(i,j) = std::pow(points.at(i).x, j);
        }
    }

    for(int i=0; i<n; i++){
        B.at<float>(i, 0) = points.at(i).y;
    }

    cv::Mat At = A.t();

    cv::solve(At*A, At*B, W, cv::DECOMP_NORMAL);

    //计算拟合点
    vector<cv::Point> result(n, cv::Point(0,0));
    for(int i=0; i<n; i++){
        result[i].x = points[i].x;
        for(int j=0; j<m+1; j++){
            result[i].y += W.at<float>(j,0) * std::pow(points[i].x, j);
        }
    }
    return result;
}

QR分解求解

LU分解/QR分解/SVD分解会在之后的文章详细叙述，此处仅给出结论和实现。使用正规方程求解的一个缺陷是需要计算作为左侧矩阵方程，但是的条件数是A矩阵条件数的平方，一个矩阵的条件数可以定义为: 条件数是衡量一个矩阵的“健康系数”，例如，一些线性高度相关的列组成的向量，其条件数一般很大，表现在其在空间的夹角非常小： 如果使用该矩阵用于求解： 那么只要系数出现微小的扰动（例如上面的1变成1.1）,W系数可能变化很大，这种现象被称为数值的不稳定性，这样的系数矩阵，一般被称为病态矩阵。

所以为求稳定的数值，必须使用那些绕过直接计算而求解线性方程的方法，对于这样的最小二乘问题，最常用的是QR分解，其只需要相同的构造，使用QR方法求解即可，反而是三种方法下最简单的： 当然SVD分解也是可行的，但其计算量稍大，稳定性也更好。

C++实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
//输入样本点集、阶数，输出拟合点
std::vector<cv::Point> PolyFit_QR(const std::vector<cv::Point>& points, int m){   ///QR分解/SVD分解方法
    int n = points.size();
    cv::Mat A = cv::Mat::zeros(n, m+1, CV_32FC1);
    cv::Mat W = cv::Mat::zeros(m+1, 1, CV_32FC1); //待求权重矩阵
    cv::Mat B = cv::Mat::zeros(n, 1, CV_32FC1);

    for(int i=0; i<n; i++){
        for(int j=0; j<m+1; j++){
            A.at<float>(i,j) = std::pow(points.at(i).x, j);
        }
    }

    for(int i=0; i<n; i++){
        B.at<float>(i, 0) = points.at(i).y;
    }

    cv::solve(A, B, W, cv::DECOMP_QR);  //或SVD

    //计算拟合点
    vector<cv::Point> result(n, cv::Point(0,0));
    for(int i=0; i<n; i++){
        result[i].x = points[i].x;
        for(int j=0; j<m+1; j++){
            result[i].y += W.at<float>(j,0) * std::pow(points[i].x, j);
        }
    }
    return result;
}

最小二乘求旋转角

待续

参考链接：

基于最小二乘法的多项式曲线拟合：从原理到c++实现