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

本文介绍基于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;
}

结果:

作为比较,我们在Matlab中采用同样的数据进行计算可得:

可见,结果还是比较相近的。

总结,本文介绍了基于Eigen库的多项式曲线拟合实现(最小二乘法)。

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值