梯度下降算法

一、梯度下降算法
本文介绍梯度下降算法，其本质是仅仅使用函数的一阶导数信息选取下降方向$d_k$。梯度下降法的方向选取非常直观，实际应用范围非常广，因此它在优化算法中的地位可相当于高斯消元法在线性方程组中的地位。此外我们也会介绍BB方法，该方法作为一种梯度法的变形，虽然理论性质目前仍不完整，但由于优秀的数值表现，也是在实际应用中使用较多的一种算法。
对于光滑函数$f(x)$，在迭代点$x_k$处，我们需要选择一个较为合理的$d_k$作为下降方向。注意到$\varphi (\alpha )=f(x_k+\alpha d_k)$有泰勒展开
$$\varphi (\alpha )=f(x_k)+\alpha \nabla f(x_k{)}^T{d}_k+O({\alpha }^2||d_k|{|}^2)$$
根据柯西不等式，当$\alpha$足够小时，取$d_k=-\nabla f(d_k)$会使得函数下降最快。因此梯度法就是选取$d_k=-\nabla f(x_k)$的算法，它的迭代格式为：
$$x_{k+1}=x_k-{\alpha }_k\nabla f(x_k).$$
其中步长${\alpha }_k$的选取可依赖于一维线搜索，也可直接选取固定的步长。

二、梯度下降法的收敛性
本节不加证明的给出梯度下降发的收敛性定理：
定理1（梯度法在凸函数上的收敛性）：设函数$f(x)$为凸的梯度L-利普希茨连续函数，$f^{\ast }=f(x^{\ast })=\underset{x}{\inf }f(x)$存在且可达。如果步长${\alpha }_k$取为常数$\alpha$且满足$0<\alpha \le \frac{1}{L}$，那么梯度下降的电联 { x k } \{x_k\} {xk​}的函数值收敛到最优值，且在函数值的意义下收敛速度为$O(\frac{1}{k})$

三、Barzilar-Borwein方法
Barzilar-Borwein方法又简称BB法，从形式上来看BB法的下降方向仍然是点$x_k$处的负梯度方向$-\nabla f(x_k)$，但步长${\alpha }_k$并不是直接由线搜索算法给出的。考虑梯度下降法的格式
$$x_{k+1}=x_k-{\alpha }_k\nabla f(x_k)$$
这种格式也可以写成
$$x_{k+1}=x_k-D_k\nabla f(x_k)$$
其中$D_k={\alpha }_{k}I$，BB法选取的${\alpha }_k$是如下两个最优问题之一的解：
$$\begin{array}{rl}& \underset{\alpha }{\min }||\alpha y_{k-1}-s_{k-1}|{|}^2,\\ & \underset{\alpha }{\min }||y_{k-1}-{\alpha }^{-1}{s}_{k-1}|{|}^2\end{array}$$
其中我们引入记号$s_{k-1}=x_k-x_{k-1}$以及$y_{k-1}=\nabla f(x_k)-\nabla f(x_{k-1})$，容易求得上述问题的解分别为${\alpha }_{BB1}^k=\frac{(s_{k-1}{)}^T{y}_{k-1}}{(y_{k-1}{)}^T{y}_{k-1}}$和${\alpha }_{BB2}^k=\frac{(s_{k-1}{)}^T{s}_{k-1}}{(s_{k-1}{)}^T{y}_{k-1}}$，因此可以得到BB方法的两种迭代格式：
$$\begin{array}{rl}{x}_{k+1}& =x_k-{\alpha }_{BB1}^k\nabla f(x_k),\\ x_{k+1}& =x_k-{\alpha }_{BB2}^k\nabla f(x_k).\end{array}$$
从两种BB步长的计算公式可知，任何一种仅仅需要函数和相邻亮部的梯度信息和迭代点信息，不需要任何线搜索算法即可选取算法步长。正因为这个特点，BB算法的使用范围特别广泛。对于一般的问题，通过BB步长计算公式得到的步长可能过大或者过小，因此我们还需要将步长做上界和下界的截断，即选取$0<{\alpha }_m<{\alpha }_M$使得${\alpha }_m\le {\alpha }_k\le {\alpha }_M$，还需要注意的是，BB方法本身是非单调方法，有时也配合非单调收敛准则使用以获得更好的实际效果。

三、梯度下降法应用举例
应用举例1
对于函数$f(x,y)=x^2+10y^2$，比较固定步长梯度下降算法和BB步长梯度下降算法的计算效率：
固定步长梯度下降法matlab实现

function [fmin, xmin, fk, xk] = gradient_descent_fix_stepsize(func, gfunc, x0, alpha, epsilon)

iIter = 1;
iterMax = 500;
xOld = x0;
xk = zeros(size(x0, 1), 66);
xk(:, 1) = x0;
fk = zeros(1, 66);

while iIter < iterMax
    dk = feval(gfunc, xOld);
    xNew = xOld - alpha * dk;
    
    if norm(xNew - xOld, 2) <= epsilon
        break;
    end
    
    iIter = iIter + 1;
    
    xk(:, iIter) = xNew;
    fk(:, iIter) = feval(func, xNew);
    
    xOld = xNew;
end


if iIter == iterMax
    fprintf('reach maximum iteration, and not found minimal x!\n');
end

xmin = xNew;
fmin = feval(func, xmin);
fprintf('iIter = %d, fmin = %f\n', iIter, fmin);   



end

BB步长梯度下降法matlab实现

function [fmin, xmin, fk, xk] = gradient_descent(func, gfunc, x0, epsilon)

iIter = 1;
iterMax = 500;
xOld = x0;

alphaMin = 1e-5;
alphaMax = 1e5;
M = 10;
%alpha = 0.5;
[~, ~, alpha] = armijo_rule(func, gfunc, x0, 0.5, -feval(gfunc, x0));

QOld = 1;
COld = feval(func, xOld);
c1 = 0.5;
eta = 0.5;
xk = zeros(size(x0, 1), 11);
fk = zeros(11, 1);
xk(:, 1) = x0;
fk(1, :) = feval(func, x0);

while iIter < iterMax
    grad = feval(gfunc, xOld);
    dk = -grad;
    
    % Zhang, Hanger nonmonotone line search
    for i = 1:M
        xNew = xOld + alpha * dk;
        fNew = feval(func, xNew); 
        if fNew <= COld + alpha * c1 * dot(dk, dk)
            break;
        end
        alpha = alpha * eta;
    end
    
    %xNew = xOld - alpha * dk;
    iIter = iIter + 1;
    
    if norm(grad, 2) < epsilon
        break;
    end
    
    % BB step-size calculation
    sk = xNew - xOld; yk = feval(gfunc, xNew) - feval(gfunc, xOld);
    if mod(iIter, 2) == 1
        alpha = dot(sk, yk) / dot(yk, yk);
    else
        alpha = dot(sk, sk) / dot(sk, yk);
    end
    
    alpha = max(min(alpha, alphaMax), alphaMin);
    
    QNew = eta * QOld + 1;
    CNew = (eta * QOld * COld + fNew) / QNew;
    COld = CNew;
    
    xOld = xNew;
    xk(:, iIter) = xNew;
    fk(iIter, :) = fNew;
end

if iIter == iterMax
    fprintf('exceed max iteration, not found  minimal point x.\n');
end

xmin = xNew;
fmin = feval(func, xmin);
fprintf('iIter = %d, fmin = %f\n', iIter, fmin);

end

这里用到的armijo准则实现见前面的线搜索部分
测试比较二者计算结果：

close all
% f(x, y) = x^2 + 10y^2;

func = @(x)(x(1)^2 + 10 * x(2)^2);
gfunc = @(x)([2*x(1); 20 * x(2)]);
x0 = [10; 1];
alpha = 0.085;
epsilon = 1e-5;
[fmin, xmin, fk, xk] = gradient_descent_fix_stepsize(func, gfunc, x0, alpha, epsilon);
x1 = -12:1e-2:12;
x2 = -10:1e-2:10;
[X1, X2] = meshgrid(x1, x2);
F = X1.^2 + 10 * X2.^2;
figure, contour(X1, X2, F, 50)
hold on
plot(xk(1,:), xk(2,:), 'LineWidth', 2)
plot(xk(1,:), xk(2,:), 'o', 'LineWidth', 2)


x0 = [-10; -1];
[fmin, xmin, fkBB, xkBB] = gradient_descent(func, gfunc, x0, epsilon);
hold on
plot(xkBB(1, :), xkBB(2, :), 'k', 'LineWidth', 2);
plot(xkBB(1, :), xkBB(2, :), 'bo', 'LineWidth', 2);

经过计算，固定步长梯度下降需要66次迭代，而BB步长梯度下降只需要12步迭代就达到最优点(0,0), fmin = 0，从等高线图来看，BB步长梯度下降是非单调线搜索一类的下降算法，并不是每一步都严格下降的。