牛顿法(newton's method)求解函数极值

来源：互联网发布：keeper是什么软件编辑：程序博客网时间：2024/04/19 21:32

// 目标函数：f(x1, x2, x3, x4) = (x1 + 10*x2)^2 + 5*(x3 - x4)^2 + (x2 - 2*x3)^4 + 10.0*(x1 - x4)^4; 求取输入向量x?// 牛顿法：  x(k+1) = x(k) - inv(F(x(k))) * g(x(k)), F是黑塞矩阵，g是导数;

double CPowell::ComputePowellFunc(float x[]){double dx[4] = {x[0], x[1], x[2], x[3]};double dResult = pow(dx[0] + 10 * dx[1], 2) + 5 * pow(dx[2] - dx[3], 2) + pow(dx[1] - 2 * dx[2], 4) + 10.0 * pow(dx[0] - dx[3], 4);return dResult;}

</pre><pre code_snippet_id="1835287" snippet_file_name="blog_20160817_4_2022945" name="code" class="html">

<span style="white-space:pre"></span>CPowell oPowell;// 初始值float x[4] = {3.0, -1.0, 0.0, 1.0};for (int i=0; i<3; i++){Eigen::VectorXf vecOld = Map<VectorXf>(x, 4);float fGraident[4] = {0.0};oPowell.ComputeGraident(x, 4, fGraident);Eigen::VectorXf vecGradient = Map<VectorXf>(fGraident, 4);float fHessian[4][4] = {0.0};oPowell.ComputeHessian(x, 4, &fHessian[0][0]);MatrixXf matHessian = Map<MatrixXf>(fHessian[0], 4, 4);   //动态矩阵，建立3行4列。Eigen::VectorXf vecNew = vecOld - matHessian.inverse() * vecGradient;cout<<"第"<<i<<"次迭代"<<endl;cout<<vecNew<<endl<<endl;for (int j=0; j<4; j++){x[j] = vecNew(j);}cout<<"函数值为："<<oPowell.ComputePowellFunc(x)<<endl<<endl;}

0 0