多变量函数极值-最速下降法

简介

基本原理：假设$f(x)\in C^1$，由于函数$f(\mathbf{x})$沿着负梯度方向$-▽f(\mathbf{x})$下降最快，因此选取它作为每次迭代的搜索方向，即${\mathbf{s}}^k=-▽f({\mathbf{x}}^k)$。故梯度法的计算步骤是：给定一个初始点${\mathbf{x}}^0$，计算${\mathbf{s}}^0=-▽f({\mathbf{x}}^0)$，从${\mathbf{x}}^0$出发沿着${\mathbf{s}}^0$方向搜索得到$f(\mathbf{x})$的近似极小点${\mathbf{x}}^1={\mathbf{x}}^0+{\lambda }_0{\mathbf{s}}^0$，在以${\mathbf{x}}^1$为起点，重复刚才的过程直至满足收敛条件为至，至于为什么沿着负梯度方向可以使得函数值下降的最快，在相关的书籍中有严格的证明(本人也看得一知半解，就不班门弄斧了，不影响使用)。
最速下降法的优点是每次迭代的计算量和存储量都较小，并且从不好的初始点也能迭代出一个收敛的解，缺点是在极小点附近，收敛速度很慢，收敛速度与变量的尺度有较大的关系，并且搜索时步长确定的不准确，可能会破坏方法的收敛性。
方法的具体实施步骤如下：

1. 给定初始点$x^0$，收敛精度$\epsilon$，令$i=0$；
2. 计算$▽f({\mathbf{x}}^i)$；
3. 若$|▽f({\mathbf{x}}^i)|<\epsilon$,则迭代结束，取极小值点${\mathbf{x}}^{\ast }={\mathbf{x}}^i$；反之，则用之前介绍的一维搜索方法，求解$$min\phi (\lambda )=f({\mathbf{x}}^{(i)}+\lambda {\mathbf{s}}^{(i)})$$的极小值点$\lambda$，使得$f({\mathbf{x}}^{(i)}+\lambda {\mathbf{s}}^{(i)})<f({\mathbf{x}}^{(i)})$；
4. 令${\mathbf{x}}^{i+1}={\mathbf{x}}^{(i)}-\lambda ▽f({\mathbf{x}}^i)$，$i=i+1$，返回步骤2

优化问题

求解$f(x_1,x_2)=x_1^2+2x_2^2-4x_1-2x_1{x}_2$的极小值点，本文采用最速下降法来确定搜索方向，并采用之前介绍的不精确直接搜索方法确定搜索步长，设置初始点$x_0=[0,0]$，收敛精度$\epsilon =1e^{-3}$，计算程序如下：

import numpy as np

def cal_fun_f(x1, x2):
    """
    函数值计算函数
    """
    return x1**2+2*x2**2-4*x1-2*x1*x2

def partial_derivative_x1(x1, x2):
    """
    求解函数对x1的偏导数
    :param x1:
    :param x2:
    """
    return 2*x1-4-2*x2

def partial_derivative_x2(x1, x2):
    """
    求解函数对x2的偏导数
    :param x1:
    :param x2:
    """
    return 4*x2-2*x1

def cal_derivative_f1(x1, x2):
    """
    求解函数一阶导数值
    :param x1:
    :param x2:
    :return:
    """
    return np.array([partial_derivative_x1(x1, x2), partial_derivative_x2(x1, x2)]).reshape(-1, 1)

def wolfe_powell_rule1(f_x, f_x_, c1, lam, d_x, s):
    """
    判断准则1
    :param f_x:k点函数值
    :param f_x_:k+1点函数值
    :param c1:给定常数
    :param lam:步长
    :param d_x:k点一阶导数
    :param s:方向向量
    :return:
    """
    if (f_x-f_x_) >= (-c1*lam*np.transpose(d_x)[0, :]@s[:, 0]):
        return True
    else:
        return False

def wolfe_powell_rule2(d_x_, s, c2, d_x):
    """
    判断准则2
    :param d_x_: k+1点一阶导数
    :param s: k点方向导数
    :param c2: 给定常数
    :param d_x: k点一阶导数
    """
    if (np.transpose(d_x_)[0, :]@s[:, 0]) >= (c2*np.transpose(d_x)[0, :]@s[:, 0]):
        return True
    else:
        return False

def get_lambda(c1, c2, a, b, lam, x, s):
    """
    采用直接法基于wolf_powell准则得到每一步的搜索步长
    :param c1:给定阐述
    :param c2:给定参数
    :param a:步长区间左端点
    :param b:步长区间右端点
    :param lam:初始步长
    :param x: 初始点
    :param s: 方向向量
    """
    while True:
        f_x = cal_fun_f(x[0, 0], x[1, 0])
        d_x = cal_derivative_f1(x[0, 0], x[1, 0])
        x_ = x + lam * s  # 下一个迭代点
        f_x_ = cal_fun_f(x_[0, 0], x_[1, 0])
        d_x_ = cal_derivative_f1(x_[0, 0], x_[1, 0])
        rule1_judge = wolfe_powell_rule1(f_x, f_x_, c1, lam, d_x, s)
        rule2_judge = wolfe_powell_rule2(d_x_, s, c2, d_x)
        if rule1_judge and rule2_judge:
            return x_, lam
        else:
            if not rule1_judge:
                b = lam
                lam = (lam + a) / 2
            if rule1_judge and rule2_judge == False:
                a = lam
                lam = min(2 * lam, (lam + b) / 2)

def most_velocity_descent_f(x, gap, c1=0.1, c2=0.5, a=0, b=float("inf"), lam=1):
    """
    采用最速下降法确定搜索方向，采用直接法不精确方法不确定搜索步长
    """
    i = 0
    while True:
        s = -cal_derivative_f1(x[0, 0], x[1, 0])  # 计算方向向量
        x, lam = get_lambda(c1, c2, a, b, lam, x, s)  # 得到迭代步长并计算下一迭代点
        i += 1
        print(f"第{i}次迭代极小值点,迭代步长", (x[0, 0], x[1, 0]), lam)
        if np.sum(np.abs(cal_derivative_f1(x[0, 0], x[1, 0])), axis=0) < gap:
            break

if __name__ == "__main__":
    most_velocity_descent_f(np.array([0, 0]).reshape(-1, 1), 1e-3)

计算结果如下：

第1次迭代极小值点,迭代步长 (2.0, 0.0) 0.5
第2次迭代极小值点,迭代步长 (2.0, 1.0) 0.25
第3次迭代极小值点,迭代步长 (2.5, 1.0) 0.25
第4次迭代极小值点,迭代步长 (3.0, 1.5) 0.5
第5次迭代极小值点,迭代步长 (3.5, 1.5) 0.5
第6次迭代极小值点,迭代步长 (3.5, 1.75) 0.25
第7次迭代极小值点,迭代步长 (3.625, 1.75) 0.25
第8次迭代极小值点,迭代步长 (3.75, 1.875) 0.5
第9次迭代极小值点,迭代步长 (3.875, 1.875) 0.5
第10次迭代极小值点,迭代步长 (3.875, 1.9375) 0.25
第11次迭代极小值点,迭代步长 (3.90625, 1.9375) 0.25
第12次迭代极小值点,迭代步长 (3.9375, 1.96875) 0.5
第13次迭代极小值点,迭代步长 (3.96875, 1.96875) 0.5
第14次迭代极小值点,迭代步长 (3.96875, 1.984375) 0.25
第15次迭代极小值点,迭代步长 (3.9765625, 1.984375) 0.25
第16次迭代极小值点,迭代步长 (3.984375, 1.9921875) 0.5
第17次迭代极小值点,迭代步长 (3.9921875, 1.9921875) 0.5
第18次迭代极小值点,迭代步长 (3.9921875, 1.99609375) 0.25
第19次迭代极小值点,迭代步长 (3.994140625, 1.99609375) 0.25
第20次迭代极小值点,迭代步长 (3.99609375, 1.998046875) 0.5
第21次迭代极小值点,迭代步长 (3.998046875, 1.998046875) 0.5
第22次迭代极小值点,迭代步长 (3.998046875, 1.9990234375) 0.25
第23次迭代极小值点,迭代步长 (3.99853515625, 1.9990234375) 0.25
第24次迭代极小值点,迭代步长 (3.9990234375, 1.99951171875) 0.5

可以看到在极小值点附近的时候，收敛的速度很慢。