数据科学和机器学习中的优化理论与算法（上）_线性规划,线性搜索和信赖域方法,无约束最优化,线性与非线性最小二乘问题,二次规划-CSDN博客

本文链接：https://blog.csdn.net/lusongno1/article/details/111475912

数据科学和机器学习中的优化理论与算法（上）

数据科学和机器学习当前越来越热，其中涉及的优化知识颇多。很多人在做机器学习或者数据科学时，对其中和优化相关的数学基础，包括随机梯度下降、ADMM、KKT 条件，拉格朗日乘数法、对偶问题等，不太了解，一知半解地用，用着用着就出错了。

本文希望从基础知识的角度，尽可能全地对最简单的优化理论和算法做一个小结。内容涵盖以下几个方面：优化简介、无约束优化、线搜索方法、信赖域方法、共轭梯度方法、拟牛顿方法、最小二乘问题、非线性方程、约束优化理论、非线性约束优化算法、二次规划、罚方法…

通过本文档的学习，相信你会掌握数据科学和机器学习中用到的优化基础知识，以后再遇到优化问题，就不会再困惑了。

建议收藏，方便时时查阅。

简介

数学背景

所谓的优化，在数学上说呢，就是最小化或者最大化一个函数，使得这个函数的变量满足一些约束。我们一般用 $x$ 向量来表示变量，称为未知数或者叫参数。 $f$ 叫做目标函数，它是我们想要最大化或者最小化的一个标量函数。 $c_i$ 我们称之为约束函数，用来定义和 $x$ 必须要满足的确定的等式或者不等式约束。
优化问题标准的数学方程一般写成下述形式：

$\begin{array}{ll} {\min _{x \in \mathbb{R}^{n}}} & {f(x)} \\ {\text { subject to }} & {c_{i}(x)=0, \quad i \in \mathcal{E}=\left\{1, \ldots, m_{e}\right\}}\\ & {c_{i}(x) \geq 0, \quad i \in \mathcal{I}=\left\{m_{e}+1, \ldots, m\right\}} \end{array} \tag{1}$

这里的 $\mathcal{E}$ 表示等式约束的指标集， $\mathcal{I}$ 表示不等式约束的指标集。规划问题又可分为线性规划和非线性规划，一般来说，只要目标函数或者约束函数当中有一个是非线性的，我们就称之为非线性规划。

凸性

一个集合 $\in \Re^{n}$ 被称为凸集，如果对于任意的两个点， $\in S,y\in S$ ，对任意的 $\alpha \in[0,1]$ ，我们有 $\alpha x+(1-\alpha) y \in S$ 。
我们说 $f$ 是个凸函数，当它的定义域是个凸集，并且下式满足。当下式对于 $x\neq y$ ，且 $\alpha$ 在开区间 $(0, 1)$ 上严格成立（取不到等号）时，我们称 $f$ 是严格凸的。若 $- f$ 是凸函数，那么 $f$ 就是一个凹函数。

$f(\alpha x+(1-\alpha) y) \leq \alpha f(x)+(1-\alpha) f(y) \quad \forall \alpha \in[0,1]$

所谓的凸优化，一般说的是在约束优化问题中，目标函数是凸的，等式约束是线性的，且不等式约束函数是凹函数。为什么不等式约束是凹函数呢？因为我们这里定义的不等式是 $c_{i}(x) \geq 0$ ，若改成 $\leq$ ，则是凸函数。

优化算法

优化算法一般都是一些迭代型的算法。所谓的迭代，就是先给定一个初值，然后通过一定的算法按照 $x_{k} \rightarrow x_{k+1}, \quad k=1,2, \ldots$ 的方式进行迭代更新。优化算法往往有一些评价指标。

鲁棒性（Robustness）：对于不同问题和初值都要是表现良好。
有效性（Efficiency）：不需要额外的时间和内存。
准确性（Accuracy）：能够准确地求解，不能对数据误差或者算法在数值计算时产生的舍入误差过度敏感。

除了以上几点，个人认为应该还添加上以下几点：收敛性、最优性、可扩展性、可靠性、可用性。

背景材料

向量和矩阵

向量的表示用 $\in \mathbb{R}^{n}: x=\left(x_{1}, \ldots, x_{n}\right)^{T}$ 。
向量內积：给定 $\in \mathbb{R}^{n}, x^{T} y=\sum_{i=1}^{n} x_{i} y_{i}$ 。
矩阵的表示用 $\in \mathbb{R}^{m \times n}$ 。
半正定矩阵： $x^{T} A x \geq 0$ 对于任意的 $\in \mathbb{R}^{n}$ 。
正交矩阵： $Q^{T} Q=Q Q^{T}=I$
特征值和特征向量： $x=\lambda x$ 。

向量范数

$\in \mathbb{R}^{n} \text { . }}$

$\begin{array}{l} {\qquad \begin{aligned} l_{1} \text { -norm: } &\|x\|_{1}=\sum_{i=1}^{n}\left|x_{i}\right| \\ l_{2} \text { -norm: } &\|x\|_{2}=\left(\sum_{i=1}^{n} x_{i}^{2}\right)^{1 / 2}=\left(x^{T} x\right)^{1 / 2} \\ l_{\infty} \text { -norm: } &\|x\|_{\infty}=\max _{i=1, \ldots, n}\left|x_{i}\right| \end{aligned}} \end{array}$

$\|x\|_{\infty} \leq\|x\|_{2} \leq \sqrt{n}\|x\|_{\infty}$ 且 $\|x\|_{\infty} \leq n\|x\|_{\infty}$ 。
柯西-斯瓦茨不等式： $\left|x^{T} y\right| \leq\|x\|_{2}\|y\|_{2}$ 。

对偶范数

$\|\cdot\|$ 的对偶范数：

$\|x\|_{D}=\max _{\|y\|=1} x^{T} y=\max _{y \neq 0} \frac{x^{T} y}{\|y\|}$

根据对偶范数的定义，容易知道 $\left|x^{T} y\right| \leq\|y\|\|x\|_{D}$ 。

矩阵范数

给定 $\in \mathbb{R}^{m \times n},$ 定义 $\|A\|=\sup _{x \neq 0} \frac{\|A x\|}{\|x\|}$ ，有

$\begin{array}{l} {\|A\|_{1}=\max _{j=1, \ldots, n} \sum_{i=1}^{m}\left|A_{i j}\right|} \\ {\|A\|_{2}=\operatorname{largest} \text { eigenvalue of }\left(A^{T} A\right)^{1 / 2}} \\ {\|A\|_{\infty}=\max _{i=1, \ldots, m} \sum_{j=1}^{n}\left|A_{i j}\right|} \end{array}$

Frobenius范数：

$\|A\|_{F}=\left(\sum_{i=1}^{m} \sum_{j=1}^{n} A_{i j}^{2}\right)^{1 / 2}$

条件数： $\kappa(A)=\|A\|\left\|A^{-1}\right\|$ 。

子空间

所谓子空间，给定 $\mathcal{S} \subset \mathbb{R}^{n}$ ，若 $\in \mathcal{S}$ ，那么 $\alpha x+\beta y \in \mathcal{S},$ 对于所有的 $\alpha, \beta \in \mathbb{R}$ 。
零空间： $\operatorname{Null}(A)=\{w | A w=0\}$ ，值空间： $\operatorname{Range}(A)=\{w | w=A v \text { for some vector } v\}$ ，零空间和转置的值空间的直和为全空间：Null $\oplus$ Range $\left(A^{T}\right)=\mathbb{R}^{n}$ 。

导数

Frechet 可微性:
$\mathbb{R}^{n} \rightarrow \mathbb{R},$ 在 $x$ 处是可微的，若 $\exists g \in \mathbb{R}^{n}$ 使得：

$\lim _{y \rightarrow 0} \frac{f(x+y)-f(x)-g^{T} y}{\|y\|}=0$

$f$ 的梯度：

$g(x)=\nabla f(x)=\left(\frac{\partial f}{\partial x_{1}}, \ldots, \frac{\partial f}{\partial x_{n}}\right)^{T} \in \mathbb{R}^{n}$

$f$ 的 Hessian 矩阵：

$H(x)=\nabla^{2} f(x)=\left[\frac{\partial^{2} f}{\partial x_{i} \partial x_{j}}\right] \in \mathbb{R}^{n \times n}$

链式法则 1：

$\frac{d \phi(\alpha(\beta))}{d \beta}=\phi^{\prime}(\alpha) \alpha^{\prime}(\beta)$

链式法则 2： $\in \mathbb{R}^{n}$ 且 $x = x (t)$ ，定义 $h (t) = f (x (t))$ ，那么

$\nabla h(t)=\sum_{i=1}^{n} \frac{\partial f}{\partial x_{i}} \nabla x_{i}(t)=\nabla x(t)^T \nabla f(x(t))$
方向导数： $p)=\lim _{\epsilon \rightarrow 0} \frac{f(x+\epsilon p)-f(x)}{\epsilon}=\nabla f(x)^{T} p$ 。

收敛率

令 $\left\{x_{k}\right\}$ 是 $\mathbb{R}^{n}$ 的一个序列，收敛到 $x^{*}$ 。

Q-线性收敛说的是存在一个常数 $r\in(0,1)$ ，使得当 $k$ 充分大时，有：

$\frac{\left\|x_{k+1}-x_{*}\right\|}{\left\|x_{k}-x^{*}\right\|} \leq r$
Q-超线性收敛指的是：

$\lim _{k \rightarrow \infty} \frac{\left\|x_{k+1}-x^{*}\right\|}{\left\|x_{k}-x^{*}\right\|}=0$

Q-二次收敛指的是，存在一个常数 $M$ ，使得当 $k$ 充分大时，

$\frac{\left\|x_{k+1}-x_{*}\right\|}{\left\|x_{k}-x^{*}\right\|^{2}} \leq M$

无约束优化

$\Re^{n} \rightarrow \Re$ 是一个光滑函数，无约束优化指的就是没有约束的优化。

$\min _{x \in \mathbb{R}^{n}} f(x)$

不同说法的解

全局极小值指的是

$f\left(x^{*}\right) \leq f(x) \quad \text { for all } x \in \mathbb{R}^{n}$

局部极小值指的是存在 $x^{*}$ 的一个领域 $\mathcal{N}$ ，

$f\left(x^{*}\right) \leq f(x) \quad \text { for all } x \in \mathcal{N}$

严格的局部极小值，

$f\left(x^{*}\right)<f(x) \quad \text { for all } x \in \mathcal{N} \text { with } x \neq x^{*}$

孤立的局部极小值，当 $x^{*}$ 是领域 $\mathcal{N}$ 中唯一的局部极小值。

识别局部极小值

泰勒定理

首先介绍一个泰勒定理。假定 $\Re^{n} \rightarrow \Re$ 是一个连续可微函数，并且 $\in \Re^{n}$ ，那么，我们有

$f(x+p)=f(x)+\nabla f(x+t p)^{T} p$

对某个 $\in(0,1)$ 。进一步，如果 $f$ 是二次可微函数，我们有

$\nabla f(x+p)=\nabla f(x)+\int_{0}^{1} \nabla^{2} f(x+t p) p d t$

并且

$f(x+p)=f(x)+\nabla f(x)^{T} p+\frac{1}{2} p^{T} \nabla^{2} f(x+t p) p$

对某个 $\in(0,1)$ 。

一、二阶必要性条件

一阶必要性条件：如果 $x^{*}$ 是一个局部极小值，并且 $f$ 在 $x^{*}$ 的某个开领域内时连续可微的，那么 $\nabla f\left(x^{*}\right)=0$ 。
二阶必要性条件：如果 $x^{*}$ 是个局部极小值点，并且 $f$ 和 $\nabla^{2} f(x)$ 存在且在 $x^{*}$ 的某个开领域中是连续的，那么 $\nabla f\left(x^{*}\right)=0$ 并且 $\nabla^{2} f\left(x^{*}\right)$ 是个半正定的矩阵。
稳定点表示满足 $\nabla f\left(x^{*}\right)=0$ 的 $x^{*}$ 点。

二阶充分性条件

假定 $\nabla^{2} f(x)$ 在 $x^{*}$ 的开领域中是连续的， $\nabla f\left(x^{*}\right)=0$ 且 $\nabla^{2} f\left(x^{*}\right)$ 正定的，那么 $x^{*}$ 是 $f$ 的一个严格极小值点。

从这个条件中可以看出，这里必要性条件，有个更高的要求，要求二阶导是正定的，得到的结论也更强，是严格的局部极小值。这个条件也只是充分的，而不是必要的条件，也就是说存在严格局部极小值点，不一定满足正定这个条件，但是半正定一定要满足，这是必要性条件告诉我们的。

如果 $f$ 是个凸函数，那么局部极小值就是全局极小值。进一步，如果 $f$ 是可微的，那么任何的稳定点都是全局极小值点。对于迭代方法，我们一般找到的是梯度消失的点，也就是稳定点。

无约束优化问题的常用算法

求解这种无约束优化问题，常用的有线搜索方法和信赖域方法，这两种方法后面会详细地提到。简单地提一句，线搜索方法其实就是每次迭代都沿着一个方向（类似于一根线）走一定的步长，等走到那个点了，再用类似的方法走到下一个点。

信赖域方法，我国的第一代留学生------孙悟空，就曾经用过此法。唐僧取经那一会儿，总归要有人出去化缘，八戒这货喜欢偷懒，经常赖着不愿意出去要饭，师傅总是要吃饭哪，所以，只能猴哥出去化缘。猴子出去了，但是又怕师傅会被妖怪吃掉，只能画个圈圈，把师徒三人圈在圈圈里面，并叮嘱他们不要出来，免得要妖怪吃掉。这个圈圈啊，其实就是我们所说的"信赖域"。孙悟空的法力值就那么大，圈圈不能画太大，画太大就没什么威力了，容易被妖怪打破。师徒三人呢，活动也不能超过这个范围，如果跑出这个范围呢，就会被妖怪吃掉。如果被妖怪吃掉了，就取不成经，取不成经，就不"收敛"了

线搜索方法简述

所谓的线搜索方法，就是说在每一次迭代我们得先找一个搜索方向 $p_k$ ，然后沿着搜索方向走一步，

$x_{k+1}=x_{k}+\alpha_{k} p_{k}$

这里的 $\alpha_k$ 我们一般称之为步长，在机器学习当中，我们一般也叫它学习率。算法描述如图。

线搜索方法

梯度下降方向

如果选择下降方向呢，一个比较流行的选择就是最速下降方向，也就是负梯度方向。 $-\nabla f_{k}$ 是当前点 $f$ 下降最快的方向，需要注意这只是当前点的下降最快的方向，所以它不一定是收敛最快的方向。为什么负梯度方向是下降最快的方向呢？从泰勒展开就可以看出来了。

$f\left(x_{k}+\alpha p\right)=f\left(x_{k}\right)+\alpha p^{T} \nabla f_{k}+\frac{\alpha^{2}}{2} p^{T} \nabla^{2} f\left(x_{k}+t p\right) p$

忽略二次项的高阶小量，欲令前后两步的差值尽量大， $p$ 自然取为负梯度方向。一般满足 $p^{T} \nabla f_{k}<0$ 的方向，我们都叫下降方向。

在机器学习中，我们经常提到随机梯度下降。我们知道，在做有监督的问题中，有一个损失函数，用来衡量模型的好坏。在回归问题中，如果我们已经知道方程的形式，我们可以做出问题的损失函数。那么我们要做的就是寻找方程的最佳参数，即寻找参数，使得损失函数函数达到最小。乍一听，这就是个优化问题。当然可以对参数的各个分量求导，最后解方程组。问题是，你能确定最后的方程组好解？一般的数值解法有梯度下降法、牛顿法、拟牛顿方法以及信赖域方法等等。

梯度下降方法简单易行（大部分情况下还是收敛挺快），可是当数据量比较大时，梯度下降的计算就特别大。为解决大数据量的优化求解，我们一般可以采用随机梯度下降。
所谓的随机梯度下降，就是相比梯度下将，我们不选择全部的数据点，而是每次随机选择一个或者若干个数据点来构建损失函数，接着再用梯度下降方法求解。看起来很不靠谱的样子，事实上多做几次这个操作，它是可以收敛的最优解的某个邻域的。

牛顿和拟牛顿方向

做二阶泰勒近似，

$f\left(x_{k}+p\right) \approx f_{k}+p^{T} \nabla f_{k}+\frac{1}{2} p^{T} \nabla^{2} f_{k} p \equiv m_{k}(p)$

我们希望找到一个 $p$ ，使得 $m_k$ 最小，求解这个二次函数的最小值，可得，

$p_{k}^{N}=-\left(\nabla^{2} f_{k}\right)^{-1} \nabla f_{k}$

这便是传说中的牛顿方向。使用牛顿方向，就称为牛顿方法。

有时候，二阶 Hessian 阵不是那么好算，我们会找一个 $B_k$ 去逼近 $\nabla^{2} f_{k}$ ，就得到了拟牛顿方向，

$p_{k}=-B_{k}^{-1} \nabla f_{k}$

那么这个近似方向 $B_k$ 怎么选呢？不要着急，后面还有专门的章节来讨论拟牛顿方法。比如说，BFGS方法， $B_k$ 就是这样迭代的，

$B_{k+1}=B_{k}-\frac{B_{k} s_{k} s_{k}^{T} B_{k}}{s_{k}^{T} B_{k} s_{k}}+\frac{y_{k} y_{k}^{T}}{y_{k}^{T} s_{k}}$

这里， $s_{k}=x_{k+1}-x_{k}, \quad y_{k}=\nabla f_{k+1}-\nabla f_{k}$ ，不知所以然，没关系，后面还会慢慢提这个事情。简单地说， $B_k$ 在满足割线方程 $B_{k+1} s_{k}=y_{k}$ 的前提下，再满足一些别的约束，就能导出不同的关于 $B_k$ 的迭代公式。

非线性共轭梯度方向

当然，除了上述提到的最速下降方向、牛顿和拟牛顿方向，还有一些其他的线搜索方向。非线性共轭梯度方法诱导的方向就是其中之一。

$p_{k}=-\nabla f\left(x_{k}\right)+\beta_{k} p_{k-1}$

这是负梯度方向和前一个搜索方向做的一个组合，这里的 $\beta_k$ 是保证 $p_k$ 和 $p_{k-1}$ 共轭的一个常数。什么叫共轭，不急，后面还会详细说。

步长选择

不管是负梯度方向，还是牛顿方向，我们只是选择了方向。步长如何选取呢？一个理想的选择，就是沿着线搜索方向，找到函数最小的点对应的步长作为搜索步长。即最小化如下函数：

$\phi(\alpha)=f\left(x_{k}+\alpha p_{k}\right), \quad \alpha>0$

因为终归是在这个方向上找到了精确的最小值点，所以一般就称为精确线搜索步长。对应于精确线搜索的，就是非精确线搜索。也就是找一堆步长值去试，直到满足一定的终止条件。
使用非精确线搜索步长，步长一定要满足一定的条件。一个简单的例子就是 $f=x^2$ ，如果步长没有选取好，即使是梯度下降方法，迭代也有可能在最小值附近一左一右跑，呈"之"字型荡来荡去，收敛极慢，或者索性就不收敛了。就不画图了，应该能想象的到。为了避免这种情况，我们需要对步长施加一些充分下降条件。充分下降条件就是保证选定的步长，使得 $f$ 值是下降的，

$f\left(x_{k}+\alpha p_{k}\right) \leq f\left(x_{k}\right)+c_{1} \alpha \nabla f_{k}^{T} p_{k}:=l(\alpha)$

这里的 $c_1 \in (0,1)$ ，一般可以取 $10^{-4}$ 。这个条件告诉我们， $f$ 的下降要和步长和方向导数成比例。满足充分下降条件还是不够的，因为只要足够小的步长 $\alpha$ 都满足这个条件，我们还需要一些条件筛掉不合理的小步长。于是便有了曲率条件，

$\nabla f\left(x_{k}+\alpha_{k} p_{k}\right)^{T} p_{k} \geq c_{2} \nabla f_{k}^{T} p_{k}$

这里 $c_{2} \in\left(c_{1}, 1\right)$ ，在牛顿和拟牛顿方法中，一般取为0.9，非线性CG方法中，一般取为0.1。 $c_1$ 来自于充分下降条件。如果定义 $\phi(\alpha)=f\left(x_{k}+\alpha p_{k}\right)$ ，那么曲率条件也可以写成 $\phi^{\prime}\left(\alpha_{k}\right) \geq c_{2} \phi^{\prime}(0)$ 。也就是说走一步之后，原下降方向就要离负梯度方向更远一些，下降得就没那么厉害一些。

所谓的 Wolfe 条件就是把充分下降条件和曲率条件放在一块，

$\begin{aligned} f\left(x_{k}+\alpha_{k} p_{k}\right) & \leq f\left(x_{k}\right)+c_{1} \alpha_{k} \nabla f_{k}^{T} p_{k} \\ \nabla f\left(x_{k}+\alpha_{k} p_{k}\right)^{T} p_{k} & \geq c_{2} \nabla f_{k}^{T} p_{k} \end{aligned}$

强 Wolfe 条件，就是在曲率条件上加个绝对值，它要求 $\alpha_k$ 满足，

$\begin{aligned} f\left(x_{k}+\alpha_{k} p_{k}\right) & \leq f\left(x_{k}\right)+c_{1} \alpha_{k} \nabla f_{k}^{T} p_{k} \\ \left|\nabla f\left(x_{k}+\alpha_{k} p_{k}\right)^{T} p_{k}\right| & \leq c_{2}\left|\nabla f_{k}^{T} p_{k}\right| \end{aligned}$

同样， $0<c_{1}<c_{2}<1$ 。说了 Wolfe 条件和强 Wolfe 条件，那么是否一定能找到这种苛刻条件的步长呢？答案是肯定的。定理表明，对于光滑有界的函数而言，总是能找到满足 Wolfe 条件和强 Wolfe 条件的步长的。

还有一个条件叫做 Goldstein 条件，它也能保证充分下降，且步长不太小。

$f\left(x_{k}\right)+(1-c) \alpha_{k} \nabla f_{k}^{T} p_{k} \leq f\left(x_{k}+\alpha_{k} p_{k}\right) \leq f\left(x_{k}\right)+c \alpha_{k} \nabla f_{k}^{T} p_{k}$

这里 $0<c<\frac{1}{2}$ ，右半部分是充分下降条件，左半部分保证了步长不太小。

我们可以利用充分下降条件，构造回溯算法。所谓的回溯法，就是开始的时候，选个比较大的步长，慢慢减小之，直到它满足了充分下降条件，立即停止，这样，也使得步长不太小了。

线搜索方法的收敛性

介绍一下 Zoutendijk 条件： $\sum_{k \geq 0} \cos ^{2} \theta_{k}\left\|\nabla f_{k}\right\|^{2}<\infty$ 。这里的 $\theta_k$ 是最速下降方向和线搜索方向的夹角，即

$\cos \theta_{k}=\frac{-\nabla f_{k}^{T} p_{k}}{\left\|\nabla f_{k}\right\|\left\|p_{k}\right\|}$

有定理表明，在一些本质条件下，若 Wolfe 条件被满足，那么 Zoutendijk 条件总是成立的。Zoutendijk 条件表明 $\cos ^{2} \theta_{k}\left\|\nabla f_{k}\right\|^{2} \rightarrow 0$ 。容易知道，不管是梯度下降方法还是牛顿方法（ $B_k$ 正定）， $\cos \theta_{k}$ 都是有下界的，也就是梯度方向和搜索方向总是不垂直，那么必然有，

$\lim _{k \rightarrow \infty}\left\|\nabla f_{k}\right\|=0$

即算法至少收敛到一个稳定点，我们说收敛到梯度为0的点，叫做"全局收敛"的。下面给出改进的线搜索牛顿法的一个算法描述，如图。

改进线搜索方法

这里保证 $B_k$ 正定，是为了保证收敛。事实上，定理表明，只要保证 $B_k$ 的条件数足够小就可以了。

最速下降方法的收敛率依赖于 $f$ 的 Hessian 的条件数，一般来说不超过线性收敛的速度，一些特殊情况下，收敛速度慢得可怜。对于牛顿方法，满足一定的不太严苛的条件，一般是全局收敛的。若满足二阶充分性条件，当初值选得足够靠近真值时，一般就是二次收敛的。同样满足二阶充分性条件，拟牛顿方法就会差一些，它是超线性收敛的。

信赖域方法

信赖域方法（Trust-region methods）又称为TR方法，它是一种最优化方法，能够保证最优化方法总体收敛。现今，信赖域算法广泛应用于应用数学、物理、化学、工程学、计算机科学、生物学与医学等学科。相信在不远将来，信赖域方法会在更广泛多样的领域有着更深远的的发展。这里简单以介绍信赖域狗腿（dogleg）方法，让大家体会一下什么叫信赖域方法，希望能为其他专业诸如机器学习的同学在使用优化工具时提供一个参考。如前所述，信赖域方法其实并不新鲜，我国的第一代留学生------孙悟空，早就用过此法，也就是我们常说的"画个圈圈诅咒你"。

原理描述

考虑模型问题：

$\min _{p \in \mathbb{R}^{n}} m(p)=f+g^{T} p+\frac{1}{2} p^{T} B p, \quad \text { s.t. }\|p\| \leq \Delta$

它具有全局最优解：

$p^{B}=-B^{-1} g$

该问题沿着负梯度方向，具有全局最优解：

$p^{U}=-\frac{g^{T} g}{g^{T} B g} g$

由于信赖域的限制，全局最优解可能落在信赖域之外。我们依据两个全局最优点是否落于信赖域中，来决定下降方向。

当B点在信赖域内部时，取方向为 $p^U$ 方向， $x_{k+1}$ 点为B点：
当B点和U点都在信赖域外时，取下一个迭代点为 $p^U$ 和信赖域边界的交点：

在这里插入图片描述

当B点在外，U点在内时，去BU连线和信赖域边界的交点：

在这里插入图片描述

我们可以用一个参数 $\tau$ 将这三种情况统一为一个表达式，即下降取值为：

$\tilde{p}(\tau)=\left\{\begin{array}{ll} {\tau p^{U},} & {0 \leq \tau \leq 1} \\ {p^{U}+(\tau-1)\left(p^{B}-p^{U}\right),} & {1 \leq \tau \leq 2} \end{array}\right.$

这里的 $\tau$ 是如何取的呢？通过简单的数学推导，我们知道它可以这样来取，刚好满足上面的三种情况。简单地说，第一种情况， $\tau$ 取2，第二种情况， $\tau$ 取为信赖域半径和 $p^U$ 长的比值，第三种情况，要求那个交点，我们就要求一元二次方程的一个根，这可以用求根公式算得。

算法步骤

总的算法框架如图。

在这里插入图片描述

$\tau$ 和 $s_k = x_{k+1} - x_k$ 的计算就是使用上面提到的狗腿方法。下降比 $\rho_k$ 的计算表达式为：

$\rho_{k}=\frac{f\left(x_{k}\right)-f\left(x_{k}+s_{k}\right)}{m_{k}(0)-m_{k}\left(s_{k}\right)}$

所以，到这事情就很明了了。一个简单的代码如下：

function [x,k] = trust_region1(varargin)
%帮助信息：
%clc
%clear
%输入参数为options = struct('f',f,'x0',[0;0],'Delta',1,'Delta0',0.1,'eta',1/4);
if nargin < 1, help(mfilename),
    f = @(x1,x2)100*(x2-x1^2)^2 + (1-x1)^2;
    options = struct('f',f,'x0',[0;0],'Delta_hat',10,'Delta0',0.1,'eta',0);
    fprintf('使用默认参数…… \n \n');
    addpath(genpath(pwd));
end

f = options.f;
syms x1 x2;
f_s = f(x1,x2);
x0 = options.x0;
Delta_hat = options.Delta_hat;
Delta0 = options.Delta0;
eta = options.eta;
df = diff_handle(f_s);
B = hessian(f_s);
Delta = Delta0;
x = x0;
step = 100;
for k = 0:step-1
    x
    fk = f(x(1),x(2));
    dfk = df(x(1),x(2));
    Bk = subs(B,{x1,x2},{x(1),x(2)});
    Bk = double(Bk);
    sk = cal_sk(dfk,Bk,Delta);
    m = @(p) fk + dfk'*p + 1/2*p'*Bk*p;
    rho_k = cal_rho_k(f,m,x,sk);
    if rho_k < 1/4
        Delta = 1/4*Delta;
    elseif rho_k > 3/4 && abs(norm(sk,2) - Delta)<1.0e-10;
        Delta = min(2*Delta,Delta_hat);
    else
        %不做任何操作；
    end

    if rho_k > eta
        x = x + sk;
    else
        %不做任何操作；
    end
end

end








function df = diff_handle(f_s)
syms x1 x2;
df = [diff(f_s,x1); diff(f_s,x2)];
df = matlabFunction(df);
end



function tau = cal_tau(pB,pU,Delta)
npB = sqrt(pB'*pB);
npU = sqrt(pU'*pU);
if npB <= Delta
    tau = 2;
elseif npU >= Delta
    tau = Delta/npU;
else
    pB_U = pB-pU;
    tau = (-pU'*pB_U+sqrt((pU'*pB_U)^2-pB_U'*pB_U*(pU'*pU-Delta^2)))/(pB_U'*pB_U);
    tau = tau + 1;
end
end
% temp = conj(dfk')*Bk*dfk;
% if  temp<= 0
%     tau = 1;
% else
%     tau = min(norm(dfk,2)/temp,1);




function sk = cal_sk(dfk,Bk,Delta)
pU = -(conj(dfk')*dfk)/(conj(dfk')*Bk*dfk).*dfk;
pB = -Bk^(-1)*dfk;
tau = cal_tau(pB,pU,Delta);
if tau >=0 && tau <=1
    sk = tau*pU;
elseif tau >= 1 && tau <=2
    sk = pU + (tau-1)*(pB-pU);
else
    error('tau的值不能为%f',tau);
end

end


function rho_k = cal_rho_k(f,m,x,sk)
rho_k = (f(x(1),x(2)) - f(x(1)+sk(1),x(2)+sk(2)))/((m([0;0])-m(sk)));
end

一个例子

使用狗腿方法来求解信赖域子问题，应用这个方法来最小化下述函数。这里的 $B_k = \nabla^2 f(x_k)$ 。使用信赖域方法的更新规则更新。这里给出前两步迭代结果。

$f(x)=100\left(x_{2}-x_{1}^{2}\right)^{2}+\left(1-x_{1}\right)^{2}$

选定的参数为：信赖域上界取为10，信赖域初始半径取为0.1， $\eta$ 取为0。在此参数下，迭代的前两个 $x$ 值为： $x_0=[0.1000,0],x_1=[0.2933,0.0513]$ 。

在一些条件的保证下，信赖域方法能够收敛到问题的一个稳定点。如果在局部能够满足二阶充分性条件，那么信赖域方法在这个局部，能有一个超线性收敛的速度。在信赖域方法，还有一个相场被提到的概念叫"柯西点"。所谓的柯西点，我们用 $s_{k}^{c}$ 来表示，它是在信赖域约束下的一个负梯度方向最小值点，

$\begin{array}{rl} {s_{k}^{c}=\arg \min } & {m_{k}(p)} \\ {\text { s.t. }} & {p=-\tau \nabla f_{k}} \\ {} & {\|p\| \leq \Delta_{k}, \tau \geq 0} \end{array}$