工程优化：牛顿法的缺陷及拟牛顿法的概念与实现

九月 21, 2019

当前设备屏幕尺寸过小，推荐使用PC模式浏览。

引用

Wikipedia:Quasi-Newton method

Preliminaries

跟上一篇一样，要读懂这些内容，需要掌握以下内容：

梯度

对于一维函数 $f(x)$ ，其导数定义为：

f'(x)=\lim \limits_{\Delta x \rightarrow 0} \frac{f(x_0+\small{\Delta} x)-f(x_0)}{\small{\Delta} x}

对于多维函数 $f(x_1,...,x_n)$ ，对 $x_i$ 求导数 $\frac{df}{dx_i}$ ，将其记为偏导数 $\frac{\partial f}{\partial x_I}$ 。特别的，记录梯度 $\triangledown f(x)$ 或简记为 $\triangledown f$ 为对 $x_i$ 求偏导后的列向量：

\triangledown f(x)=(\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2},..., \frac{\partial f}{\partial x_1})^T

海森矩阵(Hessian matrix)

若存在 $f:\R^n \rightarrow \R$ ，即一个多维输入 $x\in\R$ 到一维输出的映射，若其在任意维度上都二阶可导，则定义其海森矩阵：

H = \begin{bmatrix} \frac{\partial f^2}{\partial x_1^2} & \frac{\partial f^2}{\partial x_1\partial x_2} & \cdots &\frac{\partial f^2}{\partial x_1 \partial x_n} \\ \frac{\partial f^2}{\partial x_2 \partial x_1} & \frac{\partial f^2}{\partial x_2^2} & \cdots &\frac{\partial f^2}{\partial x_2 \partial x_n} \\ \vdots & \vdots & \ddots& \vdots \\ \frac{\partial f^2}{\partial x_n \partial x_1} & \frac{\partial f^2}{\partial x_n\partial x_2} & \cdots &\frac{\partial f^2}{ \partial x_n^2} \end{bmatrix}

显然， $H^T = H$ 、 $H$ 的尺寸为 $n\times n$ 。

雅可比矩阵(Jacobian matrix)

若存在 $f:\R^n \rightarrow \R^m$ ，即一个多维输入 $x\in\R^n$ 到多维输出 $f(x)\in\R^m$ 的映射，则 $f$ 的雅可比矩阵：

\begin{aligned} J &= [ \frac{\partial f}{\partial x_1} \frac{\partial f}{\partial x_2} \cdots \frac{\partial f}{\partial x_n}] \\ &= \begin{bmatrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \cdots & \frac{\partial f_m}{\partial x_n} \\ \end{bmatrix} \end{aligned}

显然，雅可比矩阵的尺寸为 $n \times m$ , \ $J_{ij} = \frac{\partial f_i}{\partial x_j}$ 。

泰勒定理(Taylor’s theorem)

Wikipedia:Taylor’s theorem
Tips：泰勒公式的推导和理解>见此<。

给出泰勒公式：

f(x) = f(a)+f'(a)(x-a)+\cdots+\frac{f^{(n)}(x)}{n!}(x-a)^n+R_n

$R_n$ 在此处被称为皮亚诺余项(Peano form of the remainder)。经由拉格朗日推导和柯西推导，还有拉格朗日余项和柯西余项等，但它们大部分都是等价的。

牛顿法(Newton’s method)

>上一篇<中提到了标准形式的牛顿法，并利用代码给予了实现，此处不再赘述。

正定矩阵

设有 $x^TAx$ 对 $\forall x=(x_1,\cdots,x_n)^T$ 都有 $x^TAx>0$ ，则称矩阵 $A$ 是正定矩阵；若都有 $x^TAx\ge0$ ，则称其为半正定矩阵。同时，若一个矩阵 $A$ 正定，则 $A$ 的特征值均为正数(半正定则为大于零的数)，一定存在逆矩阵 $A^{-1}$ 。

Sherman–Morrison公式

Wikipedia:Sherman–Morrison formula

设有可逆矩阵 $A \in \R^{n \times n}$ ，列向量 $u,v \in \R^n$ ，若 $A+uv^T$ 可逆，当且仅当 $1+v^TA^{-1}u =\not 0$ ：

(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}

>返回目录

拟牛顿法(Newton’s method)

牛顿法两个层面

牛顿法：迭代求根——牛顿-拉弗森求根法。
牛顿法：求取极值——二阶拟合优化。

优化问题中牛顿法的缺陷

在>上一篇<中，提到牛顿法是从泰勒公式的二阶展开形式得到：

f(x)\approx f(x_0)+(x-x_0)^Tg_{(x_0)}+\frac{1}{2}(x-x_0)^TH_{(x_0)}(x-x_0)

对其求导数：

\begin{aligned} f'(x) &\approx g_{(x_0)}+\frac{1}{2}(H_{(x_0)}+H_{(x_0)}^T)(x-x_0)\\ &\approx g_{(x_0)}+H_{(x_0)}(x-x_0) \end{aligned}

令其等于0即可得到下降根：

x_1 = x_0 - H_{(x_0)}^{-1}g_{(x_0)}

即牛顿下降公式： $x_{t+1} = x_t - H_{(x_t)}^{-1}g_{(x_t)}$ 。

这里， $H_{(x_t)}^{-1}g_{(x_t)}$ 被称为牛顿方向。
Tips:牛顿法实际上是利用泰勒公式拟合函数在 $x_0$ 处的曲线，并利用函数的驻点处导数必为0这一必要条件进行迭代的，这一做法被称为拟牛顿条件，将在下文中得到大量应用。

但是，很容易注意到牛顿法在目标函数的表面非凸的时候，海森矩阵的特征值并非全正（如当前迭代点为一个鞍点），牛顿方向实际指向了一个极大值而非我们想要得到的极小值。另外，海森矩阵的计算较为复杂，其正定条件也较为苛刻，不是随时都能满足。

总结一下它的缺点：

求海森矩阵较为复杂。(为此引入拟牛顿法)

迭代点曲率过大时有可能不收敛。(为此引入阻尼牛顿法)

海森矩阵非正定，即其可能不存在逆矩阵。(为此引入正则化方法)

牛顿法的改进

阻尼牛顿法

牛顿法仅仅是对函数的局部进行一次拟合，并不能完全说明函数的走向。因此，完全有可能出现局部曲率极大的迭代点，使得一次牛顿法下降得到的点出现大值跳变。(例如，对于 $f(x)=x^2+sin(\pi x)$ 使用牛顿-拉弗森求根法，搜寻其负根时，会发现从10和-10处下降完全不会回归到各自应该下降到的一个根，甚至完全不收敛。

因此，引入梯度下降中步长的概念，即阻尼牛顿下降公式： $x_{t+1} = x_t - \lambda H_{(x_t)}^{-1}g_{(x_t)}$ 。

此时，阻尼牛顿法一定程度上消减了牛顿法的优势，可能需要借助 Momentum 等方法以加速收敛速度。同时，该算法每次迭代仍旧需要求取一次海森矩阵的逆，然而通常情况下并不能保证其一定存在二阶导数和正定的海森矩阵。

这就极大的限制了运用范围，因此，又一改进算法被提了出来。

正则化海森矩阵

Wikipedia:Levenberg–Marquardt algorithm

针对这个问题，可以对海森矩阵引入“惩罚”，即加入正则化形式( $I$ 为一个单位矩阵)：

H=H+\alpha I

Levenberg-Marquardt 算法就是这个思想的优秀执行者，实验证明只要海森矩阵的负特征值相对接近0，效果就会很好。通常来说，当曲线的负曲率的绝对值越大，正则项超参数 $\alpha$ 就越大，才足以修正牛顿方向。

但是，需要注意的是，当 $\alpha$ 足够大时，牛顿法实际上已经退化成了梯度下降法。

因此，牛顿法依旧还有可以改进的空间，在介绍真正的拟牛顿法前，先要介绍一下拟牛顿条件。

拟牛顿条件

展开泰勒格式的二阶形式并求其驻点，即使其导数为0，在 $t$ 时刻得到:

g_{(t+1)} = f'(x_{t+1}) \approx g_{(x_t)}+H_{(x_0)}(x-x_t)

记 $\Delta x_{(t)}=(x_{t+1}-x_t)$ 、 $\Delta g_{(t)}=g(x_{t+1})-g(x_t)$ ：

\begin{aligned} \Delta g_{(t)} &\approx H_{(x_t)} \Delta x_{(t)} \\ \Delta x_{(t)} &\approx H_{(x_t)}^{-1} \Delta g_{(t)} \end{aligned}

这个公式就是拟牛顿条件，也就是搜索函数驻点的必要条件。只要满足这个条件，则牛顿法成立。

注：这个方程也被称为割线方程（Secant Equation）。

由于上文分析过函数的二阶导即海森矩阵有种种弊端，那么有没有办法得到一个替代矩阵 $B$ ，替代海森矩阵的逆即 $H_{(x_t)}^{-1}$ ，使得不用求导海森矩阵。这个矩阵被称为准海森矩阵，显然准海森矩阵应当满足拟牛顿条件。

割线方程也被记为： $\gamma_t=B_t^{-1}y_t$ ，其中： $t = 0,1\cdots n$ 、 $\gamma_t = \Delta x_{(t)}$ 、 $y_t=\Delta g_{(t)}$ 。

Wiki中选择用 $B$ 来替代 $H$ 而不是 $H^{-1}$ 。

另外，由于牛顿法是不断迭代的，我能求得某一时刻 $t$ 的拟牛顿条件及其替代矩阵 $B_t$ ，能否浮现迭代过程求得下一时刻的替代矩阵 $B_{t+1}$ ？

总结一下，拟牛顿法的思路是：
寻找满足割线方程的矩阵 $B$ 的求法，若：

$B_t$ 是正定的。

$\Delta B$ 简单可求。

因此， $B_{t+1}=B_t+\Delta B$ 简单可求且正定。

拟牛顿法：DFP方法的实现

Wikipedia:Davidon–Fletcher–Powell formula

说起来有点Tricky，DFP方法求得 $B$ 基于“待定法”，有一定的数学技巧因素。

假设有一个初始准海森矩阵 $B_0=I$ ，尝试求得 $\Delta B$ 的计算方法，以此推得每一时刻的割线方程。
DFP法基于数学上的Trick，构造了这样一个矩阵：

\Delta B_t=\alpha uu^T+ \beta vv^T

其中 $u,v \in \R^n$ 是两个待定列向量， $\alpha$ 和 $\beta$ 是两个待定标量，显然其为对称矩阵。尝试求下面的公式：

B_{t+1}=B_t+\Delta B_t

代入割线方程：

\begin{aligned} \gamma_t &= (B_t+\Delta B_t)y_t \\ &= B_ty_t+\alpha uu^Ty_t + \beta vv^Ty_t \\ &= B_ty_t+(\alpha u^Ty_t)u + (\beta v^Ty_t)v \end{aligned}

$\alpha u^Ty_t$ 和 $\beta v^Ty_t$ 为两个标量，DFP尝试找到尽量正交的两个基，因此分别令其等于1和-1，即：

\alpha = \frac{1}{u^Ty_t} , \ \beta = -\frac{1}{ v^Ty_t}

代入上式：

\gamma_t - B_ty_t = u-v

由于我们是找两个“正交基”，显然 $u=\gamma_t$ 和 $v=B_ty_t$ 是上式的一个解。
至此通过待定法，就将 $\Delta B_t$ 解出来了( $B_0=I$ )：

\begin{aligned} \Delta B_t &= \alpha uu^T+ \beta vv^T \\ &= \frac{\gamma_t\gamma_t^T}{\gamma_t^Ty_t} - \frac{B_ty_ty_t^TB_t^T}{y_t^TB_t^Ty_t} \\ &\mathop = \limits_{B_0=B_0^T} \frac{\gamma_t\gamma_t^T}{\gamma_t^Ty_t} - \frac{B_ty_ty_t^TB_t}{y_t^TB_ty_t} \\ B_{t+1}&=B_t + \frac{\gamma_t\gamma_t^T}{\gamma_t^Ty_t} - \frac{B_ty_ty_t^TB_t}{y_t^TB_ty_t} \end{aligned}

Tip：如果 $B_0$ 是一个单位矩阵（或对称矩阵），则从此往后，根据此法迭代出来的每个 $B_t$ 矩阵显然是主对角线对称矩阵，即： $B_t=B_t^T$ 。

简单对 $cost(X) = \sum_i^N(x_{0(i)}^2+\sin(\pi x_{1(i)}))$ 进行牛顿法下降和DFP法下降，部分代码见下，完整代码>见此<或>点此下载代码<：

def DFP(x_init, epsilon=1e-2, m_lambda=0.01):
    start_time = time.time()
    costs = []
    xt = x_init
    gt = grad(xt)
    Bt = np.array([[1, 0], [0, 1]]).reshape((2, 2))  # B_0
    while True:
        delta_x0, delta_x1 = gt
        if abs(delta_x0) < epsilon and abs(delta_x1) < epsilon:
            break

        xt1 = xt.T - m_lambda * Bt.dot(gt)  # 为了使其收敛，步长放小为0.01
        xt1 = xt1.T  # 得到新的迭代点

        gamma_t = (xt1 - xt).T
        gt1 = grad(xt1)  # 得到新的梯度
        yt = gt1 - gt

        # 向量化计算Bt时，有点计算技巧，即先算标量，后点乘求和
        # (gamma_t / gamma_t.T.dot(yt).T).dot(gamma_t.T)
        Bt = Bt + (gamma_t / gamma_t.T.dot(yt).T).dot(gamma_t.T) - \
             Bt.dot(yt).dot(yt.T).dot(Bt) / yt.T.dot(Bt).dot(yt)  # 计算新的替代矩阵
        gt = gt1
        xt = xt1

        current_cost = cost(xt)
        costs.append(current_cost)
    print("DFP: xt^* =\n", xt)
    print("DFP: time usage =", time.time() - start_time)
    print("DFP: iter_times =", len(costs))
    return costs, xt

效果如图：

同时，会发现在运行时间上，DFP实现的拟牛顿法要比牛顿法快上不少。

严格来说，DFP的推理及得到，是基于共轭梯度法的思想的，此处不再赘述。

拟牛顿法：BFGS方法的实现

Wikipedia:Broyden–Fletcher–Goldfarb–Shanno algorithm

DFP法尝试逼近的是海森矩阵的逆，而BFGS法的核心要义是直接逼近海森矩阵 $H_{(x_t)}$ ，即有：

\begin{aligned} \Delta g_{(t)} &\approx H_{(x_t)} \Delta x_{(t)} \\ y_t &\approx B_{(x_t)} \gamma_t \end{aligned}

$\gamma_t = \Delta x_{(t)}$ 、 $y_t=\Delta g_{(t)}$ 。

利用同样的数学Trick和待定法，得到：

\alpha=\frac{1}{u^T\gamma_t},\ \beta = -\frac{1}{ v^T\gamma_t}

y_t - B_t\gamma_t = u-v

u=y_t,\ v=B_t\gamma_t

以上三式联立，设 $B_0=I$ ，代入割线方程即可得到BFGS法的替代矩阵：

B_{t+1}=B_t+\frac{y_ty_t^T}{y_t^T\gamma_t}-\frac{B_t\gamma_t\gamma_t^TB_t}{\gamma_t^TB_t\gamma_t}

此时虽然得到了对海森矩阵直接进行逼近的替代矩阵 $B$ ，但是在搜索时还要求一次逆就非常占用时间，引入Sherman–Morrison公式，即可得到：

B_{t+1}^{-1}=(I-\frac{\gamma_ty_t^T}{y_t^T\gamma_t})B_t^{-1}(I-\frac{y_t\gamma_t^T}{y_t^T\gamma_t})+\frac{\gamma_t\gamma_t^T}{y_t^T\gamma_t}

这个公式的推导过程较为冗长，大致的思路是，基于Sherman–Morrison公式：

(A+uv^T)^{-1}=A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}

令 $A:=B_t+\frac{y_ty_t^T}{y_t^T\gamma_t}$ 、 $u:=-\frac{B_t\gamma_t}{\gamma_t^TB_t\gamma_t}$ 、 $v:=B_t\gamma_t$ 进行第一次代入；然后对 $B_t+\frac{y_ty_t^T}{y_t^T\gamma_t}$ ，令 $A:=B_t$ 、 $u:=\frac{y_t}{y_t^T\gamma_t}$ 、 $v:=y_t$ 进行第二次代入，最后进行化简。