最小二乘法拟合二维曲线的原理及实现

已知有n个数据点：$(x_1,y_1),(x_2,y_2),\dots ,(x_n,y_n)$，需要对这n个数据点进行曲线拟合，通过观察发现，它近似于抛物线，假定曲线为：$y=a_2\times x^2+a_1\times x+a_0$，其中$a_0、a_1、a_2$ 是未知的，如果把$(x_1,y_1)$代入方程，得到$y_1=a_2\times x_1^2+a_1\times x_1+a_0$，然后变形得：
$$\begin{array}{c}\left[\begin{array}{ccc}{x}_1^2& x_1& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]=y_1\end{array}$$
同理，将$(x_i,y_i),i=1,2,\dots ,n$, 代入方程，可以得到：
$$\begin{array}{c}\left[\begin{array}{ccc}{x}_n^2& x_n& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]=y_n\end{array}$$
组合成矩阵形式为：
$$\begin{array}{c}\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]=\left[\begin{array}{c}{y}_1\\ \dots \\ y_n\end{array}\right]\end{array}$$
等号两边的矩阵左边同时乘以：
$$\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\end{array}$$
可得：
$$\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]={\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{c}{y}_1\\ \dots \\ y_n\end{array}\right]\end{array}$$
等号两边的矩阵左边同时乘以：
$$\begin{array}{c}{\left[\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\end{array}\right]}^{-1}\end{array}$$
可得：
$$\begin{array}{c}{\left[\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\end{array}\right]}^{-1}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]={\left[\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\end{array}\right]}^{-1}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{c}{y}_1\\ \dots \\ y_n\end{array}\right]\end{array}$$
则：
$$\begin{array}{c}\left[\begin{array}{c}{a}_2\\ a_1\\ a_0\end{array}\right]={\left[\begin{array}{c}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]\end{array}\right]}^{-1}{\left[\begin{array}{ccc}{x}_1^2& x_1& 1\\ \dots & \dots & \dots \\ x_n^2& x_n& 1\end{array}\right]}^T\left[\begin{array}{c}{y}_1\\ \dots \\ y_n\end{array}\right]\end{array}$$

C++程序为：

#include<iostream>
#include<math.h>

using namespace std;

double a[][2] = {
					 1, 4.00,
					 2, 6.40,
					 3, 8.00,
					 4, 8.80,
					 5, 9.22,
					 6, 9.50,
					 7, 9.70,
					 8, 9.86,
					 9, 10.00,
					10, 10.20,
					11, 10.32,
					12, 10.42,
					13, 10.50,
					14, 10.55,
					15, 10.58,
					16, 10.60 
				};

void matrix_inverse(double(*a)[3], double(*b)[3]); // 声明矩阵求逆函数

int main()
{
	double a0[N][3] = {0};
	double a1[3][N] = {0};
	double a2[3][3] = {0};
	double a3[3][3] = {0};
	double a4[3][N] = {0};
	double b0[N][1] = {0};
	double c0[3][1] = {0};

	for(int i = 0; i < N; ++i)
	{
		a0[i][0] = pow(a[i][0], 2);
		a0[i][1] = a[i][0];
		a0[i][2] = 1;
		
		a1[0][i] = pow(a[i][0], 2);		// 	矩阵转置 
		a1[1][i] = a[i][0];
		a1[2][i] = 1;
		
		b0[i][0] = a[i][1];
	}

	for(int i = 0; i < 3; ++i)		// 矩阵相乘 (k的取值来源于左边矩阵的列数或右边矩阵的行数） 
		for(int j = 0; j < 3; ++j)
			for(int k = 0; k < N; ++k)
				a2[i][j] += (a1[i][k] * a0[k][j]);

	matrix_inverse(a2, a3);		// 矩阵求逆 

	for(int i = 0; i < 3; ++i)
		for(int j = 0; j < N; ++j)
			for(int k = 0; k < 3; ++k)
				a4[i][j] += (a3[i][k] * a1[k][j]);
	
	for(int i = 0; i < 3; ++i)
		for(int j = 0; j < 1; ++j)
			for(int k = 0; k < N; ++k)
				c0[i][j] += (a4[i][k] * b0[k][j]);
	
//	for(int i = 0; i < 3; ++i)
//	{
//		for(int j = 0; j < 1; ++j)
//		{
//			cout << " " << c0[i][j] << "\t";
//		}
//		cout << "\n";
//	}
	return 0;
 }
 
 void matrix_inverse(double(*a)[3], double(*b)[3])			// 矩阵求逆 
{
	using namespace std;
	int i, j, k;
	double max, temp;
	// 定义一个临时矩阵t
	double t[3][3];
	// 将a矩阵临时存放在矩阵t[n][n]中
	for (i = 0; i < 3; i++)
		for (j = 0; j < 3; j++)
			t[i][j] = a[i][j];
	// 初始化B矩阵为单位矩阵
	for (i = 0; i < 3; i++)
		for (j = 0; j < 3; j++)
			b[i][j] = (i == j) ? (double)1 : 0;
	// 进行列主消元，找到每一列的主元
	for (i = 0; i < 3; i++)
	{
		max = t[i][i];
		// 用于记录每一列中的第几个元素为主元
		k = i;
		// 寻找每一列中的主元元素
		for (j = i + 1; j < 3; j++)
		{
			if (fabs(t[j][i]) > fabs(max))
			{
				max = t[j][i];
				k = j;
			}
		}
		//cout<<"the max number is "<<max<<endl;
		// 如果主元所在的行不是第i行，则进行行交换
		if (k != i)
		{
			// 进行行交换
			for (j = 0; j < 3; j++)
			{
				temp = t[i][j];
				t[i][j] = t[k][j];
				t[k][j] = temp;
				// 伴随矩阵B也要进行行交换
				temp = b[i][j];
				b[i][j] = b[k][j];
				b[k][j] = temp;
			}
		}
		if (t[i][i] == 0)
		{
			cout << "\nthe matrix does not exist inverse matrix\n";
			break;
		}
		// 获取列主元素
		temp = t[i][i];
		// 将主元所在的行进行单位化处理
		//cout<<"\nthe temp is "<<temp<<endl;
		for (j = 0; j < 3; j++)
		{
			t[i][j] = t[i][j] / temp;
			b[i][j] = b[i][j] / temp;
		}
		for (j = 0; j < 3; j++)
		{
			if (j != i)
			{
				temp = t[j][i];
				//消去该列的其他元素
				for (k = 0; k < 3; k++)
				{
					t[j][k] = t[j][k] - temp * t[i][k];
					b[j][k] = b[j][k] - temp * b[i][k];
				}
			}
		}
	}
}