算法描述及示例引用自“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;
}