基于Eigen库的多项式曲线拟合实现（最小二乘法）

propor

已于 2024-04-26 14:46:11 修改

阅读量1.1k

点赞数 8

分类专栏：数字信号处理文章标签：最小二乘法算法多项式曲线拟合

于 2024-04-03 14:08:36 首次发布

本文链接：https://blog.csdn.net/propor/article/details/137336234

版权

数字信号处理专栏收录该内容

18 篇文章 8 订阅

订阅专栏

本文介绍基于Eigen库的多项式曲线拟合实现（最小二乘法）。

1.基础知识

1)范德蒙矩阵

范德蒙矩阵是一个n*m的矩阵，定义为

$V= \begin{pmatrix} 1 & x_1 & x_1^2 & \dots & x_1^{m-1}\\ 1 & x_2 & x_2^2 & \dots & x_2^{m-1}\\ 1 & x_3 & x_3^2 & \dots & x_3^{m-1}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & x_n^2 & \dots & x_n^{m-1} \end{pmatrix}$

其第i 行、第j 列可以表示为 $x_{i}^{j-1}$ 。范德蒙矩阵可以应用于多项式的最小二乘法。

2)最小二乘法原理

给出n个点 $(x_{i},y_{i})$ ， $i=1,2,3,\cdot \cdot \cdot ,n$ 求（n-1）次多项式

$p(x)=a_{n-1}x^{n-1}+a_{n-2}x^{n-2}+\cdot \cdot \cdot +a_{1}x^{1}+a_{0}$

满足

$\left\{\begin{matrix} p(x_{1})=a_{n-1}x_{1}^{n-1}+a_{n-2}x_{1}^{n-2}+\cdot \cdot \cdot +a_{1}x_{1}^{1}+a_{0}\\ p(x_{2})=a_{n-1}x_{2}^{n-1}+a_{n-2}x_{2}^{n-2}+\cdot \cdot \cdot +a_{1}x_{2}^{1}+a_{0} \\ \cdot \\ \cdot \\ \cdot \\ p(x_{n})=a_{n-1}x_{n}^{n-1}+a_{n-2}x_{n}^{n-2}+\cdot \cdot \cdot +a_{1}x_{n}^{1}+a_{0} \end{matrix}\right.$

将上面的线性方程组写成矩阵的形式为 $Aa=y$ ，A为n*n阶范德蒙矩阵。若点数大于待求变量个数，需对矩阵做适当变形

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

$a=(A^{T}A)^{-1}A^{T}y$

a就是我们要求的多项式的系数。

2.Eigen库

Eigen库是用C++编写的线性代数库，实现了线性代数以及矩阵分析中所有的计算方法。使用它可以很方便的实现矩阵运算。

Eigen库的下载地址：https://eigen.tuxfamily.org/index.php?title=Main_Page

Eigen库教程地址：https://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html

我这边下载的是“eigen-3.4.0”，使用时将目录下的“Eigen”目录拷贝到工程目录下，在开发环境中添加此目录即可。这里给出一个测试程序测试环境是否正常。

#include <iostream>
#include <Eigen/Dense>
 
int main()
{
    Eigen::MatrixXd m(2, 2);

    m(0, 0) = 1;
    m(1, 0) = 2;
    m(0, 1) = 3;
    m(1, 1) = 4;
    std::cout << "matrix m:\n" << m << std::endl;

    return 0;
}

程序正常编译，输出即可。

3.实现

根据上面的公式，使用Eigen提供的类可轻松实现，参考代码如下：

static bool LeastSquareMethod(Eigen::VectorXf &x, Eigen::VectorXf &y, uint8_t orders, Eigen::VectorXf &coeffs)
{
    Eigen::MatrixXf mtxVandermonde(x.size(), orders + 1);   // Vandermonde matrix of X-axis coordinate vector of sample data
    Eigen::VectorXf vectColVandermonde = x;                 // Vandermonde column
    Eigen::VectorXf vectResult;

    if ((x.size() < 2) || (y.size() < 2) || (x.size() != y.size()) || (orders < 1))
    {
        return false;
    }

    mtxVandermonde.col(0) = Eigen::VectorXf::Constant(x.size(), 1, 1);
    mtxVandermonde.col(1) = vectColVandermonde;

    // construct Vandermonde matrix column by column
    for (int32_t i = 2; i < orders + 1; i++)
    {
        vectColVandermonde = vectColVandermonde.array() * x.array();
        mtxVandermonde.col(i) = vectColVandermonde;
    }

    // calculate coefficients vector of fitted polynomial
    coeffs = (mtxVandermonde.transpose() * mtxVandermonde).inverse() * mtxVandermonde.transpose() * y;

    return true;
}

其中，

x为输入值

y为输出值

coeffs为计算的参数值

4.测试

int main()
{
    Eigen::VectorXf x(5);
    Eigen::VectorXf y(5);

    x(0) = 1;
    x(1) = 2;
    x(2) = 3;
    x(3) = 4;
    x(4) = 5;

    y(0) = 3;
    y(1) = 5;
    y(2) = 8;
    y(3) = 9;
    y(4) = 12;

    Eigen::VectorXf coeffs;
    LeastSquareMethod(x, y, 3, coeffs);

    cout << "The coefficients vector is: \n" << endl;
    for (int32_t i = 0; i < coeffs.size(); i++)
    {
        cout << "coeffs" << i << ": " << coeffs(i) << endl;
    }

    return 0;
}

结果：