chapter5.拟牛顿法

拟牛顿法

pip install optimtool来使用优化算法

optimtool源码

# 导包
import sympy as sp
import matplotlib.pyplot as plt
import optimtool as oo
​

# 可视化函数：传参接口（颜色列表，函数值列表，标题列表）
def test(colorlist, f_list, title):
    handle = []
    for j, z in zip(colorlist, f_list):
        ln, = plt.plot([i for i in range(len(z))], z, c=j, marker='o', linestyle='dashed')
        handle.append(ln)
    plt.xlabel("$Iteration \ times \ (k)$")
    plt.ylabel("$Objective \ function \ value: \ f(x_k)$")
    plt.legend(handle, title)
    plt.title("Performance Comparison")
    plt.show()
    return None
​

1. 牛顿法

  牛顿法在理论上和实践中均中均取得很好的效果。对于二次连续可微函数 f(x) ，考虑 f(x) 在迭代点 x^k 处的二阶泰勒展开：

我们的目的是根据这个二阶近似来选取合适的下降方向 d^k 。如果忽略上式中的高阶项，并将等式右边看成关于 d^k 的函数求其稳定点，可以得到：

该方程也被称为牛顿方程，容易得出当 \nabla^2 f(x^k) 非奇异时，更新方向为 d^k=-\nabla f^2(x^k)^{-1} \nabla f(x^k) 。 一般称满足该方程的 d^k 为牛顿方向。 因此经典牛顿法的更新格式为：

步长 \alpha_k 恒为1，即可以不额外考虑步长的选取。我们也称步长为1的牛顿法为经典牛顿法。

牛顿法测试函数：

def dataSet1():
    f = sp.symbols("f")
    x1, x2 = sp.symbols("x1 x2")
    f = (1 - x1)**2 + (x2 - x1**2)**2
    funcs = sp.Matrix([f])
    args = sp.Matrix([x1, x2])
    x_0 = (2, 2)
    return funcs, args, x_0
​

  修正牛顿法主要修正了牛顿方程中的海瑟矩阵 \nabla^2 f(x^k) ，使其变成正定矩阵；同时引入线搜索以改善算法稳定性。

  在经典牛顿法中，计算牛顿方向d^k依赖于求解线性方程组。当 n 较大但 \nabla^2 f(x^k) 有稀疏结构时，需要通过迭代法来求解牛顿方程。以下是非精确牛顿法的函数程序。

def train1():
    funcs, args, x_0 = dataSet1()
    f_list = []
    title = ["newton_classic", "newton_modified", "newton_CG"]
    colorlist = ["maroon", "teal", "slateblue"]
    _, _, f = oo.unconstrain.newton.classic(funcs, args, x_0, False, True)
    f_list.append(f)
    _, _, f = oo.unconstrain.newton.modified(funcs, args, x_0, False, True)
    f_list.append(f)
    _, _, f = oo.unconstrain.newton.CG(funcs, args, x_0, False, True)
    f_list.append(f)
    return colorlist, f_list, title
​

colorlist, f_list, title = train1()
test(colorlist, f_list, title)
​

  然而对于大规模问题，函数的海瑟矩阵计算代价大或者难以得到，即使得到海瑟矩阵我们还需要求解一个大规模线性方程组，那么能否使用海瑟矩阵或其逆矩阵的近似来进行牛顿迭代呢？拟牛顿法便是这样的算法，它能够在每一步以较小的计算代价生成近似矩阵，并且使用近似矩阵代替海瑟矩阵而产生的迭代序列仍具有超线性收敛的性质．

  拟牛顿方法不计算海瑟矩阵 \nabla^2 f(x) ，而是构造其近似矩阵 B^k 或其逆的近似矩阵 H^k 。我们希望 B^k 或 H^k 仍然保留海瑟矩阵的部分性质。例如使得 d^k 仍然为下降方向。那么拟牛顿矩阵应该满足一些性质？如何构造它们呢？

2. 割线方程

  设 f(x) 是二次连续可微函数，根据泰勒展开，向量值函数 \nabla f(x) 在 x^{k+1} 处的近似为：

令 x=x^k ，s^k=x^{k+1}-x^k 及 y^k=\nabla f(x^{k+1})-\nabla f(x^k) ，得到：

忽略高阶项 ||s^k||^2|| ，我们希望海瑟矩阵 \nabla^2 f(x^{k+1}) 的近似矩阵 B^{k+1} 满足方程：

或者其逆的近似矩阵 H^{k+1} 满足方程：

并称上述 (2.3) 与 (2.4) 两式为割线方程。

  注意到近似矩阵 B^k 的正定性是一个很关键的因素，在 (2.4) 式两边同时左乘 (s^k)^T 可得 (s^k)^T B^{k+1} s^k=(s^k)^T y^k ，因此条件：

为 B^{k+1} 正定的一个必要条件。我们额外要求条件 (5) 在迭代过程中始终满足，这个条件也称为曲率条件。对于一般的目标函数 f(x) ，需要使用Wolfe准则线搜索来保证曲率条件 (2.5) 成立。实际上，根据Wolfe条件有 \nabla f(x^{k+1})^T s^k \geqslant c_2 \nabla f(x^k)^T s^k ，两边同时减去 \nabla f(x^k)^T s^k ，

这是因为 c_2 < 1 以及 s^k=\alpha_k d^k 是下降方向。仅仅使用Armijo准则不能保证曲率条件成立。

  在通常情况下，近似矩阵 B^{k+1} 或 H^{k+1} 是由上一步迭代加上一个修正得到的，并且要求满足割线方程 (2.3) 。如下是拟牛顿法的框架：

给定 x^0 \in R^n，初始矩阵 B^0 \in R^{n \times n}（或 H^0），令 k=0.

while未达到停机准则do

 计算方向 d^k=-(B^k)^{-1} \nabla f(x^k) 或 d^k=-H^k \nabla f(x^k).

 通过线搜索找到合适的步长 \alpha^k > 0，令 x^{k+1}=x^k+\alpha_kd^k。

 更新海瑟矩阵的近似矩阵 B^{k+1} 或其逆矩阵的近似 H^{k+1}.

 k \leftarrow k + 1

end while

3. 拟牛顿矩阵的更新方式

3.1 秩一更新（SR1）

  秩一更新（SR1）公式是结构最简单的拟牛顿矩阵更新公式。设 B^k 是第 k 步的近似海瑟矩阵，我们通过对 B^k 进行秩一修正得到 B^{k+1} ，使其满足割线方程 (2.3) 。为此使用待定系数法求出修正矩阵，并设：$$

B^{k+1}=B^k + auu^T \tag {3.1.1}

B^{k+1} s^k=(B^k + auu^T)s^k=y^k \tag {3.1.2}

(a \cdot u^T s^k)u=y^k - B^k s^k \tag {3.1.3}

a((y^k-B^ks^k)^T s^k)(y^k - B^ks^k) = y^k - B^k s^k \tag {3.1.4}

B^{k+1}=B^k+\dfrac{(y^k-B^ks^k)(y^k-B^ks^k)^T}{(y^k-B^ks^k)^Ts^k)} \tag {3.1.5}

H^{k+1}=H^k + \dfrac{(s^k-H^ky^k)(s^k-H^ky^k)^T}{(s^k-H^ky^k)^Ty^k} \tag {3.1.6}

B^{k+1}=B^k=auu^T+bvv^T \tag {3.2.1}

B^{k+1}s^k=(B^k+auu^T+bvv^T)s^k = y^k \tag {3.2.2}

(a \cdot u^Ts^k)u + (b \cdot v^Ts^k)v = y^k - B^ks^k \tag {3.2.3}

u=y^k, a \cdot u^T s^k=1, \ v = B^ks^k, b \cdot v^Ts^k = -1 \tag {3.2.4}

B^{k+1}=B^k+\dfrac{y^k(y^k)^T}{(s^k)^Ty^k} - \dfrac{B^ks^k(B^ks^k)^T}{(s^k)^TB^ks^k} \tag {3.2.5}

H^{k+1}=(I-\rho_k y^k(s^k)^T)^T H^k (I-\rho_k y^k (s^k)^T) + \rho_ks^k(s^k)^T \tag {3.2.6}

\min_H ||H-H^k||_W, \ s.t. H=H^T, \ Hy^k=s^k \tag {3.2.7}

||H||_W=||W^{1/2} HW^{1/2}||_F \tag {3.2.8}

H^{k+1}=H^k - \dfrac{H^ky^k(H^ky^k)^T}{(y^k)^TH^ky^k} + \dfrac{s^k(s^k)^T}{(y^k)^Ts^k} \tag {3.3.1}

B^{k+1}=(I-\rho_ks^k(y^k)^T)^TB^k(I-\rho_ks^k(y^k)^T)+\rho_ky^k(y^k)^T \tag {3.2.2}

\min_B ||B-B^k||_W, \ s.t. B=B^T, \ Bs^k=y^k \tag {3.2.3}

H^{k+1}=(V^k)^TH^kV^k + \rho_ks^k(s^k)^T

\rho_k=\dfrac{1}{(y^k)^Ts^k}, V^k=I-\rho_ky^k(s^k)^T \tag {3.4.2}

\begin{aligned}
H^k&=\nabla (V^{k-m}...V^{k-1})^TH^{k-m}(V^{k-m}...V^{k-1})+ \
&\rho_{k-m}(V^{k-m+1}...V^{k-1})+ \
&\rho_{k-m+1}(V^{k-m+2}...V^{k-1})^Ts^{k-m+1}(s^{k-m+1})^T(V^{k-m+2}...V^{k-1})+ \
&...+\rho_{k-1}s^{k-1}(s^{k-1})^T
\end{aligned} \tag {3.4.5}

V^{k-1} \nabla f(x^k), V^{k-2}V^{k-1}\nabla f(x^k),...,V^{k-m}...V^{k-2}V^{k-1} \nabla f(x^k) \tag {3.4.4}

\begin{aligned}
H^k \nabla f(x^k)&=(V^{k-m}...V^{k-1})^T H^{k-m}q+ \
&(V^{k-m+1}...V^{k-1})^Ts^{k-m}\alpha_{k-m}+ \
&(V^{k-m+2}...V^{k-1})^Ts^{k-m+1}\alpha_{k-m+1}+...+s^{k-1}\alpha_{k-1}.
\end{aligned} \tag {3.4.5}

\begin{aligned}
&(V^{k-m+1}...V^{k-1})^T((V^{k-m+1})^Tr+\alpha_{k-m}s^{k-m})\
&=(V^{k-m+1}...V^{k-1})^T(r+(\alpha_{k-m}-\beta)s^{k-m}).
\end{aligned} \tag {3.4.6}

\gamma_k = \dfrac{(s^{k-1})^T y^{k-1}}{(y^{k-1})^Ty^{k-1}} \tag {3.4.7}

S^k=[s^0,s^1, ...,s^{k-1}], Y^k=[y^0, y^1, ...,y^{k-1}] \tag {3.4.8}

f(x)=\sum_{i=1}^{n/2}[(-13+x_{2i-1}+((5-x_{2i})x_{2i}-2)x_{2i})^2+(-29+x_{2i-1}+((x_{2i}+1)x_{2i}-14)x_{2i})^2]

f(x)=\sum_{i=1}^{n-1}(x_i-1)^2+(\sum_{i=1}^{n}x_j^2-0.25)^2