0 前言
自动驾驶开发中经常涉及到多项式曲线拟合,本文详细描述了使用最小二乘法进行多项式曲线拟合的数学原理,通过样本集构造范德蒙德矩阵,将一元 N 次多项式非线性回归问题转化为 N 元一次线性回归问题,并基于线性代数 C++ 模板库——Eigen 进行了实现,最后,比较了几种实现方法在求解速度与求解精度上的差异。
1 最小二乘法概述
最小二乘法(Least Square Method,LSM)通过最小化误差(也叫残差)的平方和寻找数据的最优函数匹配。
假设给定一组样本数据集 P ( x , y ) P(x, y) P(x,y), P P P 内各数据点 P i ( x i , y i ) ( i = 1 , 2 , 3 , . . . , m ) P_i(x_i, y_i) (i=1, 2, 3, ..., m) Pi(xi,yi)(i=1,2,3,...,m) 来自于对多项式
f ( x i ) = θ 0 + θ 1 x i + θ 2 x i 2 + ⋅ ⋅ ⋅ + θ n x i n f(x_i)=θ_0+θ_1x_i+θ_2x_i^2+···+θ_nx_i^n f(xi)=θ0+θ1xi+θ2xi2+⋅⋅⋅+θnxin
的多次采样,其中:
- m m m 为样本维度
- n n n 为多项式阶数
- θ j ( j = 1 , 2 , 3 , . . . , n ) θ_j (j=1, 2, 3, ..., n) θj(j=1,2,3,...,n) 为多项式的各项系数
针对样本数据集 P P P 内各数据点的误差平方和为:
S = ∑ i = 1 m [ f ( x i ) − y i ] 2 S=\sum_{i=1}^m[f(x_i)-y_i]^2 S=i=1∑m[f(xi)−yi]2
最小二乘法认为,最优函数的各项系数 θ j ( j = 1 , 2 , 3 , . . . , n ) θ_j (j=1, 2, 3, ..., n) θj(j=1,2,3,...,n) 应使得误差平方和 S S S 取得极小值。最小二乘法与极大似然估计有着密切的联系,关于最小二乘法的数学本质可参考文章《如何理解最小二乘法?》。
2 最小二乘法求解多项式曲线系数向量的数学推导
2.1 代数法
由于最优函数的各项系数 θ j ( j = 1 , 2 , 3 , . . . , n ) θ_j (j=1, 2, 3, ..., n) θj(j=1,2,3,...,n) 使得误差平方和 S S S 取得极小值,因而,对于最优函数而言,其误差平方和 S S S 对各多项式系i数 θ j ( j = 1 , 2 , 3 , . . . , n ) θ_j (j=1, 2, 3, ..., n) θj(j=1,2,3,...,n) 的偏导数应满足:
∂ S ∂ θ j = ∑ i = 1 m [ 2 ( θ 0 + θ 1 x i + θ 2 x i 2 + ⋅ ⋅ ⋅ + θ n x i n − y i ) x i j ] = 0 \frac{\partial{S}}{\partial{θ_j}}=\sum_{i=1}^{m}[2(θ_0+θ_1x_i+θ_2x_i^2+···+θ_nx_i^n-y_i)x_i^j]=0 ∂θj∂S=i=1∑m[2(θ0+θ1xi+θ2xi2+⋅⋅⋅+θnxin−yi)xij]=0
整理上式, j j j 分别取 0 , 1 , 2 , . . . , n 0,1,2,...,n 0,1,2,...,n 时,有:
{ m θ 0 + ( ∑ i = 1 m x i ) θ 1 + ( ∑ i = 1 m x i 2 ) θ 2 + ⋯ + ( ∑ i = 1 m x i n ) θ n = ∑ i = 1 m y i ( ∑ i = 1 m x i ) θ 0 + ( ∑ i = 1 m x i 2 ) θ 1 + ( ∑ i = 1 m x i 3 ) θ 2 + ⋯ + ( ∑ i = 1 m x i n + 1 ) θ n = ∑ i = 1 m ( x i y i ) ( ∑ i = 1 m x i 2 ) θ 0 + ( ∑ i = 1 m x i 3 ) θ 1 + ( ∑ i = 1 m x i 4 ) θ 2 + ⋯ + ( ∑ i = 1 m x i n + 2 ) θ n = ∑ i = 1 m ( x i 2 y i ) ⋯ ⋯ ( ∑ i = 1 m x i n ) θ 0 + ( ∑ i = 1 m x i n + 1 ) θ 1 + ( ∑ i = 1 m x i n + 2 ) θ 2 + ⋯ + ( ∑ i = 1 m x i 2 n ) θ n = ∑ i = 1 m ( x i n y i ) \begin{cases} mθ_0+(\sum\limits_{i=1}^{m}x_i)θ_1+(\sum\limits_{i=1}^{m}x_i^2)θ_2+\cdots+(\sum\limits_{i=1}^{m}x_i^n)θ_n=\sum\limits_{i=1}^{m}y_i \\ (\sum\limits_{i=1}^{m}x_i)θ_0+(\sum\limits_{i=1}^{m}x_i^2)θ_1+(\sum\limits_{i=1}^{m}x_i^3)θ_2+\cdots+(\sum\limits_{i=1}^{m}x_i^{n+1})θ_n=\sum\limits_{i=1}^{m}(x_iy_i) \\ (\sum\limits_{i=1}^{m}x_i^2)θ_0+(\sum\limits_{i=1}^{m}x_i^3)θ_1+(\sum\limits_{i=1}^{m}x_i^4)θ_2+\cdots+(\sum\limits_{i=1}^{m}x_i^{n+2})θ_n=\sum\limits_{i=1}^{m}(x_i^2y_i) \\ \cdots\cdots \\ (\sum\limits_{i=1}^{m}x_i^n)θ_0+(\sum\limits_{i=1}^{m}x_i^{n+1})θ_1+(\sum\limits_{i=1}^{m}x_i^{n+2})θ_2+\cdots+(\sum\limits_{i=1}^{m}x_i^{2n})θ_n=\sum\limits_{i=1}^{m}(x_i^ny_i) \end{cases} ⎩