多项式拟合一般方程法详细推导

多项式拟合之一般方程法

文章目录

• 多项式拟合之一般方程法
• • 1 什么是多项式拟合？
  • 2 怎样拟合？
  • 3 求 b b b
• matlab测试函数

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

1 什么是多项式拟合？

y = a 0 + a 1 x + a 2 x 2 + ⋯ + a n x n = ∑ i = 0 n a i x i \begin{aligned} y&=a_0+a_1x+a_2x^2+\cdots +a_nx^n\\ &=\sum_{i=0}^{n}{a_ix^i} \end{aligned} y​=a0​+a1​x+a2​x2+⋯+an​xn=i=0∑n​ai​xi​

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

2 怎样拟合？

考虑输入样本数据如下：
x = [ x 1 , x 2 , . . . , x m ] T y = [ y 1 , y 2 , . . . , y m ] T x=[x_1,x_2,...,x_m]^T\\ y=[y_1,y_2,...,y_m]^T x=[x1​,x2​,...,xm​]Ty=[y1​,y2​,...,ym​]T
构建方程组：
y 1 = a 0 + a 1 x 1 + a 2 x 1 2 + ⋯ + a n x 1 n y 2 = a 0 + a 1 x 2 + a 2 x 2 2 + ⋯ + a n x 2 n ⋮ y m = a 0 + a 1 x m + a 2 x m 2 + ⋯ + a n x m n \begin{aligned} y_1&=a_0+a_1x_1+a_2x_1^2+\cdots +a_nx_1^n\\ y_2&=a_0+a_1x_2+a_2x_2^2+\cdots +a_nx_2^n\\ \vdots\\ y_m&=a_0+a_1x_m+a_2x_m^2+\cdots +a_nx_m^n\\ \end{aligned} y1​y2​⋮ym​​=a0​+a1​x1​+a2​x12​+⋯+an​x1n​=a0​+a1​x2​+a2​x22​+⋯+an​x2n​=a0​+a1​xm​+a2​xm2​+⋯+an​xmn​​

令：

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

3 求 b b b

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

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

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

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

  正交矩阵

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

  QR分解

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

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

  施密特正交化

  这个请大家自行学习。

  完成QR分解后
  A b = f Q R b = f R b = Q T f b = R − 1 Q t f \begin{aligned} Ab &= f\\ QRb&=f\\ Rb&=Q^Tf\\ b&=R^{-1}Q^tf \end{aligned} AbQRbRbb​=f=f=QTf=R−1Qtf​

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;