Armijo准则的最速下降法--C++和MATLAB实现

1、Armijo准则

前言

对于一个一般最速下降问题$$\mathrm{m}inf(x)$$
选择一个下降梯度方向不断下降最终下降到最小值点
下降过程为$$x^{k+1}=x^k-\tau ▽f(x)$$

通过设计$\tau$得到不同的收敛速度。对于$\tau$有如下4种常见方法$$\tau =c$$ $$\tau =c/k$$ τ = a r g min ⁡ α { x k + α d k } \tau=\mathrm {arg} \min\limits_{\alpha}\{x^k+\alpha dk\} τ=argαmin​{xk+αdk} τ ∈ { α ∣ f ( x k ) − f ( x k + α d ) ≥ − c α d T ▽ f ( x ) } ( c ∈ ( 0 , 1 ) ) \tau\in\{\alpha |f(x^k)-f(x^k+\alpha d)\geq-c\alpha d^\mathrm{T}\bigtriangledown f(x)\}(c\in(0,1)) τ∈{α∣f(xk)−f(xk+αd)≥−cαdT▽f(x)}(c∈(0,1))

前两种不能保证充分下降，第三种是精确搜索，虽然精确搜索迭代次数少，但是将高维问题转化为低维问题求解的时间仍然很长，因此工业上大家通常采用第四种Armijo准则。

定义

设 $d=-▽f(x^k)$是梯度的负方向 若

$f(x^k)-f(x^k+\alpha d)\ge -c\alpha d^{\mathrm{T}}▽f(x)$

则称步长 $\alpha$ 满足 Armjio 准则，其中 $c\in (0,1)$是一个常数。

几何含义

Armjio 准则有非常直观的几何含义，它指的是点 $(\alpha ,\varphi (\alpha ))$ 必须在直线$l_1(\alpha )=\varphi (0)+c\alpha \nabla f(x^k{)}^T{d}^k$的下方。如下图所示：

区间 $[0,{\alpha }_1]$中的点均满足 Armijo 准则。我们注意到
$d^k$为下降方向，这说明 $l(\alpha )$的斜率为负，选取符合条件 （1） 的 $\alpha$ 确实会使得函数值下降。

MATALB代码

Armjio.m 代码如下 Armjio.m。

unction [alpha, xk, f, k] = armjio(MyFun, grad, x0, dk)
% clear;clc;close all;
%% 非精确线性搜索求步长
beta = 0.333; 		    % 步长 alpha 的迭代系数，小于 1，越大越慢
c = 0.45; 		        % 泰勒展开式补足系数，0 < c < 1/2 ,越大越快（不绝对具体问题具体分析）
alpha = 1; 			    % 初始步长为 1
k = 0; 				    % 统计迭代次数
gk = feval(grad,x0);	% x0处的梯度值
fd = feval(MyFun,x0 + alpha * dk); 	            % 函数在下一个迭代点处的目标函数值
fk = feval(MyFun,x0) + alpha * c * gk' * dk; 	%  函数在下一个迭代点处的泰勒展开值
while fd > fk
    alpha = beta * alpha;
    fd = feval(MyFun,x0 + alpha * dk);
    fk = feval(MyFun,x0) + alpha * c * gk' * dk;
    k = k + 1;
end
xk = x0 + alpha * dk;	% 下降点
f = feval(MyFun, xk);	% 下降点处函数值

调用Armjio准则的最速下降 SD.m。

function [f, xk, k] = SD(x0, MyFun, grad, eps, kmax)
%%尝试可视化，只能实现二维等高线的可视化，对于多维的不行
if(length(x0)==2)
    for j=1:length(x0)
        X(j)=sym(['X',num2str(j)]);
    end
    F=0;
    for i=1:length(x0)
        if((mod(i,2) == 0))
            F=(100*(X(i-1)^2-X(i))^2+(X(i-1)-1)^2)+F;
        end
    end
    F = matlabFunction(F); % 将符号表达式转换为函数句柄或文件
    % 画出函数的水平集，然后hold on，待会会在这个图上显示出每一个迭代步长和迭代结果噢
    figure(1); clf; ezcontour(F,[-1.5 1.5 -1.5 1.5]); axis equal; hold on
end
%%
k = 0;                  % 统计迭代次数
xk = x0;
gk = feval(grad, xk)
while k <= kmax
    if norm(gk) < eps
        break;
    end
    dk = -gk;            % 选取此处的负梯度方向作为下降方向
    % 采用 Armjio 准则求步长
    alpha = armjio("MyFun", "grad", xk, dk)
    xk = xk + alpha * dk; % 更新迭代点
    gk = feval(grad, xk) % 更新梯度值
    plot([x0(1) xk(1)],[x0(2) xk(2)],'ko-');hold on
    refresh
    x0 = xk;
    k = k + 1;
end
f =feval(MyFun,xk);
end

我的function 和梯度函数，该function是一个任意维度的Rosenbrock函数.MyFun.m。

function fun=MyFun(x)
f=zeros(1,length(x));
for i=1:length(x)
    if((mod(i,2) == 0))
        f(i/2)=(100*(x(i-1)^2-x(i))^2+(x(i-1)-1)^2);
    end
end
fun=sum(f);
end

梯度函数grad.m。

function fun=grad(x)
G=zeros(length(x),1);
for i=1:length(x)
    if((mod(i,2) == 0))
        G(i-1,1)=(-400 * x(i-1) * x(i) + 400 * x(i-1)^3 + 2 * x(i-1) - 2);
        b=G(i-1,1)
        G(i,1)=( 200 * x(i) - 200 * x(i-1)^2);
    end
end
fun=G;%返回梯度矩阵

运行MATLAB代码

1、在该代码路径下，在命令行窗口输入：

clear； x0=randn(50,1);eps=1e-5;kmax=20000;
[f, xk, k] = SD(x0, "MyFun", "grad", eps, kmax)

结果

1、x0为50维时结果如下：

2、在2维的情况下实现了等高线的显示的可视化

C++代码

代码我放在了我的GitHub里链接如下: link