梯度法

迭代算法与收敛率

现在要考虑满足在定义域内一阶可微的多元实值函数 的局部极小点, 即

在各种各样的搜索极小点的方法中, 迭代算法是较为典型且常见的一种方法。迭代法是这样的一种算法, 选定某个(接近极小点)的点 (有时会把这称为极小点的估计), 以 为初始点, 通过某种规则产生一个点 并继续以点 采用同样的规则产生 以此类推, 理想情况下(实际操作中中在达到停机条件时会停止该过程)最终会产生一个无穷序列 每一步(第 次)的迭代利用的规则实际上是序列的某个递推关系式: 理想的情况是, 该迭代序列收敛于极小点 这样我们就成功地得到了一个搜索极小点的算法. 后面介绍的算法, 大都属于迭代算法.

不管是什么算法, 它的效率/运行速度总是我们最需要关注的性质之一, 显然, 对于迭代算法而言, 这与迭代产生的序列 的收敛速度息息相关, 于是我们借此定义收敛率, 并以收敛率衡量迭代算法的收敛速度.

定义 设序列 收敛到 若 则序列 的收敛阶数为 或 阶收敛.

特别地:

• 时称 超线性收敛;
• 时称 线性收敛;
• 时称 次线性收敛;
• 时称 是二次型收敛.

一般而言, 具有超线性及以上收敛速度的算法就可以被称为"很好的算法"了. 越大, 算法/数列的收敛速度就越快.

可以证明任意序列的收敛阶数不会小于 .

收敛速度还有其他的表示方法: 称序列 收敛到 的速度是 , 如果在保证 与 的误差控制在 内时得到的迭代次数满足 , 其中 是某个单调递减且 的函数, 这种定义方式与定义 1 的总结如下:

快慢关系是 . 同时为了证明某个数列的收敛速度, 只需要列出上述的递推关系即可.

最速下降法

在迭代算法中, 我们可以这样描述第 次迭代产生的新的点 是 上的两个点, 产生点 的过程可以看作是从点 出发沿着方向 前进 的距离, 借助可行方向, 我们将向量 作为一个可行方向, 记作 为了方便用范数和方向分别描述向量, 保证 前进的距离记作 这称为步长, 这样就有 我们的任务变成了通过某种规则产生"合适"的可行方向与步长.

什么样可行方向和步长的才算是"合适"呢? 这不难想, 既然我们的目的是求极小点, 那必然要保证在迭代的过程中序列 整体呈下降的趋势, 或者干脆让条件再强一些, 序列 严格递减, 那么我们便得到了下降迭代法.

但是这还不够具体, 我们要得到可行方向和步长的确切信息.

我们希望只利用可行方向和步长得到最好的下降迭代法, 那么我们显然需要下降最快的方向以及使 在该方向上取得最小值的步长. 那么方向的选择便呼之欲出了, 即负梯度方向, 于是我们可以令 可行方向选择负梯度的算法又被称作梯度法.

固定步长

在使用最优步长之前, 我们不妨"偷一下懒", 选择固定的步长, 这样就免去了每次迭代还要找最优步长的麻烦, 于是我们得到了固定步长梯度法. 这样做的好处比较明显, 缺点是:

• 步长较大在远离极小点的地方比较快但是容易在极小点附近产生锯齿状的收敛路径, 甚至不收敛;
• 步长太小效率会偏低.

事实上, 我们可以确定固定步长梯度法中使得算法收敛至极小点的步长的范围, 这放在后面介绍.

最优步长

显然最优步长 这是一个一维函数, 记作 极小值点在驻点处取得, 于是令 有 但是在实现算法的时候我们通常使用一维搜索方法.

确定了最优可行方向和最优步长后, 我们得到了最速下降梯度法:

1. 确定初始点 并置迭代次数
2. 令 置 其中 的计算可采用一维优化方法.
3. 当满足要求的精度时停止上述过程, 否则置 转 2.

最速下降法得到的序列一定会收敛至极小点, 证明过程在后文会有提到.

停机条件的选择

理论上来讲, 我们最好把 作为停机条件, 但是在实际中我们很难恰好得到梯度为零的结果, 只能得到某个接近 的结果, 此时若给定一个较小的预设的阈值 ,便有下面三种停机条件可以使用:

• 梯度范数: 利用 作为停机条件;

• 绝对数值: 利用 或 来作为停机条件.

• 相对数值: 利用 或 作为停机条件.

利用梯度范数和绝对数值作为停机条件的缺点是它们会受到 或 的数量级的影响, 而相对数值是尺度无关的. 然而若 或 时, 相对数值反而会比绝对数值大, 为避免这一情况我们可以采用下面的修正: 这也是后面的算法常用的停机条件.

经过测试利用修正过的相对数值作为停机条件比使用绝对数值速度快了至少十倍.

二次型函数的最速下降梯度法

设 有 于是 满足 整理可得 至于 它的形式并不是很漂亮, 就没有把它写出来的必要了.

收敛性质

固定步长

对于二次型函数 , 我们能确定使迭代序列收敛的步长范围, 这表述为下面的命题:

命题 对二次型函数 使用固定步长梯度法时, 当且仅当步长满足 时迭代产生的序列 收敛至极小点 其中 是 的最大特征值.

收敛性分析考虑证明两点内容:

• 算法产生的序列能收敛到驻点
• 算法产生的序列的收敛速度

当分析性质的时候一般要加一些理想的条件.

1: 首先 满足 Lipschitz 连续: 对任意 都有 . 这个连续条件对 的变化率做了限制. 2: 为光滑函数需要满足 Lipschitz 连续, 一个等价定义是 . 同时对梯度的 Lipschitz 连续考察类似于 Lagrange 中值定理的"推广", 可以不严谨地产生一种直觉: , 这里的矩阵范数是谱范数, 还可以写作 , 其中 为单位阵, 为半正定意义下的偏序关系, 即 半正定. 这样产生了光滑的三个等价条件. 3: 强凸 (strongly convex) 需要满足 为强凸系数.

考虑强凸且光滑的情形.

首先应用光滑函数的等价定义处理 , 即令 有 类似的, 按照强凸的条件取 , 有 事实上, 这对任意 都满足, 因此当 时, 其中限制 满足 , 那么 即 , 这实际上就是线性收敛的定义, 即收敛速度为 .

需要注意的是, 梯度下降法对于光滑非凸函数和光滑凸但非强凸的函数都是次线性收敛, 收敛速度分别为 和 .

代码实现

Python

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63

import math
from numpy import *
import sympy


def change(var, val):  #构建变量与值之间的以元组为元素的列表
    l = []
    it1 = iter(var)
    it2 = iter(val)
    for i in range(var.shape[0]):
        l.append(tuple((next(it1), next(it2))))
    return l


def golden_cut(f, l, r, var, eps):
    rho = (3-math.sqrt(5))/2
    a = l+(r-l)*rho
    b = r-(r-l)*rho
    while (math.fabs(l-r) > eps):
        f_a = f.subs([(var, a)])
        f_b = f.subs([(var, b)])
        if f_a[0, 0] > f_b[0, 0]:
            l = a
            a = b
            b = r-(r-l)*rho
        else:
            r = b
            b = a
            a = l+(r-l)*rho
    return a  #取左端点


def SteepestDescent(f, x0, var, eps):
    grad_f = f.jacobian(var)  # 雅可比矩阵
    x1 = sympy.Matrix.zeros(var.shape[0], 1)
    con = sympy.Matrix.norm(x1-x0)/max(sympy.Matrix.norm(x0), 1)
    k = 0
    while con > eps:
        rel1 = change(var, x0)
        d = grad_f.subs(rel1).T  # 这里要加转置，否则会出现维数不一致的情况
        alpha = sympy.symbols('alpha', positive = True)  #步长变量
        rel2 = change(var, x0-alpha*d)
        phi = f.subs(rel2)  #以步长为变量的函数
        alpha0 = golden_cut(phi, 0, 100, alpha, 1e-4)  #黄金分割法找最短步长
        x1 = x0-alpha0*d
        tmp = x1
        x1 = x0
        x0 = tmp
        con = sympy.Matrix.norm(x1-x0)/max(sympy.Matrix.norm(x0), 1)
        k += 1
    print(x0.evalf(6))
    rel3 = change(var, x0)
    min_f = f.subs(rel3)
    print(min_f[0, 0].evalf(6))
    print("Total iterations:", k)


x, y, z, w = sympy.symbols('x,y,z,w')
w = (x-5)**2+(y+4)**2+4*(z-6)**2
var = sympy.Matrix([x, y, z])
f = sympy.Matrix([w])
x0 = sympy.Matrix([1000, 1000, 1000])
SteepestDescent(f, x0, var, 1e-6)

输出结果为:

1
2
3

Matrix([[5.00000], [-4.00000], [6.00000]])
1.17900e-11
Total iterations: 24

MATLAB

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

% %[x,f_min]中,x为极小点,f_min为极小点对应的取值
% %f为目标函数,x0为初始点,var为变量,eps为精度
% %在使用该函数时,先定义目标函数的符号变量,写出目标函数,传入目标函数,初值点,变量以及精度即可
function [x,f_min]=SteepestDescent(f,x0,var,eps)
x1=0;
k=0;%迭代次数
grad_f=jacobian(f,var);%梯度
con=norm(x1-x0)./max(1,norm(x0));%相对值作为循环结束条件
while con>eps
    g=subs(grad_f,var,x0);
    syms alpha positive;
    phi=subs(f,var,x0-alpha*g); %phi是alpha的函数,求出phi的极小值
    [a,b]=GoldenCut(phi,alpha,0,2,1.0e-3); %利用黄金分割法求出极小点,精度为0.001
    alpha0=a;%取左端点
    x1=x0-alpha0*g;%下一个点
    tmp=x0;
    x0=x1;
    x1=tmp;
    k=k+1;
    con=norm(x1-x0)./max(1,norm(x0));%停机条件
end
x=vpa(x0,6);#保留六位小数
f_min=subs(f,var,x);
f_min=vpa(f_min,6);