最优化学习笔记：无约束优化算法（6）

6.6 信赖域算法

本节介绍信赖域算法．它和线搜索算法都是借助泰勒展开来对目标函数进行局部近似，但它们看待近似函数的方式不同．在线搜索算法中，我们先利用近似模型求出下降方向，然后给定步长；而在信赖域类算法中，我们直接在一个有界区域内求解这个近似模型，而后迭代到下一个点．因此信赖域算法实际上是同时选择了方向和步长．

6.6.1 信赖域算法框架

我们对信赖域算法给一个直观的数学表达. 根据带拉格朗日余项的泰勒展开，
$$f(x^k+d)=f(x^k)+\nabla f(x^k{)}^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}{\nabla }^{2}f(x^k+td)d$$

其中 $t\in (0,1)$ 为和 $d$ 有关的正数. 和牛顿法相同，可以利用 $f(x)$ 的一个二阶近似来刻画 $f(x)$ 在点 $x=x^k$ 处的性质：
$$m_k(d)=f(x^k)+\nabla f(x^k{)}^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}{B}^{k}d(\mathrm{6.6.1})$$

其中 $B^k$ 是对称矩阵. 并且要求 $B^k$ 是海瑟矩阵的近似矩阵，如果 $B^k$ 恰好是函数 $f(x)$ 在点 $x=x^k$ 处的海瑟矩阵，那么当 $f(x)$ 充分光滑时，$m_k(d)$ 的逼近误差是 $O(\Vert d{\Vert }^3)$.

由于泰勒展开的局部性质，它仅仅对模长较小的 $d$ 有意义. 当 $d$ 过长时，模型 (6.6.1) 便不再能刻画 $f(x)$ 的特征，为此需要对模型 (6.6.1) 添加约束. 我们仅在如下球内考虑 $f(x)$ 的近似：
Ω k = { x k + d ∣ ∥ d ∥ ⩽ Δ k } \Omega_{k}=\{x^{k}+d\mid\lVert d\rVert\leqslant\Delta_{k}\} Ωk​={xk+d∣∥d∥⩽Δk​}

其中 ${\Delta }_k>0$ 是一个和迭代有关的参数. 称 ${\Omega }_k$ 为信赖域，${\Delta }_k$ 为信赖域半径. 从命名方式也可看出，信赖域就是我们相信 $m_k(d)$ 能很好地近似 $f(x)$ 的区域，而 ${\Delta }_k$ 表示了这个区域的大小．因此信赖域算法每一步都需要求解如下子问题
$$\underset{d\in 1^n}{\min }{m}_k(d),\text{s.t.}\Vert d\Vert ⩽{\Delta }_k(\mathrm{6.6.2})$$

图 6.15 显示了子问题 (6.6.2) 的求解过程. 图中实线表示 $f(x)$ 的等高线，虚线表示 $m_k(d)$ 的等高线 (这里有关系 $d=x-x^k$), $d_N^k$ 表示求解无约束问题得到的下降方向 (若 $B^k$ 为海瑟矩阵则 $d_N^k$ 为牛顿方向), $d_{TR}^k$ 表示求解信赖域子问题得到的下降方向，可以看到二者明显是不同的. 信赖域算法正是限制了 $d$ 的大小，使得迭代更加保守，因此可以在牛顿方向很差时发挥作用.
在子问题 (6.6.2) 中仍需要确定信赖域半径 ${\Delta }_k$. 实际上，选取信赖域半径非常关键，它决定了算法的收敛性. 考虑到信赖域半径是“对模型 $m_k(d)$ 相信的程度”, 如果 $m_k(d)$ 对函数 $f(x)$ 近似较好，就应该扩大信赖域半径，在更大的区域内使用这个近似，反之就应该减小信赖域半径重新计算.
引入如下定义来衡量 $m_k(d)$ 近似程度的好坏：
$${\rho }_k=\frac{f(x^k)-f(x^k+d^k)}{m_k(0)-m_k(d^k)}(\mathrm{6.6.3})$$

其中 $d^k$ 为求解信赖域子问题得到的迭代方向. 根据 ${\rho }_k$ 的定义可知，它是函数值实际下降量与预估下降量 (即二阶近似模型下降量) 的比值. 如果 ${\rho }_k$ 接近 1, 说明用 $m_k(d)$ 来近似 $f(x)$ 是比较成功的，我们应该扩大 ${\Delta }_k$；如果 ${\rho }_k$ 为非常小的正数甚至为负数，就说明我们过分地相信了二阶模型 $m_k(d)$, 此时应该缩小 ${\Delta }_k$. 使用这个机制可以动态调节 ${\Delta }_k$, 让二阶模型 $m_k(d)$ 的定义域处于一个合适的范围.

算法 6.8 给出完整的信赖域方法. 该算法虽然有一些参数，但是它对这些参数的取值并不敏感. 实际中可取 ${\bar{\rho }}_1=0.25,{\bar{\rho }}_2=0.75$ 以及 ${\gamma }_1=0.25,{\gamma }_2=2$. 要注意，信赖域半径 ${\Delta }_k$ 不会无限增长，一是因为它有上界的控制，二是如果信赖域约束不起作用 (即二次模型最优值处于信赖域内), 我们也无需增加信赖域半径. 只有当 $m_k(d)$ 近似足够好并且信赖域约束起作用时，才需要增加 ${\Delta }_k$.

# 算法 6.8 信赖域算法
import numpy as np
def m_k(x, d):
    return f(x) + np.dot(grad_f(x), d) + 0.5 * np.dot(d, np.dot(hess_f(x), d))
def trust_region_algorithm(x_0, Delta_0, Delta_max, eta, rho1_bar = 0.25, rho2_bar = 0.75, gamma1 = 0.25, gamma2 = 2, max_iter = 1000, tol = 1e-6):
    x = x_0
    Delta = Delta_0
    k = 0
    x_history = [x]
    while k < max_iter:
        dk = calculate_iterative_direction(x, Delta)  # 计算迭代方向
        mk_0 = m_k(x, np.zeros_like(x)) 
        mk_dk = m_k(x, dk)  
        rho_k = (f(x) - f(x + dk)) / (mk_0 - mk_dk)  # 计算下降率
        # 更新信赖域半径
        if rho_k < rho1_bar:
            Delta = gamma1 * Delta
        elif rho_k > rho2_bar and np.linalg.norm(dk) == Delta:
            Delta = min(gamma2 * Delta, Delta_max)
        else:
            Delta = Delta
        # 更新自变量
        if rho_k > eta:
            x = x + dk
        x_history.append(x)
        k += 1
        if np.linalg.norm(grad_f(x)) <= tol:
            break    
    return x_history, k

def calculate_iterative_direction(x_k, Delta_k):
    # 计算牛顿方向
    hessian = hess_f(x_k)
    grad = grad_f(x_k)
    # 求解线性方程组 Hessian * d_k = -grad
    d_k = np.linalg.solve(hessian, -grad)
    # 对牛顿方向进行修正，确保在信赖域范围内
    if np.linalg.norm(d_k) > Delta_k:
        d_k = (Delta_k / np.linalg.norm(d_k)) * d_k
    return d_k

def f(x):
    return x[0]**2 + 10 * x[1]**2
def grad_f(x):
    return np.array([2*x[0], 20*x[1]])
def hess_f(x):
    return np.array([[2, 0], [0, 20]])

x_0 = np.array([-10.0, -1.0])  # 初始点
Delta_0 = 1.0  # 初始半径
Delta_max = 10.0  # 最大半径
eta = 0.1
# 调用信赖域算法
trust_history, k = trust_region_algorithm(x_0, Delta_0, Delta_max, eta)
print("Optimal solution:", trust_history[-1])
'''Optimal solution: [0.00000000e+00 1.38777878e-17]'''
print("迭代次数：", k)
'''迭代次数： 5'''
# 绘制等高线图
x = np.linspace(-12, 12, 100)
y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + 10*Y**2 
plt.figure(figsize=(10, 5))
plt.contour(X, Y, Z, levels=np.logspace(-2, 3, 20), cmap='jet')  # levels 参数控制等高线的密集程度
plt.plot([x[0] for x in trust_history], [x[1] for x in trust_history], markersize = 3, marker='o', markerfacecolor='none', linewidth=1, color='orange', label='trust_region_algorithm')
plt.plot([x[0] for x in M_Newton_history], [x[1] for x in M_Newton_history], markersize = 3, marker='o', markerfacecolor='none', linewidth=1, color='r', label='modified_newton')
plt.plot([x[0] for x in quasi_Newton_history], [x[1] for x in quasi_Newton_history], markersize = 3, marker='x', markerfacecolor='none', linewidth=1, color='b', label='Quasi-Newton')
plt.legend(loc='best', markerscale=1.5)  # 设置标签位置为最佳位置，调整图例的大小
plt.show() 

6.6.2 信赖域子问题求解

在多数实际应用中，信赖域子问题的解是无法显式写出的．为了求出迭代方向 $d^k$，需要设计算法快速或近似求解子问题 (6.6.2)．

1. 迭代法
信赖域子问题是一个仅仅涉及二次函数的约束优化问题。下面的定理给出了子问题的最优解需要满足的条件：

定理 6.12 $d^{\ast }$ 是信赖域子问题
$$\begin{array}{r}\min m(d)=f+g^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}Bd,s.t.\Vert d\Vert ⩽\Delta \end{array}(\mathrm{6.6.4})$$

的全局极小解当且仅当 $d^{\ast }$ 是可行的且存在 $\lambda ⩾0$ 使得
$$(B+\lambda I)d^{\ast }=-g(\mathrm{6.6.5}a)$$

$$\lambda (\Delta -\Vert d^{\ast }\Vert )=0(\mathrm{6.6.5}b)$$

$$(B+\lambda I)⪰0(\mathrm{6.6.5}c)$$

证明
‘必要性’（利用 KKT 条件来直接写出 $d^{\ast }$ 所满足的关系）
问题(6.6.4)的拉格朗日函数为
$$L(d,\lambda )=f+g^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}Bd-\frac{\lambda }{2}({\Delta }^2-\Vert d{\Vert }^2),$$

其中乘子 $\lambda ⩾0$. 根据 KKT 条件，$d^{\ast }$ 是可行解，且
$${\nabla }_{d}L(d^{\ast },\lambda )=(B+\lambda I)d^{\ast }+g=0.$$

此外还有互补条件
$$\frac{\lambda }{2}({\Delta }^2-\Vert d^{\ast }{\Vert }^2)=0,$$

以上两式整理后就是 (6.6.5a) 式和 (6.6.5b) 式. 为了证明 (6.6.5c) 式，任取 $d$ 满足 $\Vert d\Vert =\Delta$, 根据最优性，有
$$m(d)⩾m(d^{\ast })=m(d^{\ast })+\frac{\lambda }{2}(\Vert d^{\ast }{\Vert }^2-\Vert d{\Vert }^2).$$

利用 (6.6.5a) 式消去 $g$, 代入上式整理可知
$$(d-d^{\ast }{)}^{\mathrm{T}}(B+\lambda I)(d-d^{\ast })⩾0,$$

由任意性可知 $B+\lambda I$ 半正定.

‘充分性’
定义辅助函数
$$\hat{m}(d)=f+g^{\mathrm{T}}d+\frac{1}{2}{d}^{\mathrm{T}}(B+\lambda I)d=m(d)+\frac{\lambda }{2}{d}^{\mathrm{T}}d,$$

由条件 (6.6.5c) 可知 $\hat{m}(d)$ 关于 $d$ 是凸函数. 根据条件 (6.6.5a), $d^{\ast }$ 满足凸函数一阶最优性条件，结合定理 5.5 可推出 $d^{\ast }$ 是 $\hat{m}(d)$ 的全局极小值点, 进而对任意可行解 $d$, 我们有
$$m(d)⩾m(d^{\ast })+\frac{\lambda }{2}(\Vert d^{\ast }{\Vert }^2-\Vert d{\Vert }^2).$$

由互补条件 (6.6.5b) 可知 $\lambda ({\Delta }^2-\Vert d^{\ast }{\Vert }^2)=0$, 代入上式消去 $\Vert d^{\ast }{\Vert }^2$ 得
$$m(d)⩾m(d^{\ast })+\frac{\lambda }{2}({\Delta }^2-\Vert d{\Vert }^2)⩾m(d^{\ast }).$$

定理 5.5 假设 $f$ 是适当且凸的函数，则 $x^{\ast }$ 为问题 $\underset{x\in {\mathcal{R}}^n}{\min }f(x)$ 的一个全局极小点当且仅当 $0\in \partial f(x^{\ast })$.

定理 6.12 提供了问题维数 $n$ 较小时寻找 $d^{\ast }$ 的一个方法. 根据 (6.6.5a) 式，最优解是以 $\lambda$ 为参数的一族向量. 我们定义
$$d(\lambda )=-(B+\lambda I{)}^{-1}g(\mathrm{6.6.6})$$

则只需要寻找合适的 $\lambda$ 使得 (6.6.5b) 式和 (6.6.5c) 式成立即可．根据互补条件 (6.6.5b), 当 $\lambda >0$ 时必有 $\Vert d(\lambda )\Vert =\Delta$；根据半正定条件 (6.6.5c), $\lambda$ 须大于等于 $B$ 的最小特征值的相反数.
现在研究 $\Vert d(\lambda )\Vert$ 随 $\lambda$ 变化的性质. 设 $B$ 有特征值分解 $B=Q\Lambda Q^{\mathrm{T}}$, 其中 $Q=[q_1,q_2,\cdots ,q_n]$ 是正交矩阵, $\Lambda= $Diag$( \lambda_1, \lambda_2, \cdots , \lambda_n) $ 是对角矩阵，${\lambda }_1⩽{\lambda }_2⩽\cdots ⩽{\lambda }_n$ 是 $B$ 的特征值. 为了方便，以下仅考虑 ${\lambda }_1⩽0$ 且 ${\lambda }_1$ 是单特征根的情形. $B+\lambda I$ 有特征值分解 $B+\lambda I=Q(\Lambda +\lambda I)Q^{\mathrm{T}}$. 对$\lambda >-{\lambda }_1⩾0$, 那么 $d(\lambda )$ 的表达式为：
$$d(\lambda )=-Q(\Lambda +\lambda I{)}^{-1}{Q}^{\mathrm{T}}g=-\sum _{i=1}^n\frac{q_i^{\mathrm{T}}g}{{\lambda }_i+\lambda }{q}_i(\mathrm{6.6.7})$$

这正是 $d(\lambda )$ 的正交分解，由正交性可容易求出

$$\Vert d(\lambda ){\Vert }^2=\sum _{i=1}^n\frac{(q_i^{\mathrm{T}}g{)}^2}{({\lambda }_i+\lambda {)}^2}(\mathrm{6.6.8})$$

根据 (6.6.8) 式可知当 $\lambda >-{\lambda }_1$ 且 $q_1^{\mathrm{T}}g\ne 0$ 时，$\Vert d(\lambda ){\Vert }^2$ 是关于 $\lambda$ 的严格减函数，且有

$$\underset{\lambda \to \infty }{\lim }\Vert d(\lambda )\Vert =0,\underset{\lambda \to -{\lambda }_1+}{\lim }\Vert d(\lambda )\Vert =+\infty .$$

根据连续函数介值定理，$\Vert d(\lambda )\Vert =\Delta$ 的解必存在且唯一. 一个典型的例子如图 6.16 所示.

根据上面的分析，寻找 ${\lambda }^{\ast }$ 已经转化为一个一元方程求根问题，可以使用牛顿法求解，在得到最优解 ${\lambda }^{\ast }$ 后，根据 (6.6.5a) 式即可求出迭代方向 $d^{\ast }$.

此外，在上面的分析中假定了 $q_1^{\mathrm{T}}g\ne 0$, 在实际中这个条件未必满足. 当 $q_1^{\mathrm{T}}g=0$ 时，(6.6.8) 式将没有和 ${\lambda }_1$ 相关的项. 此时未必存在 ${\lambda }^{\ast }>-{\lambda }_1$ 使得 $\Vert d({\lambda }^{\ast })\Vert =\Delta$ 成立. 记 $M=\underset{\lambda \to -{\lambda }_1+}{\lim }\Vert d(\lambda )\Vert$, 当 $M⩾\Delta$ 时，仍然可以根据介值定理得出 ${\lambda }^{\ast }(>-{\lambda }_1)$ 的存在性；而当 $M<\Delta$ 时，无法利用前面的分析求出 ${\lambda }^{\ast }$ 和 $d^{\ast }$, 此时信赖域子问题变得比较复杂.
实际上，$q_1^{\mathrm{T}}g=0$ 且 $M<\Delta$ 的情形被人们称为“困难情形(hard case)”. 此情形发生时，区间 $(-{\lambda }_1,+\infty )$ 中的点无法使得 (6.6.5b) 成立，而定理 6.12 的结果说明 ${\lambda }^{\ast }\in [-{\lambda }_1,+\infty )$, 因此必有 ${\lambda }^{\ast }=-{\lambda }_1$.

2. 截断共轭梯度法
下面再介绍一种信赖域子问题的求解方法. 既然信赖域子问题的解不易求出，则求出其近似解. Steihaug 在 1983 年对共轭梯度法进行了改造，使其成为能求解子问题的算法. 此算法能够应用在大规模问题中，是一种非常有效的信赖域子问题的求解方法. 由于子问题 (6.6.2) 和一般的二次极小化问题相差一个约束，如果先不考虑其中的约束 $\Vert d\Vert ⩽\Delta$ 而直接使用共轭梯度法求解，在迭代过程中应该能找到合适的迭代点作为信赖域子问题的近似解. 这就是截断共轭梯度法的基本思想.
为了介绍截断共轭梯度法，先简要回顾一下标准共轭梯度法的迭代过程. 对于二次极小化问题
$$\underset{s}{\min }q(s)\stackrel{\mathrm{def}}{=}{g}^{\mathrm{T}}s+\frac{1}{2}{s}^{\mathrm{T}}Bs$$

给定初值 $s^0=0,r^0=g,p^0=-g$, 共轭梯度法的迭代过程为
$$\begin{array}{rl}{\alpha }_{k+1}& =\frac{\Vert r^k{\Vert }^2}{(p^k{)}^{\text{T}}B{p}^k}\\ s^{k+1}& =s^k+{\alpha }_k{p}^k\\ r^{k+1}& =r^k+{\alpha }_{k}B{p}^k\\ {\beta }_k& =\frac{\Vert r^{k+1}{\Vert }^2}{\Vert r^k{\Vert }^2}\\ p^{k+1}& =-r^{k+1}+{\beta }_k{p}^k\end{array}$$

其中迭代序列 { s k } \{s^k\} {sk} 最终的输出即为二次极小化问题的解，算法的终止准则是判断 $\Vert r^k\Vert$ 是否足够小.

截断共轭梯度法则是给标准的共轭梯度法增加了两条终止准则，并对最后一步的迭代点 $s^k$ 进行修正来得到信赖域子问题的解. 考虑到矩阵 $B$ 不一定是正定矩阵，在迭代过程中可能会产生如下三种情况：
(1) $(p^k{)}^{\mathrm{T}}B{p}^k⩽0$, 即 $B$ 不是正定矩阵. 由于共轭梯度法不能处理非正定的线性方程组，遇到这种情况应该立即终止算法. 但根据这个条件也找到了一个负曲率方向，此时只需要沿着这个方向走到信赖域边界即可.
(2) $(p^k{)}^{\mathrm{T}}B{p}^k>0$ 但 $\Vert s^{k+1}\Vert ⩾\Delta$, 这表示若继续进行共轭梯度法迭代，则点 $s^{k+1}$ 将处于信赖域之外或边界上，此时必须马上停止迭代，并在 $s^k$ 和 $s^{k+1}$ 之间找一个近似解.
(3) $(p^k{)}^{\mathrm{T}}B{p}^k>0$ 且 $\Vert r^{k+1}\Vert$ 充分小，这表示若共轭梯度法成功收敛到信赖域内. 子问题 (6.6.2) 和不带约束的二次极小化问题是等价的.
从上述终止条件来看截断共轭梯度法仅仅产生了共轭梯度法的部分迭代点，这也是该方法名字的由来.

算法 6.9 给出截断共轭梯度法的迭代过程，其中的三个判断分别对应了之前叙述的三种情况. 在不加限制时，下一步迭代点总是为 $s^{k+1}=s^k+{\alpha }_k{p}^k$. 当情况 (1) 发生时，只需要沿着 $p^k$ 走到信赖域边界；当情况 (2) 发生时，由于 $\Vert s^k\Vert <\Delta ,\Vert s^k+{\alpha }_k{p}^k\Vert ⩾\Delta$, 由连续函数介值定理可得 $\tau$ 必定存在且处于区间 $(0,{\alpha }_k]$ 内.

# 算法 6.9 截断共轭梯度法(Steihaug-CG)
import numpy as np
# from scipy.optimize import newton
def Steihaug_CG(B, g, Delta, epsilon = 1e-6):
    s = np.zeros_like(g)
    r = g
    p = -g
    k = 0
    while np.linalg.norm(g) > epsilon:
        if np.dot(p, np.dot(B, p)) <= 0:
            # 使用牛顿法求解方程
            # tau= newton(np.linalg.norm(s + tau * p) - Delta, 0, args=(s, p, Delta))
            # 求根公式求解
            tau = np.max([0, (-2*np.dot(s, p) + np.sqrt(4*np.dot(s, p)**2 - 4*np.dot(p, p)*(np.dot(s, s) - Delta**2))) / (2*np.dot(p, p))])
            s = s + tau * p
            return s
        alpha_k = np.linalg.norm(r)**2 / np.dot(p, np.dot(B, p))
        s = s + alpha_k * p
        if np.linalg.norm(s) >= Delta:
            # 使用牛顿法求解方程
            # tau= newton(np.linalg.norm(s + tau * p) - Delta, 0, args=(s, p, Delta))
            # 求根公式求解
            tau = np.max([0, (-2*np.dot(s, p) + np.sqrt(4*np.dot(s, p)**2 - 4*np.dot(p, p)*(np.dot(s, s) - Delta**2))) / (2*np.dot(p, p))])
            s = s + tau * p
            return s
        r_next = r + alpha_k * np.dot(B, p)
        if np.linalg.norm(r) < epsilon * np.linalg.norm(g):
            return s
        beta_k = np.linalg.norm(r_next)**2 / np.linalg.norm(r)**2
        p = -r_next + beta_k * p
        r = r_next
        k += 1
    return s

def m_k(x, d):
    return f(x) + np.dot(grad_f(x), d) + 0.5 * np.dot(d, np.dot(hess_f(x), d))
def cg_trust_region_algorithm(x_0, Delta_0, Delta_max, eta, rho1_bar = 0.25, rho2_bar = 0.75, gamma1 = 0.25, gamma2 = 2, max_iter = 1000, tol = 1e-6):
    x = x_0
    Delta = Delta_0
    k = 0
    x_history = [x]
    while k < max_iter:
        dk = Steihaug_CG(hess_f(x), grad_f(x), Delta, epsilon = 1e-6)  # 计算迭代方向
        mk_0 = m_k(x, np.zeros_like(x)) 
        mk_dk = m_k(x, dk)  
        rho_k = (f(x) - f(x + dk)) / (mk_0 - mk_dk)  # 计算下降率

        if rho_k < rho1_bar:
            Delta = gamma1 * Delta
        elif rho_k > rho2_bar and np.linalg.norm(dk) == Delta:
            Delta = min(gamma2 * Delta, Delta_max)
        else:
            Delta = Delta

        if rho_k > eta:
            x = x + dk
        x_history.append(x)
        k += 1
        if np.linalg.norm(grad_f(x)) <= tol:
            break    
    return x_history, k

def f(x):
    return x[0]**2 + 10 * x[1]**2
def grad_f(x):
    return np.array([2*x[0], 20*x[1]])
def hess_f(x):
    return np.array([[2, 0], [0, 20]])

x_0 = np.array([-10.0, -1.0])  # 初始点
Delta_0 = 1.0  # 初始半径
Delta_max = 10.0  # 最大半径
eta = 0.1
# 调用信赖域算法
cg_trust_history, k = cg_trust_region_algorithm(x_0, Delta_0, Delta_max, eta)
print("Optimal solution:", cg_trust_history[-1])
'''Optimal solution: [ 0.00000000e+00 -3.33066907e-16]'''
print("迭代次数：", k)
'''迭代次数： 6'''
# 绘制等高线图
x = np.linspace(-12, 12, 100)
y = np.linspace(-4, 4, 100)
X, Y = np.meshgrid(x, y)
Z = X**2 + 10*Y**2 
plt.figure(figsize=(10, 5))
plt.contour(X, Y, Z, levels=np.logspace(-2, 3, 20), cmap='jet')  # levels 参数控制等高线的密集程度
plt.plot([x[0] for x in cg_trust_history], [x[1] for x in cg_trust_history], markersize = 4, marker='x', markerfacecolor='none', linewidth=1, color='green', label='cg_trust_region_algorithm')
plt.plot([x[0] for x in trust_history], [x[1] for x in trust_history], markersize = 4, marker='o', markerfacecolor='none', linewidth=1, color='orange', label='trust_region_algorithm')
plt.legend(loc='best', markerscale=1.5)  # 设置标签位置为最佳位置，调整图例的大小
plt.show() 

截断共轭梯度法的迭代序列 {$s^k$} 有非常好的性质，实际上我们可以证明如下定理：

定理 6.13 设 $q(s)$ 是任意外迭代步信赖域子问题的目标函数，令 { s j } \{s^j\} {sj} 是由截断共轭梯度算法产生的迭代序列，则在算法终止前 $q(s^j)$ 是严格单调递减的，即
$$q(s^{j+1})<q(s^j)(\mathrm{6.6.9})$$

并且 $\Vert s^j\Vert$ 是严格单调递增的，即
$$0=\Vert s^0\Vert <\Vert s^1\Vert <\cdots <\Vert s^{j+1}\Vert <\cdots ⩽\Delta (\mathrm{6.6.10})$$

证明
设迭代在第 $t$ 步终止. 根据算法 6.9，在终止条件之前，$(p^j{)}^{\mathrm{T}}B{p}^j>0,j<t$ 是一直成立的. 此时的算法就是共轭梯度法，而容易证明 (6.6.9) 式和 (6.6.10) 式．又 $q(s)$ 在点 $s^j$ 处的梯度是 $r^j$，根据共轭梯度迭代性质 $(r^j{)}^{\mathrm{T}}{p}^i=0,i<j$，所以
$$(r^j{)}^{\mathrm{T}}{p}^j=(r^j{)}^{\mathrm{T}}(-r^j+{\beta }_{j-1}{p}^{j-1})=-\Vert r^j{\Vert }^2<0,$$

即 $p^j$ 是下降方向. 而 ${\alpha }_j$ 的选取恰好为精确线搜索的步长，因此有 $q(s^{j+1})<$ $q(s^j).$ 此外由 $s^j$ 的定义，
$$s^j=\sum _{i=0}^{j-1}{\alpha }_i{p}^i,{\alpha }_i>0.$$

再根据共轭梯度法的性质，有
$$(p^j{)}^{\mathrm{T}}{s}^j=\sum _{i=0}^{j-1}{\alpha }_i(p^j{)}^{\mathrm{T}}{p}^i=\sum _{i=0}^{j-1}{\alpha }_i\frac{\Vert r^j{\Vert }^2}{\Vert r^i{\Vert }^2}\Vert p^i{\Vert }^2>0.$$

结合以上表达式可得
$$\Vert s^{j+1}{\Vert }^2=\Vert s^j+{\alpha }_j{p}^j{\Vert }^2=\Vert s^j{\Vert }^2+2{\alpha }_j(p^j{)}^{\mathrm{T}}{s}^j+{\alpha }_j^2\Vert p^j{\Vert }^2>\Vert s^j{\Vert }^2.$$

实际上，我们还可进一步说明截断共轭梯度算法的输出 $s$ 满足如下关系：
$$q(s)⩽q(s^t),\Vert s^t\Vert ⩽\Vert s\Vert ,$$

其中 $t$ 为算法终止时的迭代数. 这只需要分别讨论三种终止条件即可.
(1) 若 $(p^t{)}^{\mathrm{T}}B{p}^t⩽0$, 则 $p^t$ 是负曲率方向，沿着负曲率方向显然有 $q(s)⩽q(s^t)$. 注意到此时 $\Vert s\Vert =\Delta$, 因此有 $\Vert s^t\Vert ⩽\Vert s\Vert =\Delta$.
(2) 若 $(p^t{)}^{\mathrm{T}}B{p}^t>0$ 但 $\Vert s^{t+1}\Vert ⩾\Delta$, 根据最速下降法的性质，$q(s^t+\alpha p^t)$ 关于 $\alpha \in (0,{\alpha }_t]$ 单调下降，根据 $\tau$ 的取法显然有 $q(s)⩽q(s^t)$. 此时依然有 $\Vert s\Vert =\Delta$, 因此 $\Vert s^t\Vert ⩽\Vert s\Vert =\Delta$ 仍成立.
(3) 若 $(p^t{)}^{\mathrm{T}}B{p}^t>0$ 且 $\Vert r^{t+1}\Vert ⩽\epsilon \Vert r^0\Vert$, 此时算法就是共轭梯度法，结论自然成立.

定理 6.13 其实可以从以下的观点来看：记 $d(\Delta )$ 为信赖域子问题当信赖域半径取 $\Delta$ 时的解，在这里让 $\Delta$ 从 0 开始变化，则 $d(\Delta )$ 的轨迹就是一条单参数曲线. 算法 6.9 实际上是使用由共轭梯度法迭代点列 { s j } \{s^j\} {sj} 确定的折线来逼近这条曲线，并将该折线与信赖域交集中距信赖域中心的最远点作为信赖域子问题的近似解. 根据定理 6.13，由于目标函数的单调递减性，截断共轭梯度法可以保证子问题的求解精度自动满足信赖域方法的收敛条件，从而保证理论与数值上有更好的效果．

6.6.3 收敛性分析

1. 柯西点
为了估计求解每个信赖域子问题得到的函数值改善情况，引入柯西点的定义．

定义 6.6 (柯西点) 设 $m_k(d)$ 是 $f(x)$ 在点 $x=x^k$ 处的二阶近似，常数 ${\tau }_k$ 为如下优化问题的解：
$$\begin{array}{rl}\min & m_k(-\tau \nabla f(x^k))\\ s.t.& \Vert \tau \nabla f(x^k)\Vert ⩽{\Delta }_k,\tau ⩾0\end{array}$$

则称 $x_C^k\stackrel{\mathrm{def}}{=}{x}^k+d_C^k$ 为柯西点，其中 $d_C^k=-{\tau }_k\nabla f(x^k)$.

根据柯西点的定义，它实际上是对 $m_k(d)$ 进行了一次精确线搜索的梯度法，不过这个线搜索是考虑了信赖域约束的. 图 6.17 直观地解释了柯西点的含义.

实际上，给定 $m_k(d)$, 柯西点可以显式计算出来. 为了方便用 $g^k$ 表示 $\nabla f(x^k)$, 根据 ${\tau }_k$ 的定义，其表达式为

$${\tau }_k=\{\begin{array}{ll}\frac{{\Delta }_k}{\Vert g^k\Vert },& (g^k{)}^{\text{T}}{B}^k{g}^k⩽0\\ \min \left\{\frac{\Vert g^k{\Vert }^2}{(g^k{)}^{\text{T}}{B}^k{g}^k},\frac{{\Delta }_k}{\Vert g^k\Vert }\right\},& 其他\end{array}$$

以上的分析表明，柯西点是信赖域子问题 (6.6.2) 的一个可行点，但在实际计算中人们一般不会将柯西点作为下一步迭代点的近似. 这是由于求柯西点本质上为一个带截断步长的最速下降法，它并没有充分利用海瑟矩阵 $B^k$ 的信息．即便如此，人们还是可以将柯西点作为信赖域子问题算法的一个评判标准，即要求子问题算法产生的迭代点至少比柯西点要好．若迭代点取为柯西点，二次模型的目标函数值依然是下降的．实际上，有如下引理：

引理 6.3 (柯西点的下降量) 设 $d_C^k$ 为求解柯西点产生的下降方向，则
$$m_k(0)-m_k(d_C^k)⩾\frac{1}{2}\Vert g^k\Vert \min \left\{{\Delta }_k,\frac{\Vert g^k\Vert }{\Vert B^k{\Vert }_2}\right\}(\mathrm{6.6.11})$$

证明
当 $(g^k{)}^T{B}^k{g}^k\le 0$,
$$\begin{array}{rl}& m_k(0)-m_k\left(d_c^k\right)=f-m^k\left(-\frac{{\Delta }^k}{\Vert g^k\Vert }{g}^k\right)\\ & =(g^k{)}^T\frac{{\Delta }_k}{\Vert g^k\Vert }{g}^k-\frac{1}{2}\frac{{\Delta }_k}{\Vert g^k\Vert }(g^k{)}^T{B}^k\frac{{\Delta }_k}{\Vert g^k\Vert }{g}^k\\ & \ge \frac{{\Delta }_k}{\Vert g^k\Vert }(g^k{)}^T{g}^k={\Delta }_k\Vert g^k\Vert \ge \Vert g^k\Vert \min \left\{{\Delta }_k,\frac{\Vert g^k\Vert }{\Vert B^k\Vert }\right\}\end{array}$$

当 $(g^k{)}^T{B}^k{g}^k>0$ 且 $\frac{\Vert g^k{\Vert }^3}{{\Delta }_k(g^k{)}^T{B}^k{g}^k}\le 1$,
$$\begin{array}{rl}& m_k(0)-m_k(d_c^k)=f-m_k\left(-\frac{\Vert g_k{\Vert }^2}{(g^k{)}^T{B}^k{g}_k}{g}_k\right)\\ & =(g^k{)}^T\frac{\Vert g^k{\Vert }^2}{(g^k{)}^T{B}^k{g}^k}{g}^k-\frac{1}{2}\frac{\Vert g^k{\Vert }^2}{(g^k{)}^T{B}^k{g}_k}(g^k{)}^T{B}^k\frac{\Vert g^k{\Vert }^2}{(g^k{)}^T{B}^k{g}_k}{g}^k\\ & =\frac{\Vert g^k{\Vert }^4}{(g^k{)}^T{B}^k{g}^k}-\frac{1}{2}\frac{\Vert g^k{\Vert }^4}{(g^k{)}^T{B}^k{g}^k}=\frac{1}{2}\frac{\Vert g^k{\Vert }^4}{(g^k{)}^T{B}^k{g}^k}\\ & \ge \frac{1}{2}\frac{\Vert g^k{\Vert }^2}{\Vert B_k\Vert }\ge \frac{1}{2}\Vert g^k\Vert \min \left\{{\Delta }_k,\frac{\Vert g^k\Vert }{\Vert B^k\Vert }\right\}\end{array}$$

当 $(g^k{)}^T{B}^k{g}^k>0$ 且 $\frac{\Vert g^k{\Vert }^3}{{\Delta }_k(g^k{)}^T{B}^k{g}^k}>1$,
$$\begin{array}{rl}& m_k(0)-m_k(d_c^k)=f-m_k\left(-\frac{{\Delta }_k}{\Vert g^k\Vert }{g}^k\right)\\ & =(g^k{)}^T\frac{{\Delta }_k}{\Vert g^k\Vert }{g}^k-\frac{1}{2}\frac{{\Delta }_k}{\Vert g^k\Vert }(g^k{)}^T{B}^k\frac{{\Delta }_k}{\Vert g^k\Vert }{g}^k\\ & ={\Delta }_k\Vert g^k\Vert -\frac{1}{2}\frac{{\Delta }_k^2}{{\Vert g^k\Vert }^2}(g^k{)}^T{B}^k{g}^k={\Delta }_k\Vert g^k\Vert -\frac{1}{2}\left({\Delta }_k\Vert g^k\Vert \right)\frac{{\Delta }_k(g^k{)}^T{B}^k{g}^k}{{\Vert g^k\Vert }^3}\\ & \ge {\Delta }_k\Vert g^k\Vert -\frac{1}{2}{\Delta }_k\Vert g^k\Vert =\frac{1}{2}{\Delta }_k\Vert g^k\Vert \ge \frac{1}{2}\Vert g^k\Vert \min \left\{{\Delta }_k,\frac{\Vert g^k\Vert }{\Vert B^k\Vert }\right\}\end{array}$$

而前一小节介绍的子问题求解方法 (如迭代法，截断共轭梯度，Dogleg 方法) 得到的迭代方向 $d^k$ 均满足
$$m_k(0)-m_k(d^k)⩾c_2(m_k(0)-m_k(d_C^k))$$

这也就意味着估计式

$$m_k(0)-m_k(d^k)⩾\frac{1}{2}{c}_2\Vert g^k\Vert \min \left\{{\Delta }_k,\frac{\Vert g^k\Vert }{\Vert B^k{\Vert }_2}\right\}(\mathrm{6.6.12})$$

在很多情况下都会成立。

2. Dogleg 方法
对柯西点的改进，首推 Dogleg 方法。
考虑 $B_k$ 正定的情况，当 $m_k(d)$ 的极小值位于信赖域内时，步长就是牛顿法的解，我们记作 $d^B=-{B^k}^{-1}\nabla f(x^k)$ 。当这个点是信赖域子问题的可行解时，显然 $d^{\ast }=d^B$
Dogleg 方法提供了一种良好的近似，使得我们可以轻松计算出针对任意信赖域半径的最佳步长。下图展示了 Dogleg 是如何对真实的最佳步长进行近似的。

图中 “Optimal trajectory" 即为真实的最佳步长，它是一条曲线，曲线方程我们无从得知。它的特点是，当信赖域半径很小时，比较接近于柯西点方向，当信赖域半径比较大时，比较接近于无约束的柯西点与 $d^B$ 的连线方向。根据这个规律，我们只需要计算出无约束的 $d^U$ 和 $d^B$, 就能得到两个线段，这就是近似的最佳步长所在的线段。我们把这条折线叫做 Dogleg，因为它看起来像一条狗腿。
在到达 $d^B$ 之前，沿着 Dogleg 前进，步长的大小是单调递增的，同时 $m_k(d)$ 是单调下降的。因此信赖域边界与 Dogleg 至多有一个交点。在 Dogleg 方法的每次送代中，只需要计算给定信赖域半径 $\Delta$ 对应的最佳步长，然后按照式(6.6.3)判断是否满足要求，是则进入下一个迭代，否则缩小 $\Delta$ 再次计算。
子问题: $\min m_k(d),s.t.\Vert d\Vert \le {\Delta }_k$
DogLeg 方法求解 $d_k$.
$$d^B=-B^{-1}g.d^U=-\frac{g^{T}g}{g^{T}Bg}g.$$

$$d_k=\{\begin{array}{ll}\tau d^U& 0\le \tau \le 1\\ d^U+(\tau -1)(d^B-d^U)& 1\le \tau \le 2\end{array}$$

其中 $\tau \in [0,2]$

Dogleg 法的基本思想是在当前点处，通过计算最速下降方向和牛顿方向，然后根据这两个方向的信息来选择一个合适的步长和方向进行迭代。具体来说，Dogleg 法会计算一个试探点，然后根据试探点的情况来更新当前点，并继续迭代直到满足停止条件。

3. 全局收敛性
回顾信赖域算法 6.8，我们引入了一个参数 $\eta$ 来确定是否应该更新迭代点．这分为两种情况：当 $\eta$ = 0 时，只要原目标函数有下降量就接受信赖域迭代步的更新；当 $\eta$ > 0 时，只有当改善量 ${\rho }_k$ 达到一定程度时再进行更新．在这两种情况下得到的收敛性结果是不同的，我们分别介绍这两种结果.

在 $\eta$ = 0 的条件下有如下收敛性定理：

定理 6.14 (全局收敛性 1) 设近似海瑟矩阵 $B^k$ 有界，即 $\Vert B^k{\Vert }_2⩽M,\forall k,f(x)$ 在下水平集 L = { x ∣ f ( x ) ⩽ f ( x 0 ) } \mathcal{L}=\{x\mid f(x)\leqslant f(x^0)\} L={x∣f(x)⩽f(x0)} 上有下界，且 $\nabla f(x)$ 在 $\mathcal{L}$ 的一个开邻域内利普希茨连续. 若 $d^k$ 为信赖域子问题的近似解且满足 (6.6.12) 式，算法 6.8 选取参数 $\eta =0$, 则
$$\underset{k\to \infty }{\mathrm{liminf}}\Vert \nabla f(x^k)\Vert =0$$

即 $x^k$ 的聚点中包含稳定点.

在数学中，“聚点”通常指的是集合中的一点，该点周围存在无限多的集合元素。具体来说：
对于集合 $A$ 中的点 $p$, 如果对于任意给定的 $p$ 的邻域 (包含 $p$ 的某个开集), 都存在集合 $A$ 中与 $p$ 不同的点，则称 $p$ 是 $A$ 的聚点。
如果一个点是集合的聚点，那么这个点可以被集合中的序列中的点无限接近。
在拓扑学和实分析中，聚点是一个重要的概念，用于描述集合中点的密集程度。

定理 6.14 表明若无条件接受信赖域子问题的更新，则信赖域算法仅仅有子序列的收敛性，迭代点序列本身不一定收敛．下面的定理则说明选取 $\eta >0$ 可以改善收敛性结果．

定理 6.15 (全局收敛性 2) 在定理 6.14 的条件下，若信赖域算法选取参数 $\eta >0$, 且信赖域子问题近似解 $d^k$ 满足 (6.6.12) 式，则
$$\underset{k\to \infty }{\lim }\Vert \nabla f(x^k)\Vert =0.$$

和牛顿类算法不同，信赖域算法具有全局收敛性，因此它对迭代初值选取的要求比较弱. 而牛顿法的收敛性极大地依赖初值的选取.

3. 局部收敛性
在构造信赖域子问题时利用了 $f(x)$ 的二阶信息，它在最优点附近应该具有牛顿法的性质. 特别地，当近似矩阵 $B^k$ 取为海瑟矩阵 ${\nabla }^{2}f(x^k)$ 时，根据信赖域子问题的更新方式，二次模型 $m_k(d)$ 将会越来越逼近原函数 $f(x)$, 最终信赖域约束 $\Vert d\Vert ⩽{\Delta }_k$ 将会失效. 此时信赖域方法将会和牛顿法十分接近，而根据定理 6.6, 牛顿法有 Q-二次收敛的性质，这个性质很自然地会继承到信赖域算法上.
和牛顿法不同的是，在信赖域迭代算法中，我们并不知道迭代点列是否已经接近最优点. 即使信赖域半径约束已经不起作用，如果 $B^k$ 没有取为海瑟矩阵 ${\nabla }^{2}f(x^k)$ 或者信赖域子问题没有精确求解，$d^k$ 一般也不会等于牛顿方向 $d_N^k$, 但是它们的误差往往是越来越小的. 我们有如下定理：

定理 6.16 设 $f(x)$ 在最优点 $x=x^{\ast }$ 的一个邻域内二阶连续可微，且 $\nabla f(x)$ 利普希茨连续，在最优点 $x^{\ast }$ 处二阶充分条件成立，即 ${\nabla }^{2}f(x^{\ast })\succ 0$. 若迭代点列 { x k } \{x^k\} {xk} 收敛到 $x^{\ast }$, 且在迭代中选取 $B^k$ 为海瑟矩阵 ${\nabla }^{2}f(x^k)$, 则对充分大的 $k$, 任意满足 (6.6.12) 式的信赖域子问题算法产生的迭代方向 $d^k$ 均满足
$$\Vert d^k-d_N^k\Vert =o(\Vert d_N^k\Vert )(\mathrm{6.6.13})$$

其中 $d_N^k$ 为第 $k$ 步迭代的牛顿方向且满足假设 $\Vert d_N^k\Vert ⩽\frac{{\Delta }_k}{2}$.

定理 6.16 说明若信赖域算法收敛，则当 $k$ 充分大时，信赖域半径的约束终将失效，且算法产生的迭代方向将会越来越接近牛顿方向.

根据定理 6.16 很容易得到收敛速度的估计.

推论 6.3 (信赖域算法的局部收敛速度) 在定理 6.16 的条件下，信赖域算法产生的迭代序列 {$x^k$} 具有 Q-超线性收敛速度．
证明
根据牛顿法收敛性定理，对牛顿方向，$$\Vert x^k-d_N^k-x^{\ast }\Vert =\mathcal{O}(\Vert x^k-x^{\ast }{\Vert }^2),$$

因此得到估计 $\Vert d_N^k\Vert =\mathcal{O}(\Vert x^k-x^{\ast }{\Vert }^2)$, 又根据牛顿法的局部收敛性定理，$$\Vert x^k+d^k-x^{\ast }\Vert ⩽\Vert x^k+d_N^k-x^{\ast }\Vert +\Vert d^k-d_N^k\Vert =o(\Vert x^k-x^{\ast }{\Vert }^2)+o(\Vert d_N^k\Vert )=o(\Vert x^k-x^{\ast }\Vert ).$$

这说明信赖域算法是Q-超线性收敛的．

可以看出，若在迭代后期 $d^k=d_N^k$ 能得到满足，则信赖域算法是 Q-二次收敛的. 很多算法都会有这样的性质，例如前面提到的截断共轭梯度法. 因此在实际应用中，截断共轭梯度法是最常用的信赖域子问题的求解方法，使用此方法能够同时兼顾全局收敛性和局部 Q-二次收敛性.

参考教材：《最优化：建模、算法与理论》