基于最小二乘法的多项式曲线拟合：从原理到c++实现

2023-03-25 10:46| 来源: 网络整理| 查看: 265

本文首发于公众号【DeepDriving】，欢迎关注，代码获取方式见文末。

前言

在自动驾驶系统中，通常会用起点、终点和一个三阶多项式来表示一条车道线，多项式系数的求解一般用最小二乘法来实现。

本文首先介绍两种基于最小二乘法的多项式拟合方法的原理，然后基于OpenCV用c++编写了这两种拟合方法的代码，最后通过一个完整的示例来展示如何通过一个离散点集拟合出一条多项式曲线。

基于最小二乘法的多项式拟合原理推导代数方式求解

多项式曲线拟合是指基于一系列的观测点去寻找一个多项式来表示这些点的关系，最小二乘法通过最小化误差的平方和去寻找数据的最佳匹配函数。假设有点集

P={(x1，y1),(x2，y2),…,(xn，yn)}P=\left \{ (x_{1}，y_{1}),(x_{2}，y_{2}), \dots ,(x_{n}，y_{n})\right \} P={(x1，y1),(x2，y2),…,(xn，yn)}

其中，xi∈Rx_{i}\in Rxi∈R和yi∈Ry_{i}\in Ryi∈R的关系满足函数

f(xi)=yi,i=1,2,…,nf(x_{i})=y_{i}, \quad i=1, 2,\dots , nf(xi)=yi,i=1,2,…,n

假设mmm次多项式函数为

f(xi,wj)=w0+w1xi+w2xi2+⋯+wmxim=∑j=0mwjxij\begin{align*} f(x_{i},w_{j}) &= w_{0} + w_{1}x_{i} + w_{2}x_{i}^{2} + \dots +w_{m}x_{i}^{m} \\ &= \sum_{j=0}^{m} w_{j}x_{i}^{j} \end{align*}f(xi,wj)=w0+w1xi+w2xi2+⋯+wmxim=j=0∑mwjxij

其中wjw_{j}wj为多项式系数。如果要用这个mmm次多项式来表示xxx与yyy的关系，那么多项式值与真实值之间的误差为

ei=f(xi)−f(xi,wj)=yi−∑j=0mwjxij\begin{align*} e_{i} &= f(x_{i})-f(x_{i},w_{j}) \\ &= y_{i}-\sum_{j=0}^{m} w_{j}x_{i}^{j} \end{align*}ei=f(xi)−f(xi,wj)=yi−j=0∑mwjxij

采用最小二乘法进行多项式拟合的目的就是寻找一组最佳的多项式系数使得拟合后整个点集的总误差最小，而求总误差最小的问题可以转化为求误差平方和最小。整个点集的误差平方和为

E(w)=∑i=1n(∑j=0m(wjxij−yi))2E(w) = \sum_{i=1}^{n} (\sum_{j=0}^{m} (w_{j}x_{i}^{j} - y_{i}))^2E(w)=i=1∑n(j=0∑m(wjxij−yi))2

要使E(w)E(w)E(w)最小，可以对wkw_{k}wk(k=0,1,2,…,mk=0,1, 2,\dots , mk=0,1,2,…,m)求偏导并令其为零：

∂E(w)wk=2∑i=1n∑j=0m((wjxij−yi)xik)=0\frac{\partial E(w)}{w_{k}}=2\sum_{i=1}^{n} \sum_{j=0}^{m} ((w_{j}x_{i}^{j} - y_{i})x_{i}^{k}) =0wk∂E(w)=2i=1∑nj=0∑m((wjxij−yi)xik)=0

可得

∑i=1n∑j=0mwjxij+k=∑i=1nxikyik=0,1,2,…,m\sum_{i=1}^{n} \sum_{j=0}^{m}w_{j}x_{i}^{j+k} =\sum_{i=1}^{n} x_{i}^{k}y_{i} \\ \quad k=0,1, 2,\dots , mi=1∑nj=0∑mwjxij+k=i=1∑nxikyik=0,1,2,…,m

写成矩阵形式可得

[n∑i=1nxi∑i=1nxi2…∑i=1nxim∑i=1nxi∑i=1nxi2∑i=1nxi3…∑i=1nxim+1∑i=1nxi2∑i=1nxi3∑i=1nxi4…∑i=1nxim+2⋮⋮⋮⋮⋮∑i=1nxim∑i=1nxim+1∑i=1nxim+2…∑i=1nxi2m][w0w1w2⋮wm]=[∑i=1nyi∑i=1nxiyi∑i=1nxi2yi⋮∑i=1nximyi]\begin{bmatrix} n & \sum_{i=1}^{n}x_{i} & \sum_{i=1}^{n}x_{i}^{2} & \dots & \sum_{i=1}^{n}x_{i}^{m}\\ \sum_{i=1}^{n}x_{i} & \sum_{i=1}^{n}x_{i}^{2} & \sum_{i=1}^{n}x_{i}^{3} & \dots & \sum_{i=1}^{n}x_{i}^{m+1}\\ \sum_{i=1}^{n}x_{i}^{2} & \sum_{i=1}^{n}x_{i}^{3} & \sum_{i=1}^{n}x_{i}^{4} & \dots & \sum_{i=1}^{n}x_{i}^{m+2}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ \sum_{i=1}^{n}x_{i}^{m} & \sum_{i=1}^{n}x_{i}^{m+1} & \sum_{i=1}^{n}x_{i}^{m+2} & \dots & \sum_{i=1}^{n}x_{i}^{2m} \end{bmatrix} \begin{bmatrix} w_{0}\\ w_{1}\\ w_{2}\\ \vdots\\ w_{m} \end{bmatrix} = \begin{bmatrix} \sum_{i=1}^{n}y_{i} \\ \sum_{i=1}^{n}x_{i}y_{i} \\ \sum_{i=1}^{n}x_{i}^{2}y_{i}\\ \vdots\\ \sum_{i=1}^{n}x_{i}^{m}y_{i}\\ \end{bmatrix} ⎣⎡n∑i=1nxi∑i=1nxi2⋮∑i=1nxim∑i=1nxi∑i=1nxi2∑i=1nxi3⋮∑i=1nxim+1∑i=1nxi2∑i=1nxi3∑i=1nxi4⋮∑i=1nxim+2………⋮…∑i=1nxim∑i=1nxim+1∑i=1nxim+2⋮∑i=1nxi2m⎦⎤⎣⎡w0w1w2⋮wm⎦⎤=⎣⎡∑i=1nyi∑i=1nxiyi∑i=1nxi2yi⋮∑i=1nximyi⎦⎤

把上式写为AW=BAW=BAW=B，那么只要计算出

∑i=1nxij, j=0,1,2,…,2m\sum_{i=1}^{n}x_{i}^{j}, \ j=0,1,2,\dots, 2mi=1∑nxij, j=0,1,2,…,2m

和

∑i=1nxijyi, j=0,1,2,…,m\sum_{i=1}^{n}x_{i}^{j}y_{i}, \ j=0,1,2,\dots, mi=1∑nxijyi, j=0,1,2,…,m

，再代入上面的式子中构建出矩阵AAA和BBB，那么多项式系数WWW就可以通过求解线性方程组得到。

矩阵方式求解

从上一节的推导过程可以看到，用代数法求解多项式系数的方法非常繁琐而且计算量比较大，我们其实还可以用矩阵求解的方法来简化计算过程。

令

A=[1x1x12…x1m1x2x22…x2m1x3x32…x3m⋮⋮⋮⋮⋮1xnxn2…xnm],W=[w0w1w2⋮wm],B=[y1y2y3⋮yn]A=\begin{bmatrix} 1 & x_{1} & x_{1}^{2} & \dots & x_{1}^{m}\\ 1 & x_{2} & x_{2}^{2} & \dots & x_{2}^{m}\\ 1 & x_{3} & x_{3}^{2} & \dots & x_{3}^{m}\\ \vdots & \vdots & \vdots & \vdots & \vdots\\ 1 & x_{n} & x_{n}^{2} & \dots & x_{n}^{m}\\ \end{bmatrix}, \quad W=\begin{bmatrix} w_{0}\\ w_{1}\\ w_{2}\\ \vdots\\ w_{m} \end{bmatrix}, \quad B= \begin{bmatrix} y_{1} \\ y_{2} \\ y_{3}\\ \vdots\\ y_{n}\\ \end{bmatrix} A=⎣⎡111⋮1x1x2x3⋮xnx12x22x32⋮xn2………⋮…x1mx2mx3m⋮xnm⎦⎤,W=⎣⎡w0w1w2⋮wm⎦⎤,B=⎣⎡y1y2y3⋮yn⎦⎤

那么整个点集的误差平方和可以表示为

E=(AW−B)T(AW−B)=(WTAT−BT)(AW−B)=WTATAW−BTAW−WTATB+BTB\begin{align*} E &= (AW-B)^{T}(AW-B) \\ &= (W^{T}A^{T}-B^{T})(AW-B) \\ &= W^{T}A^{T}AW-B^{T}AW-W^{T}A^{T}B+B^{T}B \end{align*}E=(AW−B)T(AW−B)=(WTAT−BT)(AW−B)=WTATAW−BTAW−WTATB+BTB

对WWW求导可得

∂E∂W=ATAW+WTATA−ATB−ATB+0=2ATAW−2ATB\begin{align*} \frac{\partial E}{\partial W} &= A^{T}AW+W^{T}A^{T}A - A^{T}B - A^{T}B + 0 \\ &= 2A^{T}AW-2A^{T}B \end{align*}∂W∂E=ATAW+WTATA−ATB−ATB+0=2ATAW−2ATB

令导数为零，那么可以得到

ATAW=ATBA^{T}AW = A^{T}BATAW=ATB

解得

W=(ATA)−1ATBW = (A^{T}A)^{-1}A^{T}BW=(ATA)−1ATB

上面求导过程中用到了几个矩阵求导的性质：

\frac{\partial x^{T}ax}{\partial x} = ax+a^{T}x

> >$$ \frac{\partial x^{T}a}{\partial x} = \frac{\partial a^{T}x}{\partial x}=a 如果aaa是对称矩阵，那么

ax+aTx=2axax+a^{T}x=2axax+aTx=2ax

所以有

\frac{\partial W^{T}A^{T}AW}{\partial W} = A^{T}AW+W^{T}A^{T}A=2A^{T}AW

> >$$ \frac{\partial B^{T}AW}{\partial W} =\frac{\partial W^{T}A^{T}B}{\partial W}= A^{T}B 多项式曲线拟合的C++代码实现

学习了前面的理论知识，我们用c++来编码实现一下，毕竟光看一堆数学公式还是不够直观。从前面的理论知识可以知道，用最小二乘法做多项式曲线拟合其实质就是一个矩阵求解的问题，我们可以借助一些常用的数学库比如OpenCV或Eigen来求解。本文将采用OpenCV来实现。

1. 代数方式求解的C++代码实现

如果理解了前面的理论，写代码其实就比较简单了。对照前面理论部分的公式，构建好矩阵A和B，然后调用OpenCV的solve()函数来求解即可：

// 代数方式求解 void PolyFit(const std::vector &points, const int order, cv::Mat *coeff) { const int n = points.size(); cv::Mat A = cv::Mat::zeros(order + 1, order + 1, CV_64FC1); cv::Mat B = cv::Mat::zeros(order + 1, 1, CV_64FC1); // 构建A矩阵 for (int i = 0; i < order + 1; ++i) { for (int j = 0; j < order + 1; ++j) { for (int k = 0; k < n; ++k) { A.at(i, j) += std::pow(points.at(k).x, i + j); } } } // 构建B矩阵 for (int i = 0; i < order + 1; ++i) { for (int k = 0; k < n; ++k) { B.at(i, 0) += std::pow(points.at(k).x, i) * points.at(k).y; } } (*coeff) = cv::Mat::zeros(order + 1, 1, CV_64FC1); // 求解 if (!cv::solve(A, B, *coeff, cv::DECOMP_LU)) { std::cout

【本文地址】

公司简介

联系我们