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

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