数值优化（5）——信赖域子问题的求解，牛顿法及其拓展

学弱猹

发布于 2021-08-09 16:41:13

1.6K0

发布于 2021-08-09 16:41:13

文章被收录于专栏：学弱猹的精品小屋

上一节笔记：数值优化（4）——非线性共轭梯度法，信赖域法

————————————————————————————————————

大家好！俗话说得好，DDL是唯一生产力……在DDL的逼迫下，高产自然就来了2333。

在上一节最后我们介绍了信赖域方法，但是遗留了一个问题：二次模型怎么解？毕竟如果柯西点就足够的话，似乎也有点太没含金量了点。所以这一节我们会介绍解二次模型的具体的算法。

那么我们开始吧

Source

Nocedal, Wright, Numerical Optimization
厦门大学课堂笔记，教授主页：https://www.math.fsu.edu/~whuang2/index.html
厦门大学课堂笔记，教授主页：http://math.xmu.edu.cn/group/nona/kuidu.html
P.E. Gill, W. Murray, and Wright M. H. Practical Optimization

信赖域方法的子问题求解

上一节我们留下了二次问题怎么解的问题。事实上，我们有如下的这个性质

Proposition 1:

p^*

是信赖域子问题

\min_{\|p\| \le \Delta}f(x) + \nabla f(x)^T p + \frac12 p^T B p

的解当且仅当存在

\lambda \ge 0

，使得

(B + \lambda I)p^* = -\nabla f(x)

，且

\lambda (\Delta - \|p^*\|) = 0

，并且

B + \lambda I

为对称半正定的矩阵。

求解可以采用Lagrange乘数法和KKT条件的思路。你可以通过KKT条件得到前两个结论，然后通过反证法得到

B+ \lambda I

为一个对称半正定矩阵。当然了其实你可以通过对信赖域框架的理解，也即

p^*

是一个下降方向，来得到这个结论。不过这不是这一节的重点，所以我们不在这里给出这个性质的详细证明。

事实上，虽然这样的子问题是存在一个解析解的形式的，但是有一个问题是我们需要依赖一个线性方程组

(B + \lambda I)p^* = -\nabla f(x)

，我们需要求解这个

\lambda

才能够得到我们的解，如果对于一个规模很大的问题来说，这样的问题的求解就不再那么友好了。因此我们这里再介绍一个可以用来求解信赖域子问题的近似解的方法。

逼近信赖域子问题的狗腿法

下面这张图给出了一个狗腿法（dogleg method）的图示

其中

g_k

表示梯度。这里的

y = x^k + d^k

表示的是在当前迭代点上（注意这里的

x^k

是指第

步迭代点的含义，书上的记号和这一节记号不太一样，注意区分）对应的信赖域子问题的极小值点。这里的

表示的是信赖域子问题上线搜索取精确步长得到的解，

是指无约束的情况下信赖域子问题的最优解，那么这样的话，连接

，就可以得到一个和信赖域相交的点

，这个点就是所谓的狗腿点，这个方法就是狗腿法。可以看出在这个情况下，狗腿点算是一个真实的极小值点的还算不错的近似。

下图给出了一条狗的图示，如果有人养过狗或者被狗追过的话，可以看出狗腿和

x^kuv

这一条折线的形状类似，可能这也就是狗腿法得名的原因吧。

如果推广一下，其实狗腿法的迭代公式可以写成这么一个形式

\tilde p(\tau) = \begin{cases}\tau p^U & 0 \le \tau \le 1 \\ p^U + (\tau - 1) (p^B - p^U) & 1 \le \tau \le 2\end{cases}

这里的

p^U = - \frac{g_k^Tg_k}{g_k^TB_kg_k}g_k

，

p^B = -B_k^{-1}g_k

，至于这个

\tau

怎么选取，就会取决于我们的信赖域半径

\Delta

。事实上，可以说明的事情是

\|\tilde p(\tau)\|

是严格递增的，子问题也是严格递减的，但是同样的，因为不是重点，我们也不再这里给出证明。

牛顿法

牛顿法（Newton Method）也是一个很经典的迭代算法，它的思路非常简单：我们直接找方程

\nabla f(x) = 0

的根。通过泰勒展开，容易得到

\nabla f(x + p) \simeq \nabla f(x) + \nabla^2 f(x) p

如果我们希望迭代收敛，自然也就希望

\nabla f(x + p) \simeq 0

，也就是说可以求得

p = -(\nabla ^2f(x))^{-1}\nabla f(x)

，所以不难得到我们的迭代公式

x_{k + 1} = x_k - (\nabla^2 f(x_k))^{-1} \nabla f(x_k)

注意这个迭代公式是没有步长选取的过程的。一个很有意思的结论是它在局部上具有二次收敛性。

Theorem 1: 设

x^*

为

的最小值，设

f \in C^2, \nabla f(x^*) = 0, \nabla ^2 f(x)

正定，并且在

x^*

的一个邻域上Lipschitz连续，那么如果说

x_0

与

x^*

足够靠近，则牛顿法会使得

x_0

收敛于

x^*

，且具有二次收敛速度。

我们证明一下这个结论。既然要看收敛速度，那么自然要计算迭代点与最小值点的差值。设当前的迭代点为

，那么注意到

\|x+p - x^*\| = (\nabla^2 f(x))^{-1}[\nabla ^2 f(x) ( x - x^*) - (\nabla f(x) - \nabla f(x^*))]

我们希望利用到海塞矩阵的信息，这样就可以利用连续性做放缩。所以我们考虑Taylor展开，得到

\nabla f(x) - \nabla f(x^*) = \int_0^1 \nabla f^2(x + t(x^* - x))(x - x^*)dt

只需要再额外注意到

\nabla^2 f(x) (x - x^*) = \int_0^1 \nabla ^2f(x)(x - x^*)dt

，就可以得到我们的式子为

\|\nabla ^2 f(x)\|\|x + p -x^*\| = \|\int_0^1 [\nabla ^2 f(x) - \nabla ^2 f(x + t(x^* - x))](x - x^*) dt\|\\ \le \int_0^1 \|\nabla ^2 f(x) - \nabla ^2 f(x + t(x^* - x))\|dt\|x - x^*\| \le \frac 1 2L\|x - x^*\|^2

现在我们要利用一个性质，注意到

\nabla ^2 f(x^*)

是正定的，也就是说在这个极小值的局部上，它应该是凸的。那么如果足够的靠近，比方说设

\|x - x^*\| \le \gamma

，那么也应该具有一定的凸性。换句话说，存在

\gamma

，使得对于任意的

\|x - x^*\| \le \gamma

，我们有

\|\nabla ^2 f(x)\|^{-1} \le 2 \|\nabla ^2 f(x^*)\|^{-1}

把这个不等式代进去，就可以得到有

\|x + p - x^*\| \le L \|\nabla ^2 f(x^*)\|^{-1}\|x - x^*\|^2

所以如果我们证明了收敛，也就可以通过这个式子得到二次收敛的速度。而这只需要通过拆解，看出

\|x + p - x^*\| \le L \|\nabla ^2 f(x^*)\|^{-1}\|x - x^*\|\|x -x^*\|

所以只要令

\|x _0 - x^*\| \le \frac1 {2LL \|\nabla ^2 f(x^*)\|^{-1}}

，根据归纳法直接就可以得到

\|x_{k+1} - x^*\| \le \frac12 \|x_k -x^*\|

，也就说明了收敛性。到此也就算是说清楚了这个证明。

事实上通过这个性质，不难看出牛顿法的快速收敛性，但是存在的问题也很多。如果我们的迭代点并不靠近最小值，怎么保证收敛性呢？解方程的时候，如何保证搜索方向是下降方向呢？事实上，单纯的牛顿法并无法解决这些问题，这就引出了下面的一些修正的方法。

非精确牛顿法

非精确牛顿法（Inexact Newton Method）是在牛顿法的基础上，针对它无法解决的那些问题进行修正得到的方法。其本质上来说，就是解下面的这个方程组

(\nabla ^2 f(x) + E_i)p = - \nabla f(x)

归根到底，我们都是希望使得修正后的矩阵为对称正定阵，这样才能保证方程组有解，并且搜索方向为下降方向。

正如岭回归的想法一样，我们可以考虑通过设置

E_i = kI, k > 0

，使得我们的矩阵的特征值都大于0。活着可以考虑通过特征值分解的方式，将所有的非正的特征值都“强制修改为”某一个正数。这是可以通过矩阵的谱分解得到的，也就是说可以修改矩阵为

A = \sum_{i = 1}^n \max(\lambda_i, \delta) u_i u_i^T

这里的

\lambda_i, u_i

就是对应的

\nabla^2 f(x)

的特征值与特征向量。但是事实上谱分解是一个很耗时的方法，实际情况下会有很多麻烦。即使是考虑直接加上

的这样的方式，也需要估计最小特征值。因此这样的方法一般都不建议采纳。

有没有更好的方法呢？事实上对于矩阵分解而言，想要知道它的特征值，最好的方式就是Cholesky分解。具体来说，就是考虑下面的这个分解

\begin{bmatrix}a_{11} & a_{21} & a_{31} \\ a_{21} & a_{22} & a_{32} \\ a_{31} & a_{32} & a_{33}\end{bmatrix} = \begin{bmatrix}1 & & \\ l_{21} & 1 & \\l_{31} & l_{32} & 1\end{bmatrix}\begin{bmatrix}d_1 & & \\ & d_2 & \\& & d_3\end{bmatrix}\begin{bmatrix}1 & l_{21} & l_{31}\\ & 1 & l_{32}\\ & & 1\end{bmatrix}

简写为

A = L D L^T

。这个分解中的

d_i

就是特征值。所以如果说这个矩阵不是对称正定的，自然

d_i

就需要做修改。在提修改的思路之前，我们先把Cholesky分解的算法实现放在下面。

我们做如下的调整。

d_j = c_{jj}

\to

d_j = \max(|c_{jj}|, (\frac{\max_{j < i \le n}(|c_{ij}|)}{\beta})^2, \delta)

这里的

\delta, \beta

均提前给定。所以事实上你会发现修改之后，有如下的结论

Proposition 2:

d_j \ge \delta

，并且

|l_{ij}\sqrt{d_j}| \le \beta

第二个小结论直接根据算法中的

l_{ij} = c_{ij}/d_j

即可得到。这个结论是为了保证分解得到的新矩阵的元素不会显得太过于浮夸，这样一定程度上解决了不稳定的问题。当然了，这样的修改也确实使得方法具备全局收敛性了，分解的矩阵的条件数也有了上界，Source中有提供参考的文献，这里我们就不列证明了。

最后提一个稀疏矩阵的处理技巧。如果说海塞矩阵比较的稀疏的话，可以考虑做一个行列变换，使得它非零元素的部分都堆积在右下角，这样的话会使得Cholesky分解得到的矩阵

依然保持稀疏性。这样会给求解这个方程组带来方便，因为

LDL^T y = -\nabla f(x)

这个方程组的求解非常简单，因为

L, D

的结构均比较特殊，数值分析中有介绍专门的方法，这里就不细说了。

牛顿CG方法

牛顿CG方法是一种截断共轭梯度（Truncated Conjugate Gradient）方法，它的思想是通过带截断的共轭梯度法来求解牛顿法中涉及到的方程组。为什么说是“带截断的”呢？是因为我们在上一节提过，线性共轭梯度法的计算复杂度主要在于

，

p^TAp

等这样的东西，如果对应的矩阵

是一个稠密的矩阵，这个复杂度是

2n^3

（如果我们考虑CG算法做

步迭代的话，这是因为对于一个线性系统而言，算法最多

步就会收敛）。这是比Cholesky分解要慢的（Cholesky分解是

\frac13 n^3

）。所以实质上，我们会加一个限制使得算法能够具有“提前终止”的特性，这样一方面保证了较低的计算复杂度，另一方面也可以提高收敛速度。

我们先把带截断的牛顿CG方法的算法放在下面。

我们一步一步来解释这个方法

Bp = -g

对应的就是最开始的方程组

\nabla ^2 f(x) p = -\nabla f(x)

。因为我们的搜索方向就需要满足这个方程组，而这个方程组的求解其实就符合线性共轭梯度法，所以整个算法看起来就很像共轭梯度法，只不过这里迭代点是

p_k

，而对应的CG中的搜索方向就变成了

d_k

。

除去一般的过程和CG相同以外，算法中还有两个标红的地方。第一个地方说如果

d_j^T B d_j \le 0

就会终止。这是因为在方程组无解的情况下，这个矩阵是非正定的，也就可能会存在一个方向使得方向生成的二次型非正。那么这样的话，就说明我们CG的搜索方向已经无法逼近我的解了，所以就需要停止，返回当前的迭代点

，也就是对应的我们要的搜索方向。

第二个地方说满足某一个条件就返回，这个部分从表面上来看是为了保证迭代下降，提前终止使得算法复杂度降低，但是我们到后面会看到它也是为了保证收敛速度。

到这里，相信你已经明白了这个算法的运作原理，但是为什么它能保证我们得到的

是一个下降方向呢？注意到如果

d_j^T B d_j \le 0

就会终止，换句话说我们只会使用满足

d_i^T B d_i > 0, i = 0, 1, \cdots, j - 1

的

d_i

。那么如果我们设

s_i = Bd_i

，就会发现，我们会有

s_i^T d_i > 0

，且

s_j^T d_j \le 0

。所以可以看到，虽然我们在之前一直在考虑

\min_{x} \frac12 x^T Bx + g^T x

（这也是解

Bp = -g

的原因），但是实质上，我们只用到了前

步的信息，只考虑了前

步CG。为什么说这很重要呢？是因为我们在第3节

数值优化（3）——线搜索中的步长选取方法，线性共轭梯度法

有说明过这么一个性质：

Proposition 3:

x_k

是函数

f(x) = \frac12 x^T B x + g^T x

在

x_0 + \operatorname{span}(p_0, \cdots, p_{k-1})

这个空间上的最小值。

在这里因为我们的

p_0 = 0

。对应

x_0 = 0

，所以事实上，如果我们可以构造一个对称正定矩阵

\tilde B

，使得它落在这个线性空间，那么这个情况下，虽然我们的CG算法提前终止，但是得到的解依然是某一个函数的最小值。那么自然的，也就不用怀疑说它是一个下降方向了。

这样的矩阵是存在的，也就是

\tilde B = \sum_{i = 0}^{j - 1}\frac{s_is_i^T}{s_i^T d_i}

那么可以得到，这个时候的得到的方向

，其实就是函数

\frac12 x^T \tilde B x + g^T x

的极小值点，并且在搜索方向的时候，有

s_i = \tilde B d_i

，所以其实对于CG算法来说，它并无法察觉矩阵发生了改变。

有一点美中不足的地方在于它是对称半正定的，但是没有关系，因为对于一个

n \times n

的矩阵，如果它不满秩，那么一定可以通过谱分解，找到它的非零的特征向量。这些非零的特征向量可以组合成一组像空间的基。那么我们要做的就是补上核空间的基即可，也就是说考虑设

\hat B = \tilde B + Q_\perp Q_\perp^T

这里

Q_\perp

就是

\operatorname{span}(d_0, \cdots, d_{j - 1})

的正交补空间，补上这些之后，根据谱分解，不难得到

\hat B

就是正定对称阵，并且其实我们得到的

可以写成更好的表达式

p = -\hat B ^{-1} g

根据我们第二节提到的下降方向的充分必要条件，就可以发现

是一个下降方向。

有了这个东西之后，下面的就是正常的线搜索了，也就是说

x_{k + 1} = x_k + \alpha_k p_k

这里的步长选取条件为B-N条件即可，初始步长取1。

这一部分的理解其实还是有点难度的，需要线性代数的基的相关理论理解的比较到位。

好的，下面我们来解释一下标红的第二个地方究竟是什么含义。

Theorem 2: 设

\theta > 0, \kappa \in (0, 1)

，设

\nabla ^2 f(x)

在极小值

x^*

附近存在且连续，

\nabla^2f(x^*)

正定，且

\{x_k\}

收敛于

x^*

，那么牛顿CG算法具有超线性收敛速度。并且若

\nabla^2f(x)

在

x^*

附近Lipschitz连续，那么收敛阶为

\min(1 + \theta, 2)

。

我们证明一下这个结论。我们的想法是利用梯度的性质，因为梯度的本质就是两个函数的差的极限，所以根据中值定理，如果梯度是有限的数，那么它必然会与迭代点的差有相关的关系。

注意到我们有

\|r_k\| = \|\nabla ^2 f(x_k)p_k + \nabla f(x_k)\| \le \|r_0\|\min(\|r_0\|^\theta, \kappa)

那么根据

\|r_0\| = \|\nabla f(x)\|

，我们容易得到

\|r_k\| \le \|\nabla f(x_k)\|^{1 + \theta}

这样的话，因为我们希望看到的是关于梯度的变化情况，所以我们根据Taylor展开，会有

\nabla f(x_{k + 1}) = \nabla f(x_k) + \int_0^1 \nabla^2 f(x_k + tp_k)p_k dt \\ = r_k + \int_0^1 [\nabla^2 f(x_k + tp_k) - \nabla^2 f(x_k)]p_kdt

根究数分的知识，因为有连续性在，所以后面的一项就可以写成是

p_k

的小量，所以可以得到

\nabla f(x_{k + 1}) = r_k + o(\|p_k\|)

这样的话，综合在一起，就会得到

\|\nabla f(x_{k + 1})\| \le \|\nabla f(x_k)\|^{1 + \theta} + o(\|p_k\|)

现在我们相当于得到了一个

p_k

的表达式，所以还需要进一步化简它。注意到

\|p_k\| = \|\nabla^2 f(x_k)^{-1}(r_k - \nabla f(x_k))\|

因为我们要求函数足够靠近极值点，所以要求

\|\nabla f(x_k)\| \le 1

是不过分的。同时因为靠近，所以

\|\nabla^2 f(x_k)\|^{-1}

也必然存在上界。也就是说

\|p_k\| \le 2 \|\nabla^2 f(x_k)^{-1}\| \|\nabla f(x_k)\| \le 2L\|\nabla f(x_k)\|

因此代入这个

\|p_k\|

的估计，我们会有

\|\nabla f(x_{k + 1})\| \le (\|\nabla f(x_k)\|^\theta + o(1))\|\nabla f(x_k)\|

到此，整个表达式都已经只剩梯度这个变量了，只需要利用中值定理得到存在

\alpha, \beta > 0

，使得

\alpha\|x_k - x^*\| \le \|\nabla f(x_k)\| \le \beta \|x_k - x^*\|

成立，就可以得到最终的收敛速度的结论。

信赖域框架下的牛顿CG算法

我们在之前介绍了线搜索方法下实现的牛顿CG算法，但是有一个问题在于，我们虽然可以得到

p^*

为

\frac12 x^T \tilde B x + g^T x

在某一个线性空间下的极小值，但是如果说

\tilde B

是一个病态的矩阵，那么这个时候我们会发现

p^* = -(\tilde B)^+g

（表示

\tilde B

的广义逆）是不稳定的，且因为矩阵病态，矩阵的逆的条件数就会很大，就会导致这个

p^*

取得很长。很长就会导致线搜索的步长选取出现困难。要解决这个问题，我们在前一节提到的信赖域框架自然就派上用场了。也就是说我们牛顿CG方法的那个部分，其实就是希望求解下面这个问题

\min_{\|p\| \le \Delta}m(p) = f(x) + \nabla f(x)^T p + \frac12 p^T \nabla^2 f(x) p

这个问题和之前的二次凸问题没有本质上的区别，并且因为有了半径的约束，所以我们不会使得

p^*

取得太过分。

具体的算法形式可以见下面这几张图。

可以看出来，有变化的地方只有两个。一是一开始，如果在CG算法中的搜索方向（也就是

d_i

）出了问题，那么我们就会返回柯西点，而不再是单纯的负梯度。在半径过大的时候，我们会缩小我们的

\tau

，使得最终的我们的子问题满足半径的约束。

既然要说这样的方法好，自然我们需要了解一下它的收敛性。我们在之前有说，即使是柯西点，我们也能够得到一个信赖域方法的全局收敛性。所以如果我们能够说明函数的每一个点的效力都比柯西点要大（也就是说新的点会使得函数在二次模型的意义下取得更优的值），那么也就可以得到全局收敛性（这是因为更好的点也会满足那个不等式，具体的见上一节的Proposition 3）。为此我们先来观察迭代点的性质。

Proposition 4: 设

\{p_i\}

为牛顿CG方法中得到的迭代点，那么

m(p_i)

严格递减，并且

m(p^*) \le m(p_i), i = 1, \cdots, j

，假设到第

步迭代终止。并且有

\|p_i\|

严格递增，并且

\|p^*\| \ge \|p_i\|

我们证明一下这个结论。注意到我们的Proposition 2，也即

p_i

为

m(p) = \frac12 p^T B p + g^T p

在

\operatorname{span}(d_0, \cdots, d_{i - 1})

处的极小值，那么因为我们之后迭代的时候，会导致我的这个线性空间增大，所以函数值必然会下降（如果不下降，不妨设这一步为第

步，那么如果设

p_k = \sum_{i= 0 }^{k - 1}\alpha_i d_i

，那么会得到

\alpha_k d_k

这一步的权重是无效的，在线性CG算法中这只有可能对应着

\alpha_k = 0

，也就对应着

r_k = 0

，这就意味着最小值已经找到，所以迭代已经终止了，不会再有所谓的扩大线性空间了）。这就证明了结论。

接下来看一下我们迭代点模长的性质，根据我们的迭代公式，可以得到

\|p_{i + 1}\|^2 = \|p_i\|^2 + \alpha_i^2 \|d_i\|^2 + 2 \alpha_i d_i^T p_i

要说明我们的模长严格递增，只需要说明

d_i^T p_i > 0, i >= 1

即可。

根据初始值

p_0 = 0

，我们可以得到

p_0^T d_0 = 0

。现在我们假设

p_i^T d_i > 0

（有的人可能会问为什么不是假设

p_i^T d_i \ge 0

，这个我们后面再解释），那么考察一下

p_{i+1}^T d_{i + 1}

，可以得到

p_{i+1}^T d_{i + 1} = p_{i+1}^T (-r_{i + 1} + \beta_{i + 1}d_i)

注意到因为我们有

p_j = \sum_{i = 0}^{j - 1}\alpha_i d_i

，潜在的含义就是我们的

p_j

就是之前所有步搜索方向的线性组合，而每一步搜索方向都是与残差相互正交的（共轭性），所以可以得到

p_{i}^Tr_i = 0

，所以容易得到

p_{i+1}^T d_{i + 1} = p_{i+1}^T\beta_{i + 1}d_i = \beta_{i + 1}(p_i + \alpha_i d_i)^T d_i

再利用一下归纳假设和

\beta_{i + 1} > 0

即可得到结论。事实上从这里也可以看出，即使我们这里的第一步是

p_0^T d_0 = 0

，因为每一步都会有一个正的加成

\alpha_i d_i^T d_i

，所以之后的每一步的归纳假设依然是成立的，也就不影响我们的推导。

事实上，通过这样的推论，我们会发现，如果第一步就返回了，就会返回柯西点，如果没有在第一步返回，就会返回一个比柯西点更好的点。那自然不难证明全局收敛性了。

最后我们简单说两句有关局部收敛性的结论。事实上对于信赖域方法，也有完全相同的局部收敛性结果。一个重要的观察是在迭代点靠近极值点的时候，对应的信赖域方法会使得最终的

\rho_k

取得很大，这就意味着半径不会收缩，也就意味着这个时候的表现已经和线搜索方法无异。当然这样的说法是不严格的，只是严格的证明略显复杂，这里我们就略去不说了。

到此，我们就算介绍好了牛顿法相关的细节和内容。事实上针对牛顿法而言，使用牛顿CG加上信赖域的框架即可。而且事实上因为牛顿CG单独提供了求解子问题的框架，所以也不会用上开头提到的狗腿法啥的内容。

小结

本节主要介绍的是信赖域子问题的求解，当然重点放在了牛顿法和它的一系列推广上，我们通过介绍牛顿CG方法，连带着介绍了信赖域方法使用的一个具体例子。事实上，子问题的思路在拟牛顿法中得到了极为广泛的应用，这个我们下一节再说。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2020-05-25，如有侵权请联系 cloudcommunity@tencent.com 删除

编程算法

本文分享自学弱猹的精品小屋微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

编程算法

登录后参与评论

0 条评论

热度