【最优化(四)】无约束优化问题的牛顿法及其代码实现

BlaXuan

梯度法仅仅依赖于零阶信息（函数值）和一阶信息（梯度）来对目标函数做优化。如果目标函数充分光滑，则可以利用二阶导数的信息来构造下降方向。牛顿类算法就是利用二阶导数的信息来构造迭代方向的算法之一。由于利用到了更高阶的信息，牛顿法的实际效果一般要优于最速下降法，但牛顿法对目标函数的光滑性的要求更高，这也就决定了牛顿法的使用前提会更加地苛刻。本篇笔记介绍牛顿法和修正牛顿法两大经典算法，并且附上牛顿法，修正牛顿法和最速下降法的代码实现（方便大家教学使用），最后再给出牛顿法的收敛性证明。
1 牛顿法

无约束优化问题的表达式：
\underset{x\in R^n}{\min\text{ }}f\left( x \right) \quad (1.1) \\
在本节中我们要求 f(x) 是一个二次连续可微函数。
牛顿法的算法整体框架依然是三步走的框架，和上一节中最速下降法是同一个框架：

• 每一步走的方向（搜索方向）；
  
• 第二个是每一步走多远（搜索步长）；
  
• 第三个是何时停下来（最优性条件）；
  

我们将目标函数 f(x) 在第 k 步的迭代点处采用二阶泰勒展开：
f\left( x_k+d_k \right) =f\left( x_k \right) +\nabla f\left( x_k \right) ^Td_k+\frac{1}{2}d_{k}^{T}\nabla ^2f\left( x_k \right) d_k+o\left( \lVert d_k \rVert ^2 \right)  \\
上述二阶泰勒展开实际上是对目标函数 f(x) 在迭代点 x_k 初采用二次函数来近似。我们希望通过极小化上式可以选取合适的下降方向 d_k。如果忽略掉上式的高阶无穷小量，并将等式右端看成对于搜索方向 d_k 的极小化，由此可得：
\nabla ^2f\left( x_k \right) d_k=-\nabla f\left( x_k \right) \quad (1.2) \\
式(1.2)也被称为牛顿方程，我们假设 \nabla ^2f\left( x_k \right)是可逆的，由此可得：
d_k=-\nabla ^2f\left( x_k \right) ^{-1}\nabla f\left( x_k \right) \quad (1.3) \\
在得到上述搜索方向之后我们需要验证一下这个方向到底是不是一个下降方向？对于可微函数 f，在某一点 x_k\in R^n处，如果搜索方向满足(1.4)就是一个下降方向：
\nabla f\left( x_k \right) ^Td_k<0 \quad (1.4) \\
进一步我们将式(1.3)带入到(1.4)可得：
-\nabla f\left( x_k \right) ^T\nabla ^2f\left( x_k \right) ^{-1}\nabla f\left( x_k \right) <0 \quad (1.5) \\
显然 \nabla ^2f\left( x_k \right) ^{-1} 是正定矩阵的时候，式(1.5)可以满足。由此可知，若 \nabla ^2f\left( x_k \right) ^{-1} 正定则牛顿法的搜索方向 d_k 是一个下降方向。至此我们通过牛顿法可以得到迭代公式为：
x_{k+1}=x_k-\nabla ^2f\left( x_k \right) ^{-1}\nabla f\left( x_k \right) \quad (1.6) \\
在经典的牛顿法中我们一般把步长选取为 1，所以我们这里暂时就不考虑其它的步长选取方法，而最优性条件依然可以采用和之前最速下降法相同的条件。至此算法的关键三要素：搜索方向，搜索步长和最优性条件都已经确定了，下面我们直接给出经典的牛顿法的算法流程：
Step 1: 任意给出 x_0 \in R^n, \epsilon>0, k:=0
Step 2: 若满足终止条件 \lVert \nabla f\left( x_k \right) \rVert ^2<\epsilon，则停止迭代 输出x_k为最优解；若不满足终止条件则继续
Step 3: 令搜索方向 d_k=-\nabla ^2f\left( x_k \right) ^{-1}\nabla f\left( x_k \right)
Step 4: x_{k+1}:=x_k+d_k，k:=k+1，跳转到 Step 2
2 修正牛顿法

从上面牛顿法的推导过程中，我们不难发现牛顿法有三大缺点：

• 每步迭代都需要计算一个 n \times n 的逆矩阵 \nabla ^2f\left( x_k \right) ^{-1}，我们知道求逆矩阵的计算量是很大的，而原来的最速下降法就没有这个求逆的操作。
  
• 当 \nabla ^2f\left( x_k \right) ^{-1} 不是正定矩阵的时候，我们就无法保证牛顿法得到的搜索方向是一个下降方向，这个我们在式(1.5)中已经论证过了。
  
• 当初始的迭代点距离最优值较远的时候，选取步长为1会使得迭代极其不稳定，在某些时候会让算法发散。（关于这一缺点的具体论证我们放在牛顿法收敛性分析中去探讨，这里只需要暂时知道有这么一回事即可）。
  

为了克服上述缺陷，我们对经典的牛顿法进行了改进，主要是改进搜索方向和搜索步长。如下所示是修正牛顿法的算法流程：
Step 1: 任意给出 x_0 \in R^n, \epsilon>0, k:=0
Step 2: 若满足终止条件 \lVert \nabla f\left( x_k \right) \rVert ^2<\epsilon，则停止迭代 输出x_k为最优解；若不满足终止条件则继续
Step 3: 令搜索方向 d_k=-\left( \nabla ^2f\left( x_k \right) +\mu I \right) ^{-1}\nabla f\left( x_k \right)，其中 \mu>0, I 是单位矩阵。
Step 4: 由线搜索方法确定出步长 \alpha_k
Step 5: x_{k+1}:=x_k+\alpha_kd_k，k:=k+1，跳转到 Step 2
修正牛顿法和经典的牛顿法的主要区别是 Step 3 中搜索方向 多加了一个 单位矩阵 \mu I ，另外就是步长的计算是通过线搜索方法来确定的而不是像之前直接就是取1，这两点改进可以解决我们前面所述的 缺点2和3，关于缺点1需要计算 Hessian 的逆矩阵的问题我们将会在下一节 拟牛顿法去解决这个问题。
3 代码实现

代码的主干架构主要由三部分构成：
main.py：代码的主干结构
model.py：定义优化问题的决策变量，目标函数，目标函数一阶和二阶导数等信息
optimizer.py：定义各种优化方法，包括梯度法，牛顿法，修正牛顿法和拟牛顿法等。
以下仅展示部分代码，完整代码可见：https://github.com/WenYuZhi/DifferentiableOptimizer/tree/master/unconstraint
from model import Model
from optimizer import Gradient
from optimizer import Netwon
from optimizer import ModifiedNetwon
import numpy as np
import pandas as pd

dim = 10 # 决策变量维数
model = Model(dim = dim) # 优化模型
model.set_objective(name = 'quad_fun1') # 设置目标函数
model.diff() # 求一阶导数
model.diff_second() # 求二阶导数

n_iter, eps = 100, 10**(-6) # 迭代次数和收敛误差
optimizer = ModifiedNetwon(model, eps) # 指定修正牛顿法
optimizer.set_init(lb=0, ub=100) # 生成初始解

for k in range(n_iter):
    optimizer.search_direction(miu = 0.1) # 搜索方向
    optimizer.get_step() # 搜索步长
    optimizer.update()   
    optimizer.moniter(iter_time = k)
    if optimizer.is_stop():
        break

# optimizer.save()


class Optimizer:
    def __init__(self, model, eps) -> None:
        self.m, self.eps = model, eps
        self.trace_obj, self.trace_x = [], []
        self.trace_step, self.trace_d = [], []
   
    def set_init(self, lb = 0, ub = 1)-> None:
        assert(lb <= ub), print("ub is larger than lb")
        self.x = (ub - lb) + np.random.rand(self.m.dim) + lb * np.ones(self.m.dim)
   
    def first_order(self):
        self.grad, self.obj_values = self.m.diff_eval(self.x), self.m.obj_fun_eval(self.x)
        return self.grad, self.obj_values
   
    def second_order(self):
        self.hessian = self.m.diff_second_eval(self.x)
        return self.hessian
   
    def get_step(self):
        # self.step = 0.01
        self.step = self.__armijo_step() # 默认采用 Armijo 方法获取步长
   
    def __armijo_step(self):
        alpha, beta, pho = 10, 0.1, 0.8
        for _ in range(100):
            if self.m.obj_fun_eval(self.x + alpha * self.d) <= self.obj_values + alpha * pho * np.dot(self.grad.T, self.d):
                return alpha
            alpha *= beta
        return False

    def __exact_step(self):
        pass

    def moniter(self, iter_time): # 显示迭代过程中的结果
        self.trace_obj.append(self.obj_values)
        self.trace_x.append(self.x)
        self.trace_d.append(self.d)
        self.trace_step.append(self.step)
        if iter_time % 10 == 9:
            print("iterations {} - objective function: {}".format(iter_time, self.m.obj_fun_eval(self.x)))
   
    def is_stop(self): # 判断收敛准则
        if np.linalg.norm(self.grad, 1.0) <= self.eps:
            return True
        return False
   
    def update(self):
        assert(self.step != False), print("objective function is not decreased")
        self.x += self.step * self.d
   
    def save(self):
        self.file_route = "./results//"
        pd.DataFrame(self.trace_obj).to_csv(self.file_route + "trace_obj.csv")
        pd.DataFrame(self.trace_x).to_csv(self.file_route + "trace_x.csv")
        pd.DataFrame(self.trace_d).to_csv(self.file_route + "trace_d.csv")
        pd.DataFrame(self.trace_step).to_csv(self.file_route + "trace_step.csv")

class Netwon(Optimizer):
    def __init__(self, model, eps) -> None:
        super().__init__(model, eps)
   
    def search_direction(self):
        self.grad, self.obj_values = self.first_order()
        self.hessian = self.second_order()
        self.d = -1.0 * np.dot(np.linalg.inv(self.hessian), self.grad)

class ModifiedNetwon(Optimizer):
    def __init__(self, model, eps) -> None:
        super().__init__(model, eps)
   
    def search_direction(self, miu):
        self.grad, self.obj_values = self.first_order()
        self.hessian = self.second_order()
        self.d = -1.0 * np.dot(np.linalg.inv(self.hessian + miu*np.eye(self.hessian.shape[0])), self.grad)

class Gradient(Optimizer):
    def __init__(self, model, eps) -> None:
        super().__init__(model, eps)
   
    def search_direction(self):
        self.grad, self.obj_values = self.first_order()
        self.d = -1.0 * self.grad

4 牛顿法收敛性分析

最后我们对牛顿法的收敛性进行一个分析。首先我们假设

• 目标函数 f(x) 是二阶连续可微的函数；
  
• 二阶导函数在最优值点 x^* 的邻域内 N_{\delta}\left( x^* \right) 是李普希兹连续的，即存在常数 L>0 使得：
  

\lVert \nabla ^2f\left( x \right) -\nabla ^2f\left( y \right) \rVert \le L\lVert x-y \rVert ,\ \forall x,y\in N_{\delta}\left( x^* \right) \quad (4.1) \\
如果目标函数 f(x) 在点 x^* 处满足 \nabla f\left( x^* \right) =0,\nabla ^2f\left( x^* \right) \succ 0，则对于牛顿法有如下结论：

• 如果初始点离 x^* 足够接近，则牛顿法产生的迭代点序列 \{x_k\} 收敛到 x^*；
  
• 序列 \{x_k\} 收敛到 x^* 的速度是 Q 二次的；
  
• 梯度的范数 \lVert \nabla f_k \lVert 序列收敛到0的速度是二次的；
  

证明：根据牛顿法的迭代公式(1.6)可得
x_{k+1}-x^*=x_k-\nabla ^2f\left( x_k \right) ^{-1}\nabla f\left( x_k \right) -x^* \quad (4.2) \\
由于 \nabla f\left( x^* \right) =0 进一步整理可得：
x_{k+1}-x^*=\nabla ^2f\left( x_k \right) ^{-1}\left[ \nabla ^2f\left( x_k \right) \left( x_k-x^* \right) -\left( \nabla f\left( x_k \right) -\nabla f\left( x^* \right) \right) \right] \quad (4.3) \\
根据泰勒公式有
\nabla f\left( x_k \right) -\nabla f\left( x^* \right) =\int_0^1{\nabla ^2f\left( x_k+t\left( x^*-x_k \right) \right) \left( x^*-x_k \right) dt} \quad (4.4) \\
因此我们有
\lVert \nabla ^2f\left( x_k \right) \left( x_k-x^* \right) -\left( \nabla f\left( x_k \right) -\nabla f\left( x^* \right) \right) \rVert \quad (4.5) \\=\lVert \int_0^1{\left[ \nabla ^2f\left( x_k+t\left( x^*-x_k \right) \right) -\nabla ^2f\left( x_k \right) \right] \left( x_k-x^* \right) dt} \rVert  \quad (4.6) \\\le \int_0^1{\lVert \nabla ^2f\left( x_k+t\left( x^*-x_k \right) \right) -\nabla ^2f\left( x_k \right) \rVert \lVert x_k-x^* \rVert dt} \quad (4.7) \\\le \lVert x_k-x^* \rVert ^2\int_0^1{Ltdt} \quad (4.8) \\=\frac{L}{2}\lVert x_k-x^* \rVert ^2 \quad (4.9) \\
从(4.5)到(4.6)是因为(4.4)满足，从(4.6)到(4.7)是利用积分的性质，从(4.7)到(4.8)是因为假设条件2中的李普希兹连续，从(4.8)到(4.9)是显然的。

IT圈老男孩1

第4节收敛性分析中关于牛顿法收敛性的三个结论中，1和3好像重复了。

JamesB

追本溯源，牛顿法一开始用来求解f（x）=0的根，当时的主流是几何图形分析。对于优化问题，等价于用牛顿法求解grad（f）=0，有可能会收敛到驻点。

BlaXuan

谢谢提醒已更正。