多项式拟合之一般方程法

之前经常会用到多项式拟合来拟合函数关系，十分好用，得出的表达式使用起来超级方便。在此对多项式拟合原理进行分析。

1 什么是多项式拟合？

$$\begin{array}{rl}y& =a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n\\ & =\sum _{i=0}^n{a}_i{x}^i\end{array}$$

这是多项式拟合的输出表达式。

2 怎样拟合？

考虑输入样本数据如下：
$$\begin{gathered} x=[x_1,x_2,...,x_m{]}^T \\ y=[y_1,y_2,...,y_m{]}^T \end{gathered}$$
构建方程组：
$$\begin{array}{rl}{y}_1& =a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n\\ y_2& =a_0+a_1{x}_2+a_2{x}_2^2+\cdots +a_n{x}_2^n\\ ⋮& \\ y_m& =a_0+a_1{x}_m+a_2{x}_m^2+\cdots +a_n{x}_m^n\end{array}$$

令：

则上面的方程组可以写为：
$$Ab=f$$
拟合过程就是求$b$的过程

3 求$b$

现在一个很直观的想法，就是：
$$\begin{array}{r}Ab=f\\ b=A^{-1}f\end{array}$$
但是，没这么简单。考虑一下矩阵维度
$$A_{[m,n+1]}{b}_{[n+1,1]}=f_{[m,1]}$$
强行规定输入的$x=[x_1,x_2,...,x_m]$不能重复，则$A$的各个列向量都线性无关

• 当$m=n+1$时，$A$是可逆矩阵，可以直接计算

• 当$m<n+1$时，方程数少于未知数个数，方程不可解，拟合失败

• 当$m>n+1$时，怎么算呢？

  正交矩阵

  对于一个矩阵A，如果有$A{A}^T=E$，则称之为正交矩阵

  QR分解

  将一个矩阵A，分解为一个正交矩阵和一个非奇异上三角阵的过程
  $$A_{[m,n+1]}=Q_{[m,n+1]}{R}_{[n+1,n+1]}$$
  其中$Q$是正交矩阵，$R$是非奇异上三角阵

  那么，如何进行QR分解呢？

  施密特正交化

  这个请大家自行学习。

  完成QR分解后
  $$\begin{array}{rl}Ab& =f\\ QRb& =f\\ Rb& =Q^{T}f\\ b& =R^{-1}{Q}^{t}f\end{array}$$

matlab测试函数

测试数据

%% 准备测试数据
x=1:0.1:50;
y=0.001*x.^4+100*rand(size(x));

plot(x,y,'b')
title('测试数据')

拟合与测试

%% 拟合
n=4;
c = mypolyfit(x',y',n);
%% 测试效果
y1 = mypolyval(c,x');
plot(x,y,'b',x,y1','r')
title('测试拟合效果')

mypolyfit函数

function c = mypolyfit(x,f,n)
A = x.^(0:n);
c = A\f;

mypolyval函数

function y = mypolyval(c,s)
n = length(c)-1;
B = s.^(0:n);
y = B*c;