曲线拟合的最小二乘原理

文章目录

1 什么是最小二乘
2 最小二乘原理
3 最小二乘应用示例
4 法方程到底是什么

1 什么是最小二乘

在科学实验的统计方法研究中，往往要从一组实验数据 ( x i , y i ) ( i = 0 , 1 , 2 , … , m ) (x_i,y_i)(i=0,1,2,…,m) (xi,yi)(i=0,1,2,…,m) 中寻找自变量 x x x 与因变量 y y y 之间的函数关系 y = F ( x ) y=F(x) y=F(x). 由于观测数据往往不准确，因此不要求 y = F ( x ) y=F(x) y=F(x) 经过所有点 ( x i , y i ) (x_i,y_i) (xi,yi)，而只要求在给定点 x i x_i xi 上误差 δ i = F ( x i ) − y i ( i = 0 , 1 , 2 , … , m ) δ_i=F(x_i )-y_i (i=0,1,2,…,m) δi=F(xi)−yi(i=0,1,2,…,m)按某种标准最小。若记 δ = ( δ 0 , δ 1 , … , δ m ) T δ=(δ_0,δ_1,…,δ_m)^T δ=(δ0,δ1,…,δm)T，就是要求向量 δ δ δ 的范数最小，通常采用计算较为简单的欧式范数 ‖ δ ‖ 2 ‖δ‖_2 ‖δ‖2 作为误差衡量的标准。

关于最小二乘的一般提法是：对给定的一组数据 ( x i , y i ) ( i = 0 , 1 , 2 , … , m ) (x_i,y_i)(i=0,1,2,…,m) (xi,yi)(i=0,1,2,…,m)，要求在函数类 φ { φ 0 , φ 1 , … , φ n } φ\{φ_0,φ_1,…,φ_n\} φ{φ0,φ1,…,φn}中找到一个函数 y = S ∗ ( x ) y=S^* (x) y=S∗(x)，使误差平方和最小，即
∥ δ ∥ 2 2 = ∑ i = 0 m δ i 2 = ∑ i = 0 m [ S ∗ ( x i ) − y i ] 2 = m i n S ( x ) ∈ φ ⁡ ∑ i = 0 m [ S ( x i ) − y i ] 2 − − − − ( F o r m u l a . 1 ) \Vert δ \Vert_2^2=∑_{i=0}^m{δ_i^2}=∑_{i=0}^m[S^* (x_i )-y_i ]^2=min_{S(x)∈φ}⁡∑_{i=0}^m[S(x_i )-y_i ]^2----(Formula.1) ∥δ∥22=i=0∑mδi2=i=0∑m[S∗(xi)−yi]2=minS(x)∈φ⁡i=0∑m[S(xi)−yi]2−−−−(Formula.1)

其中，
S ( x ) = a 0 φ 0 ( x ) + a 1 φ 1 ( x ) + ⋯ + a 1 φ 1 ( x ) , ( n < m ) − − − − ( F o r m u l a . 2 ) S(x)=a_0 φ_0 (x)+a_1 φ_1 (x)+⋯+a_1 φ_1 (x) ,(n<m)----(Formula.2) S(x)=a0φ0(x)+a1φ1(x)+⋯+a1φ1(x),(n<m)−−−−(Formula.2)

这就是一般的最小二乘逼近，用几何语言说，就称为曲线拟合的最小二乘法。

2 最小二乘原理

用最小二乘法求拟合曲线时，首先要确定 S ( x ) S(x) S(x) 的形式。这不单纯是数学问题，还与所研究问题的运动规律及所得观测数据 ( x i , y i ) (x_i,y_i) (xi,yi) 有关；通常要从问题的运动规律及给定数据描图，确定 S ( x ) S(x) S(x) 的形式，并通过实际计算选出较好的结果。 S ( x ) S(x) S(x)的一般表达式为 F o r m u l a . 2 Formula.2 Formula.2 式表示的线性形式。为了使问题的提法更有一般性，通常把最小二乘法中 ∥ δ ∥ 2 2 \Vert δ \Vert_2^2 ∥δ∥22都考虑为加权平方和
∥ δ ∥ 2 2 = ∑ i = 0 m ω ( x i ) [ S ( x i ) − f ( x i ) ] 2 − − − − ( F o r m u l a . 3 ) \Vert δ \Vert_2^2=\sum_{i=0}^m\omega(x_i)\left[S(x_i)-f(x_i)\right]^2----(Formula.3) ∥δ∥22=i=0∑mω(xi)[S(xi)−f(xi)]2−−−−(Formula.3)
这里， ω ( x ) ≥ 0 ω(x)≥0 ω(x)≥0是 [ a , b ] [a ,b] [a,b]上的权函数，它表示不同点 ( x i , f ( x i ) ) (x_i,f(x_i)) (xi,f(xi))处的数据比重不同，通常表示在点 ( x i , f ( x i ) ) (x_i,f(x_i)) (xi,f(xi)) 处重复观测的次数。用最小二乘法求拟合曲线的问题，就是在形如 ( F o r m u l a . 2 ) (Formula.2) (Formula.2) 式的 S ( x ) S(x) S(x) 中求一函数 y = S ∗ ( x ) y=S^* (x) y=S∗(x)，使 F o r m u l a . 3 Formula.3 Formula.3 式取得最小。可以将这个问题转化为求多元函数
I ( a 0 , a 1 , … , a n ) = ∑ i = 0 m ω ( x i ) [ ∑ j = 0 n a j φ j ( x i ) − f ( x i ) ] 2 − − − − − ( F o r m u l a . 4 ) I(a_0,a_1,…,a_n )=∑_{i=0}^mω(x_i)\left[∑_{j=0}^na_j φ_j (x_i )-f(x_i)\right]^2-----(Formula.4) I(a0,a1,…,an)=i=0∑mω(xi)[j=0∑najφj(xi)−f(xi)]2−−−−−(Formula.4)
的极小值点 ( a 0 ∗ , a 1 ∗ , … , a n ∗ ) (a_0^*,a_1^*,…,a_n^*) (a0∗,a1∗,…,an∗) 的问题。

由求多元函数极值的必要条件，有
∂ I ∂ a k = 2 ∑ i = 0 m ω ( x i ) [ ∑ j = 0 n a j φ j ( x i ) − f ( x i ) ] φ k ( x i ) = 0 , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 5 ) \cfrac{∂I}{∂a_k}=2∑_{i=0}^mω(x_i )\left[∑_{j=0}^na_jφ_j(x_i)-f(x_i)\right] φ_k (x_i )=0,(k=0,1,…,n)----(Formula.5) ∂ak∂I=2i=0∑mω(xi)[j=0∑najφj(xi)−f(xi)]φk(xi)=0,(k=0,1,…,n)−−−−(Formula.5)
若记
( φ j , φ k ) = ∑ i = 0 m ω ( x i ) φ j ( x i ) φ k ( x i ) − − − − ( F o r m u l a . 6 ) (φ_j,φ_k )= ∑_{i=0}^mω(x_i ) φ_j (x_i ) φ_k (x_i )----(Formula.6) (φj,φk)=i=0∑mω(xi)φj(xi)φk(xi)−−−−(Formula.6)
( f , φ k ) = ∑ i = 0 m ω ( x i ) f ( x i ) φ k ( x i ) ≡ d k , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 7 ) (f,φ_k )=∑_{i=0}^mω(x_i ) f(x_i ) φ_k (x_i )≡d_k,(k=0,1,…,n)----(Formula.7) (f,φk)=i=0∑mω(xi)f(xi)φk(xi)≡dk,(k=0,1,…,n)−−−−(Formula.7)
则 F o r m u l a . 4 Formula.4 Formula.4 可改写为
∑ j = 0 n ( φ k , φ j ) a j = d k , ( k = 0 , 1 , … , n ) − − − − ( F o r m u l a . 8 ) ∑_{j=0}^n(φ_k,φ_j)a_j=d_k,(k=0,1,…,n)----(Formula.8) j=0∑n(φk,φj)aj=dk,(k=0,1,…,n)−−−−(Formula.8)
F o r m u l a . 8 Formula.8 Formula.8 式称为法方程，矩阵形式为
G a = d − − − − ( F o r m u l a . 9 ) Ga=d----(Formula.9) Ga=d−−−−(Formula.9)
其中， a = ( a 0 , a 1 , . . . , a n ) T a=(a_0,a_1,...,a_n)^T a=(a0,a1,...,an)T， d = ( d 0 , d 1 , . . . , d n ) T d=(d_0,d_1,...,d_n)^T d=(d0,d1,...,dn)T
G = [ ( φ 0 , φ 0 ) ( φ 0 , φ 1 ) ⋯ ( φ 0 , φ n ) ( φ 1 , φ 0 ) ( φ 1 , φ 1 ) ⋯ ( φ 1 , φ n ) ⋮ ⋮ ⋱ ⋮ ( φ n , φ 0 ) ( φ n , φ 1 ) ⋯ ( φ n , φ n ) ] G=\begin{bmatrix} {(φ_0,φ_0)}&{(φ_0,φ_1)}&{\cdots}&{(φ_0,φ_n)}\\ {(φ_1,φ_0)}&{(φ_1,φ_1)}&{\cdots}&{(φ_1,φ_n)}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {(φ_n,φ_0)}&{(φ_n,φ_1)}&{\cdots}&{(φ_n,φ_n)}\\ \end{bmatrix} G=⎣⎢⎢⎢⎡(φ0,φ0)(φ1,φ0)⋮(φn,φ0)(φ0,φ1)(φ1,φ1)⋮(φn,φ1)⋯⋯⋱⋯(φ0,φn)(φ1,φn)⋮(φn,φn)⎦⎥⎥⎥⎤
由于 φ 0 , φ 1 , . . . , φ n φ_0,φ_1,...,φ_n φ0,φ1,...,φn 线性无关，因此 ∣ G ∣ ≠ 0 |G|≠0 ∣G∣=0，方程组 F o r m u l a . 8 Formula.8 Formula.8 存在唯一解
a k = a k ∗ ，（ k = 0 , 1 , . . . , n ） − − − − ( F o r m u l a . 10 ) a_k=a_k^*，（k=0,1,...,n）----(Formula.10) ak=ak∗，（k=0,1,...,n）−−−−(Formula.10)
从而得到函数 f ( x ) f(x) f(x) 的最小二乘解为
S ∗ ( x ) = a 0 ∗ φ 0 ( x ) + a 1 ∗ φ 1 ( x ) + . . . + a n ∗ φ n ( x ) − − − − ( F o r m u l a . 11 ) S^*(x)=a_0^*φ_0(x)+a_1^*φ_1(x)+...+a_n^*φ_n(x)----(Formula.11) S∗(x)=a0∗φ0(x)+a1∗φ1(x)+...+an∗φn(x)−−−−(Formula.11)
可以证明 S ∗ ( x ) S^*(x) S∗(x) 对于任何形如 F o r m u l a . 2 Formula.2 Formula.2 式的 S ( x ) S(x) S(x) ，都有
∑ i = 0 m ω ( x i ) [ S ∗ ( x i ) − f ( x i ) ] 2 ≤ ∑ i = 0 m ω ( x i ) [ S ( x i ) − f ( x i ) ] 2 − − − − ( F o r m u l a . 12 ) ∑_{i=0}^mω(x_i ) [S^* (x_i )-f(x_i )]^2≤∑_{i=0}^mω(x_i ) [S(x_i)-f(x_i )]^2 ----(Formula.12) i=0∑mω(xi)[S∗(xi)−f(xi)]2≤i=0∑mω(xi)[S(xi)−f(xi)]2−−−−(Formula.12)

也就是说，只要对一组数据的法方程进行求解，就可以得到唯一一组多项式的系数解。如何求解方程组，将会在后续的博客中给出。

3 最小二乘应用示例

下面通过一个例题进一步理解曲线的最小二乘

4 法方程到底是什么

相信不少人对于法方程 G a = d Ga=d Ga=d 中 G G G 的元素到底是什么存在疑问，那么 G G G 中的 φ φ φ 到底是什么呢？下面通过分析一阶、二阶、三阶多项式拟合的法方程，帮助大家理解这个问题。

回顾一下 G G G
G = [ ( φ 0 , φ 0 ) ( φ 0 , φ 1 ) ⋯ ( φ 0 , φ n ) ( φ 1 , φ 0 ) ( φ 1 , φ 1 ) ⋯ ( φ 1 , φ n ) ⋮ ⋮ ⋱ ⋮ ( φ n , φ 0 ) ( φ n , φ 1 ) ⋯ ( φ n , φ n ) ] G=\begin{bmatrix} {(φ_0,φ_0)}&{(φ_0,φ_1)}&{\cdots}&{(φ_0,φ_n)}\\ {(φ_1,φ_0)}&{(φ_1,φ_1)}&{\cdots}&{(φ_1,φ_n)}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {(φ_n,φ_0)}&{(φ_n,φ_1)}&{\cdots}&{(φ_n,φ_n)}\\ \end{bmatrix} G=⎣⎢⎢⎢⎡(φ0,φ0)(φ1,φ0)⋮(φn,φ0)(φ0,φ1)(φ1,φ1)⋮(φn,φ1)⋯⋯⋱⋯(φ0,φn)(φ1,φn)⋮(φn,φn)⎦⎥⎥⎥⎤

需要明确的一点是：若以 x x x 为自变量，则 φ 0 ( x ) = x 0 = 1 ， φ 1 ( x ) = x 1 ， φ 2 ( x ) = x 2 ， … ， φ n ( x ) = x n φ_0 (x)=x^0=1，φ_1 (x)=x^1，φ_2 (x)=x^2，… ，φ_n (x)=x^n φ0(x)=x0=1，φ1(x)=x1，φ2(x)=x2，…，φn(x)=xn.

已知一组数据点 P = { p 1 ( x 1 , y 1 ) , p 2 ( x 2 , y 2 ) , … , p m ( x m , y m ) } P=\{p_1 (x_1,y_1 ),p_2 (x_2,y_2 ),…,p_m (x_m,y_m )\} P={p1(x1,y1),p2(x2,y2),…,pm(xm,ym)}，每个点只观测一次，即 ω ( x ) ≡ 1 \omega(x)≡1 ω(x)≡1. 分别对其进行一阶、二阶、三阶多项式拟合，对应的拟合函数与法方程如下：

拟合阶数	拟合函数	法方程
1	S ( x ) = a 0 + a 1 x S(x)=a_0+a_1 x S(x)=a0+a1x	[ ∑ i = 1 m x i 0 ∑ i = 1 m x i 1 ∑ i = 1 m x i 1 ∑ i = 1 m x i 2 ] \begin{bmatrix}\sum_{i=1}^mx_i^0&\sum_{i=1}^mx_i^1&\\\sum_{i=1}^mx_i^1&\sum_{i=1}^mx_i^2&\\\end{bmatrix} [∑i=1mxi0∑i=1mxi1∑i=1mxi1∑i=1mxi2] [ a 0 a 1 ] \begin{bmatrix}a_0\\a_1\\\end{bmatrix} [a0a1]= [ ∑ i = 1 m x i 0 y i 1 ∑ i = 1 m x i 1 y i 1 ] \begin{bmatrix}\sum_{i=1}^mx_i^0y_i^1\\\sum_{i=1}^mx_i^1y_i^1\\\end{bmatrix} [∑i=1mxi0yi1∑i=1mxi1yi1] (详细)
1	S ( x ) = a 0 + a 1 x S(x)=a_0+a_1 x S(x)=a0+a1x	[ m ∑ x ∑ x ∑ x 2 ] \begin{bmatrix}m&\sum x\\\sum x&\sum x^2\\\end{bmatrix} [m∑x∑x∑x2] [ a 0 a 1 ] \begin{bmatrix}a_0\\a_1\\\end{bmatrix} [a0a1]= [ ∑ y ∑ x y ] \begin{bmatrix}\sum y\\\sum xy\\\end{bmatrix} [∑y∑xy]
2	S ( x ) = a 0 + a 1 x + a 2 x 2 S(x)=a_0+a_1 x+a_2 x^2 S(x)=a0+a1x+a2x2	[ m ∑ x ∑ x 2 ∑ x ∑ x 2 ∑ x 3 ∑ x 2 ∑ x 3 ∑ x 4 ] \begin{bmatrix}m&\sum x&\sum x^2 \\\sum x&\sum x^2&\sum x^3\\\sum x^2&\sum x^3&\sum x^4\\\end{bmatrix} ⎣⎡m∑x∑x2∑x∑x2∑x3∑x2∑x3∑x4⎦⎤ [ a 0 a 1 a 2 ] \begin{bmatrix}a_0\\a_1\\a_2\\\end{bmatrix} ⎣⎡a0a1a2⎦⎤= [ ∑ y ∑ x y ∑ x 2 y ] \begin{bmatrix}\sum y\\\sum xy\\\sum x^2y\\\end{bmatrix} ⎣⎡∑y∑xy∑x2y⎦⎤
3	S ( x ) = a 0 + a 1 x + a 2 x 2 + a 3 x 3 S(x)=a_0+a_1 x+a_2 x^2+a_3 x^3 S(x)=a0+a1x+a2x2+a3x3	[ m ∑ x ∑ x 2 ∑ x 3 ∑ x ∑ x 2 ∑ x 3 ∑ x 4 ∑ x 2 ∑ x 3 ∑ x 4 ∑ x 5 ∑ x 3 ∑ x 4 ∑ x 5 ∑ x 6 ] \begin{bmatrix}m&\sum x&\sum x^2&\sum x^3 \\\sum x&\sum x^2&\sum x^3&\sum x^4\\\sum x^2&\sum x^3&\sum x^4&\sum x^5\\\sum x^3&\sum x^4&\sum x^5&\sum x^6\\\end{bmatrix} ⎣⎢⎢⎡m∑x∑x2∑x3∑x∑x2∑x3∑x4∑x2∑x3∑x4∑x5∑x3∑x4∑x5∑x6⎦⎥⎥⎤ [ a 0 a 1 a 2 a 3 ] \begin{bmatrix}a_0\\a_1\\a_2\\a_3\\\end{bmatrix} ⎣⎢⎢⎡a0a1a2a3⎦⎥⎥⎤= [ ∑ y ∑ x y ∑ x 2 y ∑ x 3 y ] \begin{bmatrix}\sum y\\\sum xy\\\sum x^2y\\\sum x^3y\\\end{bmatrix} ⎣⎢⎢⎡∑y∑xy∑x2y∑x3y⎦⎥⎥⎤
…	…	…

为了便于理解，首先给出一阶多项式详细的法方程。

对比不同阶数的法方程可以看出， G G G 中的每个元素以一定规律对自变量的 k k k 次幂进行求和得到，低阶的法方程是高阶的法方程的一个 “子集”，且矩阵 G G G 为对称正定矩阵。

相关链接

基于Eigen库的线性方程组/矩阵方程求解(方法汇总)

最小二乘曲线拟合的C++实现

参考资料：

《数值分析》李庆扬