求解非线性最小二乘法 Eigen

[html] view plain copy

print?

1. // 利用Eigen 求解非线性最小二乘;  

[html] view plain copy

print?

1. // 示例：<span style="font-family: Arial, Helvetica, sans-serif;">y = 10*(x0+3)^2 + (x1-5)^2</span><pre name="code" class="html"><span style="font-family: Arial, Helvetica, sans-serif;">#include "math.h"</span>  

[html] view plain copy

print?

1. #include "iostream"  
2. #include "vector"  
3. #include "list"  
4.   
5. using namespace std;  
6.   
7. #include "Eigen/Dense"  
8. #include "Eigen/Core"  
9. #include <unsupported/Eigen/NonLinearOptimization>  
10. #include <unsupported/Eigen/NumericalDiff>  
11.   
12. using namespace Eigen;  

[html] view plain copy

print?

1. // Generic functor  
2. template<typename _Scalar, int NX = Eigen::Dynamic, int NY = Eigen::Dynamic>  
3. struct Functor  
4. {  
5.     typedef _Scalar Scalar;  
6.     enum {  
7.         InputsAtCompileTime = NX,  
8.         ValuesAtCompileTime = NY  
9.     };  
10.     typedef Eigen::Matrix<Scalar,InputsAtCompileTime,1> InputType;  
11.     typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,1> ValueType;  
12.     typedef Eigen::Matrix<Scalar,ValuesAtCompileTime,InputsAtCompileTime> JacobianType;  
13.   
14.     int m_inputs, m_values;  
15.   
16.     Functor() : m_inputs(InputsAtCompileTime), m_values(ValuesAtCompileTime) {}  
17.     Functor(int inputs, int values) : m_inputs(inputs), m_values(values) {}  
18.   
19.     int inputs() const { return m_inputs; }  
20.     int values() const { return m_values; }  
21.   
22. };  

[html] view plain copy

print?

1. struct my_functor : Functor<double>  
2. {  
3.     // 输出个数必须大于输入个数, 故用2不用1;  
4.     my_functor(void): Functor<double>(2, 2) {}  
5.     int operator()(const Eigen::VectorXd &x, Eigen::VectorXd &fvec) const  
6.     {  
7.         // Implement y = 10*(x0+3)^2 + (x1-5)^2  
8.         fvec(0) = 10.0*pow(x(0)+3.0,2) +  pow(x(1)-5.0,2);  
9.         fvec(1) = 0;  
10.   
11.         return 0;  
12.     }  
13. };  

[html] view plain copy

print?

1. int main(int argc, char *argv[])  
2. {  
3.     Eigen::VectorXd x(2);  
4.     x(0) = 1.0;  
5.     x(1) = 3.0;  
6.   
7.     my_functor functor;  
8.     Eigen::NumericalDiff<my_functor> numDiff(functor);  
9.     Eigen::LevenbergMarquardt<Eigen::NumericalDiff<my_functor>,double> lm(numDiff);  
10.   
11.     Eigen::VectorXd y(2);  
12.     functor.operator()(x, y);  
13.   
14.     std::cout << "x first input: \n" << x << std::endl;  
15.     std::cout<<"y first outpout: \n" << y << std::endl;  
16.   
17.     lm.parameters.maxfev = 1000;  
18.     lm.parameters.xtol = 1.0e-8;  
19.   
20.     int iRet = lm.minimize(x);  
21.     std::cout << "迭代次数：\n"<<lm.iter << std::endl;  
22.     std::cout << "计算标志：\n" << iRet << std::endl;  
23.   
24.     std::cout << "x finnal: \n" << x << std::endl;  
25.   
26.     functor.operator()(x, y);  
27.     std::cout<<"y outpout((minimized): \n" << y << std::endl;  
28.   
29.     return 0;  
30. }  

// 最终看到输出了x = [-3.0, 5.0]. 使得目标最小!