最小二乘法的多项式拟合
给定一系列点,运用最小二乘法拟合出这些点的多项式曲线函数。假设拟合出的多项式(k项)函数为:
又假设真实曲线函数为:
则对于给定的n个点(xi,yi ),0<i≤n,拟合曲线与实际偏差要最小,即:
则问题转换为:求一组解(a0,a1,a2,…,ak)^T,使得函数L的值最小。由极值条件可知,让L取最小值时的一组解,L对其的偏导数为0;为了求解,对L求ai的偏导并令其为0:
化简可得:
表示成矩阵为:
明显上式是一个解线性方程组问题:A*X=B。
在编程时,对于A的获取不必计算(k+1)^2次,可以发现A中的元素仅有2k+1种:
因此编程时仅需先算出这2k+1项,再将其按规律填入A即可。下面的示例代码演示了用Eigen库的Qr分解算法求解A*X=B问题。
#include <Eigen/Dense>
//--Input
//---x:数据点的x坐标值数组
//---y:数据点的y坐标值数组
//---n:数据点的个数
//---k:多项式的阶
//--Output
//---返回A*X=B中的解X
Eigen::VectorXd polynomialFitting(double*x, double*y, unsigned n, unsigned k)
{
Eigen::MatrixXd A=Eigen::MatrixXd::Zero(k+1,k+1);
Eigen::VectorXd B=Eigen::VectorXd::Zero(k+1);
double*x_pow_sum=new double[2*k+1];//存放A中2k+1项不重复元素
for(unsigned i=0;i<2*k+1;i++)
{
double sum=0.;
for(unsigned j=0;j<n;j++)
{
sum+=pow(x[j],i);//计算每项幂的和
}
x_pow_sum[i]=sum;//
}
for(unsigned i=0;i<k+1;i++)
{
for(unsigned j=0;j<k+1;j++)
{
A(i,j)=x_pow_sum[i+j];//按规律将2k+1项不重复元素填入矩阵A
}
}
for(unsigned i=0;i<k+1;i++)
{
double sum=0;
for(unsigned j=0;j<n;j++)
{
sum=sum+pow(x[j],i)*y[j];//计算B的每个元素
}
B[i]=sum;
}
//利用Eigen库中的Qr分解求解X
Eigen::VectorXd X=A.colPivHouseholderQr().solve(B);
return X;
}