【Matlab算法】拟牛顿法(Quasi-Newton Methods)(附MATLAB完整代码)

发布时间:2023年12月20日

拟牛顿法(Quasi-Newton Methods)

前言

拟牛顿法是一类迭代优化算法,用于求解无约束优化问题。与牛顿法类似,拟牛顿法的目标是通过迭代逼近目标函数的最优解,但是它不显式计算目标函数的二阶导数
(Hessian矩阵) 。相反,它通过逐步构建一个拟牛顿矩阵 (Quasi-Newton Matrix) 来模拟Hessian矩阵的逆。

以下是拟牛顿法的基本思想和步骤:

  1. 初姶伙参数:选择一个初始点 x ( 0 ) x^{(0)} x(0) 作为优化的起始点,并初始化一个初始的拟牛顿矩阵 B ( 0 ) B^{(0)} B(0) (通常选择单位矩阵)。
  2. 选代过程:
  • 对于每次迭代 t t t :
  • 计算目标函数 f ( x ( t ) ) f\left(x^{(t)}\right) f(x(t)) 在当前点 x ( t ) x^{(t)} x(t) 处的梯度 ? f ( x ( t ) ) \nabla f\left(x^{(t)}\right) ?f(x(t))
  • 使用拟牛顿矩阵 B ( t ) B^{(t)} B(t) 近似Hessian矩阵的逆,得到搜索方向 δ x ( t ) \delta x^{(t)} δx(t) :
    δ x ( t ) = ? B ( t ) ? f ( x ( t ) ) \delta x^{(t)}=-B^{(t)} \nabla f\left(x^{(t)}\right) δx(t)=?B(t)?f(x(t))
  • 选择一个合适的步长 α ( t ) \alpha^{(t)} α(t) ,通常通过线搜索来确定。
  • 更新参数: x ( t + 1 ) = x ( t ) + α ( t ) δ x ( t ) x^{(t+1)}=x^{(t)}+\alpha^{(t)} \delta x^{(t)} x(t+1)=x(t)+α(t)δx(t)
  • 更新拟牛顿矩阵 B ( t + 1 ) B^{(t+1)} B(t+1) ,以使其逼近目标点 x ( t + 1 ) x^{(t+1)} x(t+1) x ( t ) x^{(t)} x(t) 处的梯度变化。拟牛顿矩阵的更新规则是拟合Hessian矩阵逆的变化。
    B ( t + 1 ) = B ( t ) + y ( t ) ( y ( t ) ) T ( y ( t ) ) T s ( t ) ? B ( t ) s ( t ) ( B ( t ) s ( t ) ) T ( s ( t ) ) T B ( t ) s ( t ) B^{(t+1)}=B^{(t)}+\frac{y^{(t)}\left(y^{(t)}\right)^T}{\left(y^{(t)}\right)^T s^{(t)}}-\frac{B^{(t)} s^{(t)}\left(B^{(t)} s^{(t)}\right)^T}{\left(s^{(t)}\right)^T B^{(t)} s^{(t)}} B(t+1)=B(t)+(y(t))Ts(t)y(t)(y(t))T??(s(t))TB(t)s(t)B(t)s(t)(B(t)s(t))T?

其中, s ( t ) s^{(t)} s(t) 是参数的更新量, y ( t ) y^{(t)} y(t) 是梯度的更新量。
3. 停止条件:检查是否满足停止条件。可能的停止条件包括:

  • 达到预定的迭代次数。
  • 梯度的范数小于某个容许误差。
  1. 输出结果:输出最终的参数 x ( t ) x^{(t)} x(t) ,以及在最优点的目标函数值 f ( x ( t ) ) f\left(x^{(t)}\right) f(x(t))

拟牛顿法的优点在于它避免了显式计算目标函数的Hessian矩阵,从而减少了计算的复杂性。不同的拟牛顿法算法可能使用不同的矩阵更新规则,例如DFP算法和BFGS算法是两种常用的拟牛顿法算法。这些方法通过维护一个逐步逼近的Hessian逆矩阵来实现参数的更新。

正文

使用拟牛顿法计算目标函数 f ( x ) = x ( 1 ) 2 + x ( 2 ) 2 ? 2 x ( 1 ) x ( 2 ) + sin ? ( x ( 1 ) ) + f(x)=x(1)^2+x(2)^2-2 x(1) x(2)+\sin (x(1))+ f(x)=x(1)2+x(2)2?2x(1)x(2)+sin(x(1))+ cos ? ( x ( 2 ) ) \cos (x(2)) cos(x(2)) 的过程如下:

  1. 初始化参数: 选择一个初始点 x ( 0 ) x^{(0)} x(0) 作为优化的起始点,并初始化一个初始的拟牛顿矩阵 B ( 0 ) B^{(0)} B(0) (通常选择单位矩阵)。
  2. 选代这程:
  • 对于每次迭代 t t t :
  • 计算目标函数在当前点 x ( t ) x^{(t)} x(t) 处的梯度 ? f ( x ( t ) ) \nabla f\left(x^{(t)}\right) ?f(x(t))
    ? f ( x ( t ) ) = [ 2 x 1 ? 2 x 2 + cos ? ( x 1 ) 2 x 2 ? 2 x 1 ? sin ? ( x 2 ) ] \nabla f\left(x^{(t)}\right)=\left[\begin{array}{l} 2 x_1-2 x_2+\cos \left(x_1\right) \\ 2 x_2-2 x_1-\sin \left(x_2\right) \end{array}\right] ?f(x(t))=[2x1??2x2?+cos(x1?)2x2??2x1??sin(x2?)?]
  • 使用拟牛顿矩阵 B ( t ) B^{(t)} B(t) 近似Hessian矩阵的逆,得到捜索方向 δ x ( t ) \delta x^{(t)} δx(t) :
    δ x ( t ) = ? B ( t ) ? f ( x ( t ) ) \delta x^{(t)}=-B^{(t)} \nabla f\left(x^{(t)}\right) δx(t)=?B(t)?f(x(t))
  • 选择一个合适的步长 α ( t ) \alpha^{(t)} α(t) ,通常通过线搜索来确定。
  • 更新参数: x ( t + 1 ) = x ( t ) + α ( t ) δ x ( t ) x^{(t+1)}=x^{(t)}+\alpha^{(t)} \delta x^{(t)} x(t+1)=x(t)+α(t)δx(t)
  • 计算参数的更新量 s ( t ) s^{(t)} s(t) 和梯度的更新量 y ( t ) y^{(t)} y(t)
    s ( t ) = x ( t + 1 ) ? x ( t ) y ( t ) = ? f ( x ( t + 1 ) ) ? ? f ( x ( t ) ) \begin{aligned} & s^{(t)}=x^{(t+1)}-x^{(t)} \\ & y^{(t)}=\nabla f\left(x^{(t+1)}\right)-\nabla f\left(x^{(t)}\right) \end{aligned} ?s(t)=x(t+1)?x(t)y(t)=?f(x(t+1))??f(x(t))?
  • 更新拟牛顿矩阵 B ( t + 1 ) B^{(t+1)} B(t+1) 使用拟牛顿法的矩阵更新规则:
    B ( t + 1 ) = B ( t ) + y ( t ) ( y ( t ) ) T ( y ( t ) ) T s ( t ) ? B ( t ) s ( t ) ( B ( t ) s ( t ) ) T ( s ( t ) ) T B ( t ) s ( t ) B^{(t+1)}=B^{(t)}+\frac{y^{(t)}\left(y^{(t)}\right)^T}{\left(y^{(t)}\right)^T s^{(t)}}-\frac{B^{(t)} s^{(t)}\left(B^{(t)} s^{(t)}\right)^T}{\left(s^{(t)}\right)^T B^{(t)} s^{(t)}} B(t+1)=B(t)+(y(t))Ts(t)y(t)(y(t))T??(s(t))TB(t)s(t)B(t)s(t)(B(t)s(t))T?
  1. 停止条件: 检查是否满足停止条件。可能的停止条件包括:
  • 达到预定的迭代次数。
  • 梯度的范数小于某个容许误差。
  1. 输出结果:输出最终的参数 x ( t ) x^{(t)} x(t) ,以及在最优点的目标函数值 f ( x ( t ) ) f\left(x^{(t)}\right) f(x(t))
    这个过程中,拟牛顿矩阵 B ( t ) B^{(t)} B(t) 不断地被更新以逼近Hessian矩阵的逆,从而改善参数的捜索方向。这使得拟牛顿法具有更广泛的适用性,尤其是在高维问题中。

代码实现

可运行代码

% 定义目标函数
f = @(x) x(1)^2 + x(2)^2 - 2*x(1)*x(2) + sin(x(1)) + cos(x(2));

% 定义目标函数的梯度
grad_f = @(x) [2*x(1) - 2*x(2) + cos(x(1)); 2*x(2) - 2*x(1) - sin(x(2))];

% 设置参数
max_iterations = 20;
tolerance = 1e-6;

% 初始化起始点和拟牛顿矩阵
x = [0; 0];
B = eye(2);  % 初始拟牛顿矩阵为单位矩阵

% 存储迭代过程中的参数和目标函数值
history_x = zeros(2, max_iterations);
history_f = zeros(1, max_iterations);

% 拟牛顿法迭代
for iteration = 1:max_iterations
    % 计算梯度
    gradient = grad_f(x);
    
    % 使用拟牛顿矩阵近似Hessian矩阵的逆,得到搜索方向
    delta_x = -B * gradient;
    
    % 选择合适的步长,通常通过线搜索确定
    alpha = 0.1;  % 这里简化为常数步长,实际应用中需要进行线搜索
    % 更新参数
    x_new = x + alpha * delta_x;
    
    % 计算参数的更新量和梯度的更新量
    s = x_new - x;
    y = grad_f(x_new) - gradient;
    
    % 更新拟牛顿矩阵
    B = B + (y * y') / (y' * s) - (B * s * (B * s)') / (s' * B * s);
    
    % 更新参数
    x = x_new;
    
    % 存储迭代过程中的参数和目标函数值
    history_x(:, iteration) = x;
    history_f(iteration) = f(x);
    
    % 检查停止条件
    if norm(gradient) < tolerance
        break;
    end
end

% 可视化迭代过程
figure;
subplot(2, 1, 1);
plot(1:iteration, history_x(1, 1:iteration), '-o', 'LineWidth', 1.5);
hold on;
plot(1:iteration, history_x(2, 1:iteration), '-o', 'LineWidth', 1.5);
title('参数迭代过程');
legend('x(1)', 'x(2)');
xlabel('迭代次数');
ylabel('参数值');

subplot(2, 1, 2);
plot(1:iteration, history_f(1:iteration), '-o', 'LineWidth', 1.5);
title('目标函数值迭代过程');
xlabel('迭代次数');
ylabel('目标函数值');

% 显示最终结果
fprintf('最优解: x = [%f, %f]\n', x(1), x(2));
fprintf('f(x)的最优值: %f\n', f(x));
fprintf('迭代次数: %d\n', iteration);

迭代结果

请添加图片描述

文章来源:https://blog.csdn.net/AlbertDS/article/details/135095506
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。