数值分析：多项式的最小二乘法拟合
多项式拟合-最小二乘法
最小二乘法使用多项式拟合样本数据，其基本代数原理以下。
假设使用的是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
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;
}
效果： 
参考链接：