基于最小二乘法的多项式曲线拟合（附C++代码）

智驾人在路上

于 2024-08-31 18:08:39 发布

阅读量1.5k

点赞数 30

文章标签：最小二乘法算法机器学习自动驾驶

本文链接：https://blog.csdn.net/qq_40875526/article/details/141753641

版权

1.最小二乘法的一般原理

将测量或实验得到的数据点（ $x_{i}$ ， $y_{i}$ ），拟合成最佳曲线 $\varphi$ （x），即y = $\varphi$ （x）-> y = f（x）

，数据特点如下：

数据量较大
数据是通过测量得到，数据本身有一定的误差。

分析：

若用插值法，通过这几个已知点所求得的插值多项式必定是高次插值多项式，而高次插值是数值不稳定的。
由于数据本身存在的误差，利用插值所得到的插值多项式必保留了所有测量误差，所得结果与实际问题误差较大。
对于数据拟合问题，一般来讲，并不要求所得到的近似解表达式通过所有已知点，而只要求尽可能通过它们附近，这样可抵消原数据组中测量误差。

问题：

用 $\varphi(x)$ 拟合数据（ $x_{i}$ ， $y_{i}$ ），（i = 1,2,...,n），标准应是什么？

（1） $\sum_{i=1}^{n}[\varphi (x_{i}) - y_{i}]$ 偏差之和最小，但偏差有正有负，在求和时可能互相抵消

（2） $\sum_{i=1}^{n}|\varphi (x_{i}) - y_{i}|$ 偏差的绝对值之和最小，但有绝对值符号，不便于分析

（3） $\sum_{i=1}^{n}[\varphi (x_{i}) - y_{i}]^{2}$ 偏差平方和最小，便于分析。

所以，根据偏差的平方和为最小的条件来选择 $\varphi(x)$ 的方法称为最小二乘法（最小二乘曲线拟合法）。

2.拟合函数

用 $\varphi(x)$ 去拟合数据（ $x_{i}$ ， $y_{i}$ ）（i = 1,2,...,n）， $\varphi(x)$ 一般都是由m个线性无关函数 $\varphi_{1}(x)$ ， $\varphi_{2}(x)$ ，...， $\varphi_{m}(x)$ （基函数）线性组合而成，即:

$\varphi(x)$ = $\varphi (x) = a_{1}\varphi _{1}(x) + a_{2}\varphi _{2}(x) + ... + a_{m}\varphi _{m}(x) (m < n-1)$

常见的 $\varphi_{1}(x)$ ， $\varphi_{2}(x)$ ，...， $\varphi_{m}(x)$ 如下：

1, x, $x^{2}$ ,..., $x^{m}$

sinx, sin2x, sin3x, ... , sinmx

$e^{\lambda _{1}x}$ , $e^{\lambda _{2}x}$ , $e^{\lambda _{3}x}$ ,..., $e^{\lambda _{m}x}$

3.用最小二乘法求解矛盾方程组

求解线性方程组时，通常要求未知数的个数与方程式的个数相等。若方程式的个数多于未知数的个数，方程无解，称为矛盾（超定）方程组。

（1）矛盾方程组求解过程

最小二乘法的基本思想：矛盾方程组无精确解，只能寻求某种意义下的近似解，这种近似解并非对精确解之近似，而是寻求各未知数的一组值，使方程中各式能近似相等。

令 $R_{i} = \sum_{j=1}^{m}a_{ij}x_{j} - b_{i} (i = 1,2,...,n)$

按最小二乘法原则，求各个方程式误差平方和

$Q = \sum_{i= 1}^{n}R_{i}^2 = \sum_{i=1}^{n}[ \sum_{j=1}^{m}a_{ij}x_{j} - b_{i} ]^2$

若 $x_{j}(j = 1,2,...,m)$ 取值，使误差平方和Q达到最小，则称这组值是矛盾方程组的最优近似解。

Q可看做m个x自变量的二次函数，且连续，故存在一组数 $x_{1},x_{2},...,x_{m},$ 使Q达到最小（极值问题），即

$\frac{\partial Q}{\partial x_{k}} = 0 (k = 1,2,...,m)$

$\frac{\partial Q}{\partial x_{k}} = \sum_{i=1}^{n}2[\sum_{j=1}^{m}a_{ij}x_{j} - b_{i}]a_{ik}$

$\frac{\partial Q}{\partial x_{k}} = 2\sum_{i=1}^{n}[\sum_{j=1}^{m}a_{ij}a_{ik}x_{j} - a_{ik}b_{i}]$

$\frac{\partial Q}{\partial x_{k}} = 2\sum_{i=1}^{n}(\sum_{j=1}^{m}a_{ij}a_{ik})x_{j} - 2\sum_{i=1}^{n}a_{ik}b_{i} = 0$

$\sum_{j=1}^{m}(\sum_{i=1}^{n}a_{ij}a_{ik})x_{j} = \sum_{i=1}^{n}a_{ik}b_{i} (k = 1,2,...,m)$ （2）

这是m个未知量，m个方程的线性方程组，对应于矛盾方程组（1）的正规方程组。

显然，（2）的解是（1）的最优近似解。记

所以，最小二乘法解矛盾方程组的步骤为：

1.计算 $A^{T}A$ 和 $A^{T}b$ ，得正规方程组

$A^{T}Ax = A^{T}b$

2.求解正规方程组，得矛盾方程组的最优近似解。

4.用多项式作最小二乘曲线拟合

取基函数为： $\varphi_{0}(x) = 1$ ， $\varphi_{1}(x) = x$ ， $\varphi_{2}(x) = x^{2}$ ,..., $\varphi_{m}(x) = x^{m}$ ，

则拟合多项式为：

$P(x) = a_{0} + a_{1}x + a_{2}x^{2} + ... + a_{m}x^{m} (m<n-1)$

分析：

由最小二乘的定义，即要通过给定的数据（ $x_{i}$ ， $y_{i}$ ）（i = 1,2,...,n）确定系数aj，使得在各个点上的误差的平方和为最小。

将n个数据点代入多项式P(x)，就得到一个具有m+1个未知数aj的n个方程的矛盾方程组：

记 $a = (a_{0},a_{1},...,a{n})^T, y = (y_{0},y_{1},...,y{n})^T$

所以， Aa = y

它对应的正规方程组为：

$A^{T}Aa = A^{T}y$

只需要计算：

所以，最小二乘法的拟合步骤为：

1.计算正规方程组的系数矩阵和常数项各元素

2.利用迭代法求正规方程组的解 $a_{0},a_{1},...,a_{m}$

$P(x) = a_{0} + a_{1}x + a_{2}x^{2} + ... + a_{m}x^{m}$

5.用高斯消元法求正规方程组的解

高斯消元法的具体示例：

假设我们有这样一个方程：

可以得到对应的增广矩阵为：

我们将使用列主元的方法

（1）找主元。主元就是左边矩阵的第一列绝对值最大的数。

（2）如果需要，进行行交换，这样保证主元在第一行。

如上图所示的增广矩阵，第一列的主元为 1，因此我们需要进行行交换。交换后的矩阵如下图。

（3）对第一列进行高斯消去，这样我们可以得到下图的增广矩阵。

（4）查找下一列（第二列）的主元。注意已经完成的列不需要参与。

（5）交换行，如果需要。我们必须保证这个主元在这列的对角线位置。

（6）对第二列进行高斯消去，这样我们可以得到下图的增广矩阵。

（7）查找下一列（第三列）的主元。注意已经完成的列不需要参与，交换行，如果需要。

（8）对第三列进行高斯消去，这样我们可以得到下图的增广矩阵。

（9）到这里，我们就形成了上三角矩阵，下面我们用向后替换算法求解。

6.C++代码实现

#include <iostream>
#include <cmath>
#include <cstdint>

#define Point_Num_Max 20

class Fit
{
	float ssr;
	float sse;
	float rmse;
	float fitedYs[Point_Num_Max];

public:
	Fit() : ssr(0), sse(0), rmse(0) {}
	~Fit() {}

	template <typename T>
	void polyfit(const T* x, const T* y, T* factor, uint8_t arrayLength, uint8_t poly_n, bool isSaveFitYs = true)
	{
		if (arrayLength > Point_Num_Max)
		{
			return;
		}
		polyfit1(x, y, factor, arrayLength, poly_n, isSaveFitYs);
	}

	template <typename T>
	void polyfit1(const T* x, const T* y, T* factor, size_t length, uint8_t poly_n, bool isSaveFitYs = true)
	{
		if (length > Point_Num_Max)
		{
			//std::cerr << "Error: Input length exceeds maximum allowed size." << std::endl;
			return;
		}

		uint8_t        i, j;
		float tempx[Point_Num_Max] = { 0 };
		float tempy[Point_Num_Max] = { 0 };
		float sumxx[8] = { 0 }; // Adjust size to 2*poly_n+1
		float ata[16] = { 0 }; // Adjust size to (poly_n+1)^2
		float sumxy[8] = { 0 }; // Adjust size to poly_n+1

		// Calculate sums for matrix and vector
		for (i = 0; i < length; i++)
		{
			tempx[i] = 1.0f;
			tempy[i] = y[i];
		}

		for (i = 0; i < 2 * poly_n + 1; i++) //sumxx[0] = length  summxx[1] = sum(xi)   summxx[2] = sum(xi^2)
		{                                    // sumxx[2n] = sum(xi^2n)
			for (sumxx[i] = 0, j = 0; j < length; j++)
			{
				sumxx[i] += tempx[j];
				tempx[j] *= x[j];
			}
		}

		for (i = 0; i < poly_n + 1; i++) //sumxy[0] = sum(yi)   sumxy[1] = sum(xi*yi// sumxy[n] = sum(xi^n*yi)
		{
			for (sumxy[i] = 0, j = 0; j < length; j++)
			{
				sumxy[i] += tempy[j];
				tempy[j] *= x[j];
			}
		}

		for (i = 0; i < poly_n + 1; i++)
			for (j = 0; j < poly_n + 1; j++)
				ata[i * (poly_n + 1) + j] = sumxx[i + j]; //构建系数矩阵

		gauss_solve(poly_n + 1, ata, factor, sumxy);
		calcError(x, y, factor, length, this->ssr, this->sse, this->rmse, isSaveFitYs);
	}

	template <typename T>
	float getY(const T x, T* factor, uint8_t poly_n) const
	{
		float ans(0);
		for (size_t i = 0; i <= poly_n; ++i)
		{
			ans += factor[i] * pow((float)x, (uint8_t)i);
		}
		return ans;
	}

	template <typename T>
	static T Mean(const T* v, size_t length)
	{
		T total(0);
		for (size_t i = 0; i < length; ++i)
		{
			total += v[i];
		}
		return (total / length);
	}

private:
	template <typename T>
	void calcError(
		const T*        x,
		const T*        y,
		T*              factor,
		size_t          length,
		float& r_ssr,
		float& r_sse,
		float& r_rmse,
		bool            isSaveFitYs = true)
	{
		T mean_y = Mean<T>(y, length);
		T yi(0);
		r_ssr = 0;
		r_sse = 0;

		for (uint8_t i = 0; i < length; ++i)
		{
			yi = getY(x[i], factor, 3); // Assume poly_n is 3 for cubic
			r_ssr += ((yi - mean_y) * (yi - mean_y));
			r_sse += ((yi - y[i]) * (yi - y[i]));
			if (isSaveFitYs)
			{
				fitedYs[i] = yi;
			}
		}
		r_rmse = sqrt(r_sse / float(length));
	}

	template <typename T>
	void gauss_solve(uint8_t n, T* A, T* x, T* b)
	{
		gauss_solve1(n, A, x, b);
	}

	template <typename T>
	void gauss_solve1(uint8_t n, T* A, T* x, T* b) //如果是三次的话，则有四个系数，n = 4
	{
		uint8_t        i, j, k, r;
		float max;
		for (k = 0; k < n - 1; k++) //遍历A矩阵的每一列
		{
			max = abs(A[k * n + k]); /*find maximum*/
			r = k;
			for (i = k + 1; i < n; i++) //计算第一列的最大值是在第几行
			{
				if (max < abs(A[i * n + k])) // Compare with A[i * n + k]
				{
					max = abs(A[i * n + k]);
					r = i;
				}
			}
			if (r != k) //交换增广矩阵的行，将第一列最大值的行放到第一行
			{
				for (i = 0; i < n; i++) /*change array:A[k]&A[r] */
				{
					max = A[k * n + i];
					A[k * n + i] = A[r * n + i];
					A[r * n + i] = max;
				}
				max = b[k]; /*change array:b[k]&b[r] */
				b[k] = b[r];
				b[r] = max;
			} //将第一列消0
			for (i = k + 1; i < n; i++)
			{
				for (j = k + 1; j < n; j++)
				{
					A[i * n + j] -= A[i * n + k] * A[k * n + j] / A[k * n + k];
				}
				b[i] -= A[i * n + k] * b[k] / A[k * n + k];
			}
		}

		//采用后向替换算法求解
		for (i = n - 1; i >= 0 && i < n; i--)
		{
			x[i] = b[i];
			for (j = i + 1; j < n; j++)
			{
				x[i] -= A[i * n + j] * x[j];
			}
			x[i] /= A[i * n + i];
		}
	}
};

int main()
{
	float m_LongDist_f[20] = { -2.39526,-4.5901,-6.78494,-8.97879,-11.1726,-13.366,-15.5595,-17.7523,-19.9452,-22.1374,-24.3296,
		-26.5202,-28.7108,-30.9001,-33.0894,-35.2786,-37.4679,-39.6572,-41.8465,-44.0357 };
	float m_LatDist_f[20] = { -0.0449439,-0.047939,-0.0509341,-0.0550729,-0.0592116,-0.0665418,-0.073872,-0.0798808,-0.0858895,-0.0882831,
		-0.0906766,-0.0910919,-0.0915072,-0.101326,-0.111145,-0.120964,-0.130783,-0.140602,-0.150421,-0.16024 };

	Fit fit;
	float factors[4] = { 0 }; // To store polynomial coefficients: [a0, a1, a2, a3]

	fit.polyfit(m_LongDist_f, m_LatDist_f, factors, 20, 3,false); // 3 for cubic polynomial

	std::cout << "Polynomial coefficients are:" << std::endl;
	for (int i = 0; i < 4; ++i)
	{
		std::cout << "a" << i << " = " << factors[i] << std::endl;
	}

	return 0;
}

运算结果为：