Levenberg-Marquardt算法实现高斯曲线拟合（qt creator）

基于qt creator开发环境下的高斯曲线拟合实现过程：

空气VOCs色谱图得到的一系列离散数据，色谱峰处符号高斯分布，故采用高斯函数对其进行曲线拟合。开发环境为qt creator，拟合算法选用Levenberg-Marquardt，结果与origin拟合结果一致。Matlab中具有强大的矩阵运算功能，容易实现。本文软件在qt creator下实现则调用了Eigen矩阵库，进行矩阵操作。

1、qt creator下加载Eigen库过程

Eigen官网http://eigen.tuxfamily.org/index.php?title=Main_Page下载Eignen库，解压并命名eigen3，将其放在工程文件下，在头文件中包含即可。

#include <eigen3/Eigen/Dense>

#include <eigen3/Eigen/Sparse>

2、Levenberg-Marquardt算法

参考链接：https://blog.csdn.net/stihy/article/details/52737723

.h代码如下：

#ifndef MAINWINDOW_H

#define MAINWINDOW_H

#include <QMainWindow>

#include <eigen3/Eigen/Dense>

#include <eigen3/Eigen/Sparse>

#include <math.h>

#include <iomanip>

#include <QDebug>

using namespace Eigen;

namespace Ui {

class MainWindow;

class MainWindow : public QMainWindow

    Q_OBJECT

public:

    explicit MainWindow(QWidget *parent = 0);

    double maxMatrixDiagonale(const MatrixXd &Hessian);

    double linerDeltaL(const VectorXd& step, const VectorXd& gradient, const double u);

    void levenMar(const VectorXd &input, const VectorXd &output, VectorXd &params);

    double Func(const VectorXd& obj);

    VectorXd objF(const VectorXd& input, const VectorXd& output, const VectorXd& params);

    double func(const VectorXd& input, const VectorXd& output, const VectorXd& params, double objIndex);

    MatrixXd Jacobin(const VectorXd& input, const VectorXd& output, const VectorXd& params);

    double Deriv(const VectorXd& input, const VectorXd& output, int objIndex, const VectorXd& params,

                     int paraIndex);

    ~MainWindow();

private slots:

    void on_pushButton_matrix_clicked();

private:

    Ui::MainWindow *ui;

};

#endif // MAINWINDOW_H

cpp代码如下：

#include "mainwindow.h"

#include "ui_mainwindow.h"

#include<iostream>

using namespace std;

const double DERIV_STEP = 1e-5;

const int MAX_ITER = 400;

MainWindow::MainWindow(QWidget *parent) :

    QMainWindow(parent),

    ui(new Ui::MainWindow)

    ui->setupUi(this);

MainWindow::~MainWindow()

    delete ui;

void MainWindow::on_pushButton_matrix_clicked()  //矩阵计算

    /***********函数模型y=x1+x2*exp(-0.5 *((x-x3)/x4)*((x-x3)/x4))***********/

        VectorXd  params_levenMar(4);            //高斯函数

        params_levenMar(0) = 494.5472;

        params_levenMar(1) = 44.3528;

        params_levenMar(2) = 17;

        params_levenMar(3) = 2;

        double m[27] = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27};

        double y[27] = {494.5472,494.6084,494.8052,495.2096,495.872,497.08,499.3508,502.1232,505.9028,510.85,516.218,521.6704,527.1972,532.0016,535.6004,537.9836,538.9,538.7624,537.606,535.7508,533.7564,531.7208,529.8448,528.3388,527.2388,526.498,526.2584 };

        //generate random data using these parameter

        int total_data = 27;

        VectorXd input(total_data);

        VectorXd output(total_data);

        //load observation data

        for (int i = 0; i < total_data; i++)

              input(i) = m[i];

              output(i) = y[i];

           levenMar(input, output, params_levenMar);

           cout << "Levenberg-Marquardt parameter: " << endl << params_levenMar << endl << endl << endl;

double MainWindow::maxMatrixDiagonale(const MatrixXd &Hessian)

    int max = 0;

    for(int i = 0; i < Hessian.rows(); i++)

        if(Hessian(i,i) > max)

            max = Hessian(i,i);

    return max;

double MainWindow::linerDeltaL(const VectorXd &step, const VectorXd &gradient, const double u)

    double L = step.transpose() * (u * step - gradient);

    return L/2;

void MainWindow::levenMar(const VectorXd &input, const VectorXd &output, VectorXd &params)

   // int errNum = input.rows();      //error num

    int paraNum = params.rows();    //parameter num

    //initial parameter

    VectorXf obj = objF(input,output,params);

    MatrixXf Jac = Jacobin(input, output, params);  //jacobin

    MatrixXf A = Jac.transpose() * Jac;             //Hessian

    VectorXd gradient = Jac.transpose() * obj;      //gradient

    //initial parameter tao v epsilon1 epsilon2

    double tao = 1e-3;

    long long v = 2;

    double eps1 = 1e-12, eps2 = 1e-12;

    double u = tao * maxMatrixDiagonale(A);

    bool found = gradient.norm() <= eps1;

    if(found) return;

   // double last_sum = 0;

    int iterCnt = 0;

   while (iterCnt < MAX_ITER)

        VectorXd obj = objF(input,output,params);

        MatrixXd Jac = Jacobin(input, output, params);  //jacobin

        MatrixXd A = Jac.transpose() * Jac;             //Hessian

        VectorXd gradient = Jac.transpose() * obj;      //gradient

        if( gradient.norm() <= eps1 )

            cout << "stop g(x) = 0 for a local minimizer optimizer." << endl;

            break;

        cout << "A: " << endl << A << endl;

        VectorXd step = (A + u * MatrixXd::Identity(paraNum, paraNum)).inverse() * gradient; //negtive Hlm.

        cout << "step: " << endl << step << endl;

        if( step.norm() <= eps2*(params.norm() + eps2) )

            cout << "stop because change in x is small" << endl;

            break;

        VectorXd paramsNew(params.rows());

        paramsNew = params - step; //h_lm = -step;

        //compute f(x)

        obj = objF(input,output,params);

        //compute f(x_new)

        VectorXd obj_new = objF(input,output,paramsNew);

        double deltaF = Func(obj) - Func(obj_new);

        double deltaL = linerDeltaL(-1 * step, gradient, u);

        double roi = deltaF / deltaL;

        cout << "roi is : " << roi << endl;

        if(roi > 0)

            params = paramsNew;

            u *= max(1.0/3.0, 1-pow(2*roi-1, 3));

            v = 2;

        else

            u = u * v;

            v = v * 2;

        cout << "u = " << u << " v = " << v << endl;

        iterCnt++;

        cout << "Iterator " << iterCnt << " times, result is :" << endl << endl;

double MainWindow::Func(const VectorXd &obj)

      return obj.squaredNorm()/2;

//return vector make up of func() element.

VectorXd  MainWindow::objF(const VectorXd& input, const VectorXd& output, const VectorXd& params)

    VectorXd obj(input.rows());

    for(int i = 0; i < input.rows(); i++)

        obj(i) = func(input, output, params, i);

    return obj;

double  MainWindow::func(const VectorXd& input, const VectorXd& output, const VectorXd& params, double objIndex)

    double x1 = params(0);

    double x2 = params(1);

    double x3 = params(2);

    double x4 = params(3);

    double t = input(objIndex);

    double f = output(objIndex);

   // return x1*exp(-x2 * t)+x3 - f;       //返回误差

    return x2*exp(-0.5 *((t-x3)/x4)*((t-x3)/x4))+x1 -f;

   // return x1 * sin(x2 * t) + x3 * cos( x4 * t) - f;

/***************求雅克比矩阵*********************/

MatrixXd  MainWindow::Jacobin(const VectorXd& input, const VectorXd& output, const VectorXd& params)

    int rowNum = input.rows();

    int colNum = params.rows();

    MatrixXd Jac(rowNum, colNum);

    for (int i = 0; i < rowNum; i++)

        for (int j = 0; j < colNum; j++)

            Jac(i,j) = Deriv(input, output, i, params, j);

    return Jac;

double MainWindow::Deriv(const VectorXd& input, const VectorXd& output, int objIndex, const VectorXd& params,

                 int paraIndex)

     VectorXd para1 = params;

     VectorXd para2 = params;

     para1(paraIndex) -= DERIV_STEP;

     para2(paraIndex) += DERIV_STEP;

     double obj1 = func(input, output, para1, objIndex);

     double obj2 = func(input, output, para2, objIndex);

     return (obj2 - obj1) / (2 * DERIV_STEP);

Levenberg-Marquardt算法实现高斯曲线拟合（qt creator）

猜你喜欢