// 目标函数: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;
}