Levenberg–Marquardt算法

算法描述及示例引用自“Machine Learning An Algorithm Perspective 2nd Edition”

使用C++进行浮点数的比较时，一定要注意容差问题。

#include <iostream>#include <math.h>#include <Eigen/Core>#include <Eigen/Dense>using namespace std;double rx(const Eigen::MatrixXd &x,int i){    switch(i)    {        case 0: return 10*(x(1,0) - x(0,0) * x(0,0));        case 1: return 1-x(0,0);     }    return INT_MAX;}double fx(const Eigen::MatrixXd &x){    return 100 *(x(1,0) - x(0,0) * x(0,0)) * (x(1,0) - x(0,0) * x(0,0)) + (1-x(0,0)) * (1-x(0,0));}double Jac(const Eigen::MatrixXd &x,int i,int j){    int index = 2*i+j;    switch(index)    {    case 0: return -20 * x(0,0);    case 1: return 10.0;    case 2: return -1.0;    case 3: return 0.0;    }}void computeRelatedValues(const Eigen::MatrixXd& x,double &fval,Eigen::MatrixXd &r,Eigen::MatrixXd &J,Eigen::MatrixXd &grad){    fval = fx(x);    for (int i = 0;i<2;++i)    {        r(i,0) = rx(x,i);    }    for(int i = 0;i<2;++i)        for(int j = 0;j<2;++j)        {            J(i,j) = Jac(x,i,j);        }    grad = J.transpose() * r;}void LM(){    Eigen::MatrixXd x(2,1);    x(0,0) = -1.92;    x(1,0) = 2.0;    int max_iter = 100;    double func_tol = pow(10.0,-5);    int iter = 0;    double fval = INT_MAX;    Eigen::MatrixXd r(2,1);    Eigen::MatrixXd J(2,2);    Eigen::MatrixXd grad;    computeRelatedValues(x,fval,r,J,grad);    double nu = 0.01;    cout<<fval<<" ["<<x(0,0)<<" "<<x(1,0)<<"] "<<grad.norm()<<" "<<nu<<endl;    while(iter < max_iter && grad.norm() > func_tol)    {        ++iter;        computeRelatedValues(x,fval,r,J,grad);        Eigen::MatrixXd I = Eigen::MatrixXd::Identity(2,2);        Eigen::MatrixXd H = J.transpose() *J + nu*I;        Eigen::MatrixXd x_new = Eigen::MatrixXd::Zero(2,1);        int iter1 = 0;        double diff = (x-x_new).norm();        //注意如何比较两个向量是否相等        while(diff > pow(10.0,-5) && iter1 < max_iter)        {            ++iter1;            Eigen::MatrixXd delta_x = H.colPivHouseholderQr().solve(grad);            cout<<delta_x<<endl;            Eigen::MatrixXd x_new = x - delta_x;            cout<<x_new<<endl;            double fval_new = INT_MAX;            Eigen::MatrixXd r_new(2,1);            Eigen::MatrixXd J_new(2,2);            Eigen::MatrixXd grad_new;            computeRelatedValues(x_new,fval_new,r_new,J_new,grad_new);            double rho = abs(fval - fval_new) / (grad.transpose() * (x_new - x)).norm();            if (rho > 0)            {                x = x_new;                diff = (x - x_new).norm();                cout<<diff<<endl;                if (rho > 0.25)                {                    nu = nu /10;                }            }            else            {                nu = nu *10;            }        }         cout<<fval<<" ["<<x(0,0)<<" "<<x(1,0)<<"] "<<grad.norm()<<" "<<nu<<endl;    }}int main(){    LM();    return 0;}