Newton法程序设计

一、实验目的
掌握Hesse矩阵的计算方法和Newton法的基本思想及其迭代步骤；学会运用MATLAB编程实现常用优化算法；能够正确处理实验数据和分析实验结果及调试程序。
二、实验内容

（1）求解无约束优化问题： ；
（2）终止准则取 ；
（3）完成Newton法（牛顿法）的MATLAB编程、调试；
（4）选取初始点，并给出相关实验结果的对比及分析（从最优解、最优值、收敛速度（迭代次数）等方面进行比较）；
（5）按照模板撰写实验报告，要求规范整洁。
三、算法步骤、代码、及结果
1. 算法步骤
1. 初始化参数：
定义初始点集合 initial_points。
定义迭代终止准则 tol = 1e-6。
2. 遍历所有初始点：
对于每个初始点 x0，调用 newton_method 函数来进行牛顿法迭代。
3. 定义目标函数 f(x)：
计算给定 x 处的目标函数值。
4. 定义梯度函数 grad_f(x)：
计算给定 x 处的目标函数的梯度。
5. 定义 Hessian 矩阵函数 hessian_f(x)：
计算给定 x 处的目标函数的 Hessian 矩阵。
6. 牛顿法迭代 newton_method：
初始点为 x0。
当梯度的范数大于终止准则 tol 时，重复以下步骤：
计算当前 x 处的 Hessian 矩阵 hess。
计算当前 x 处的梯度 grad。
更新 x 的值：x = x - hess \ grad，即利用 Hessian 矩阵的逆矩阵乘以梯度进行更新。
增加迭代计数器 k。
当梯度的范数小于或等于 tol 时，停止迭代，输出当前 x 作为最优解，同时输出最优函数值 fval 和迭代次数 k。
7. 输出结果：
对每个初始点，输出最优解 x_opt，最优函数值 fval 和迭代次数 iterations。

2. 代码
function newton_optimization()
% 定义初始点
initial_points = [
1, 1, 1, 1;
-2, 2, -1, 2;
0, 0, 0, 0
];

% 迭代终止准则
tol = 1e-6;

% 遍历所有初始点
for i = 1:size(initial_points, 1)
x0 = initial_points(i, :)';
[x_opt, fval, iterations] = newton_method(@f, @grad_f, @hessian_f, x0, tol);
fprintf('Initial Point %d: \n', i);
fprintf('Optimal x: %s\n', mat2str(x_opt));
fprintf('Optimal function value: %.6f\n', fval);
fprintf('Number of iterations: %d\n\n', iterations);
end
end

function val = f(x)
val = (x(1) + 10 * x(2))^2 + 5 * (x(3) - x(4))^2 + (x(2) - 2 * x(3))^4 + 10 * (x(1) - x(4))^4;
end

function grad = grad_f(x)
grad = zeros(4, 1);
grad(1) = 2 * (x(1) + 10 * x(2)) + 40 * (x(1) - x(4))^3;
grad(2) = 20 * (x(1) + 10 * x(2)) + 4 * (x(2) - 2 * x(3))^3;
grad(3) = 10 * (x(3) - x(4)) - 8 * (x(2) - 2 * x(3))^3;
grad(4) = -10 * (x(3) - x(4)) - 40 * (x(1) - x(4))^3;
end

function hess = hessian_f(x)
hess = zeros(4, 4);
hess(1, 1) = 2 + 120 * (x(1) - x(4))^2;
hess(1, 2) = 20;
hess(1, 4) = -120 * (x(1) - x(4))^2;

hess(2, 1) = 20;
hess(2, 2) = 200 + 12 * (x(2) - 2 * x(3))^2;
hess(2, 3) = -24 * (x(2) - 2 * x(3))^2;

hess(3, 2) = -24 * (x(2) - 2 * x(3))^2;
hess(3, 3) = 10 + 48 * (x(2) - 2 * x(3))^2;
hess(3, 4) = -10;

hess(4, 1) = -120 * (x(1) - x(4))^2;
hess(4, 3) = -10;
hess(4, 4) = 10 + 120 * (x(1) - x(4))^2;
end

function [x, fval, k] = newton_method(f, grad_f, hessian_f, x0, tol)
x = x0;
k = 0;
while norm(grad_f(x)) > tol
hess = hessian_f(x);
grad = grad_f(x);
% 更新步骤 x = x - H^-1 * grad
x = x - hess \ grad;
k = k + 1;
end
fval = f(x);
end

3. 结果
Initial Point 1: Iterations = 14, Optimal Value = 0.000000, Optimal Point = [0.002333, -0.000233, 0.001596, 0.001596]
Initial Point 2: Iterations = 19, Optimal Value = 0.000000, Optimal Point = [-0.002578, 0.000258, -0.000773, -0.000773]
Initial Point 3: Iterations = 0, Optimal Value = 0.000000, Optimal Point = [0.000000, 0.000000, 0.000000, 0.000000]