非线性最优化方法概述

非线性最优化算法

算法一般性流程

在这里插入图片描述
在非线性优化理论中,迭代算法通常遵循 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk 的形式,其中:

  • x k x_k xk 是当前迭代点
  • d k d_k dk 是搜索方向
  • α k α_k αk 是移动步长

不同的优化算法采用不同的策略来确定搜索方向 d k d_k dk。一个好的搜索方向应该至少是一个下降方向,即在该方向上目标函数值能够减小。如果 ∇ f ( x k ) T d < 0 \nabla f(x_k)^T d < 0 f(xk)Td<0,则方向 d d d 是函数 f ( x ) f(x) f(x) 在点 x k x_k xk 处的下降方向。

搜索方向的选择

以下是一些主要的、不同的计算搜索方向的方法:

最速下降法

  • 思想: \textbf{思想:} 思想: 沿负梯度方向搜索,因为负梯度方向是函数值在局部下降最快的方向。
  • 计算: \textbf{计算:} 计算: d k = − ∇ f ( x k ) d_k = -\nabla f(x_k) dk=f(xk)
  • 特点: \textbf{特点:} 特点:
    • 简单直观,易于实现
    • 保证是下降方向(梯度非零)
    • 收敛速度慢,尤其是在目标函数等值线呈狭长山谷状时,会出现“锯齿形”现象
    • 属于一阶方法(只使用一阶导数)
  • 算法收敛特点: \textbf{算法收敛特点:} 算法收敛特点:
    最速下降法在每一步选择的是当前点局部下降最快的方向,即负梯度方向,然而,这个局部最快的下降方向,并不一定是通向全局最优解的最短或最直接的路径。 其算法收敛路径如下图所示:
    在这里插入图片描述
    在当前迭代点 x k x_k xk沿着负梯度方向 d k = − ∇ f ( x k ) d_k = -\nabla f(x_k) dk=f(xk) 进行线搜索时,那么在新的迭代点 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk 处,目标函数 φ ( α ) = f ( x k + α d k ) \varphi(α) = f(x_k + α d_k) φ(α)=f(xk+αdk) 关于 α α α 的导数应为零,即
    φ ′ ( α k ) = ∇ f ( x k + α k d k ) T d k = ∇ f ( x k + 1 ) T d k = 0 \varphi'(α_k) = \nabla f(x_k + α_k d_k)^T d_k = \nabla f(x_{k+1})^T d_k = 0 φ(αk)=f(xk+αkdk)Tdk=f(xk+1)Tdk=0 d k = − ∇ f ( x k ) d_k = -\nabla f(x_k) dk=f(xk) 代入,得到:
    ∇ f ( x k + 1 ) T ( − ∇ f ( x k ) ) = 0 \nabla f(x_{k+1})^T (-\nabla f(x_k)) = 0 f(xk+1)T(f(xk))=0
    ∇ f ( x k + 1 ) T ∇ f ( x k ) = 0 \nabla f(x_{k+1})^T \nabla f(x_k) = 0 f(xk+1)Tf(xk)=0 这意味着连续两次迭代的梯度方向是正交的(基于精确线搜索)。在实际求解过程中,一般使用非精确线搜索寻找 α k α_k αk ,这些方法不追求 ∇ f ( x k + 1 ) T d k = 0 \nabla f(x_{k+1})^T d_k = 0 f(xk+1)Tdk=0,因此连续的搜索方向通常不是严格正交的。

牛顿法

  • 思想: \textbf{思想:} 思想: 在当前迭代点 x k x_k xk 处,利用二阶泰勒展开近似目标函数 f ( x ) f(x) f(x),然后取该近似函数的极小点(或驻点)作为下一个迭代点 x k + 1 x_{k+1} xk+1 的搜索方向 d k d_k dk和步长 α k α_k αk( α k α_k αk通常取1)。

  • 计算: \textbf{计算:} 计算: d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) d_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k) dk=[2f(xk)]1f(xk),其计算步骤如下:
    首先将目标函数 f ( x ) f(x) f(x) 在点 x k x_k xk 附近二阶泰勒展开为:
    f ( x k + s ) ≈ q ( s ) = f ( x k ) + ∇ f ( x k ) T s + 1 2 s T ∇ 2 f ( x k ) s f(x_k + s) \approx q(s) = f(x_k) + \nabla f(x_k)^Ts + \frac{1}{2}s^T\nabla ^2 f(x_k)s f(xk+s)q(s)=f(xk)+f(xk)Ts+21sT2f(xk)s 其中:

    • s = x − x k s = x - x_k s=xxk 是从 x k x_k xk 出发的位移(即搜索步);
    • ∇ f ( x k ) \nabla f(x_k) f(xk) f ( x ) f(x) f(x) x k x_k xk 处的梯度向量;
    • ∇ 2 f ( x k ) \nabla ^2 f(x_k) 2f(xk) f ( x ) f(x) f(x) x k x_k xk 处的Hessian矩阵,记为 H k H_k Hk.

    将二次模型 q ( s ) q(s) q(s) 取极值时的解 s s s 作为搜索方向,对 q ( s ) q(s) q(s) 关于 s s s 求导并令其等于零:
    q ′ ( s ) = ∇ f ( x k ) + ∇ 2 f ( x k ) s = 0 q'(s) = \nabla f(x_k) + \nabla ^2 f(x_k)s = 0 q(s)=f(xk)+2f(xk)s=0 如果Hessian矩阵 ∇ 2 f ( x k ) \nabla ^2f(x_k) 2f(xk) 是非奇异的(可逆的),则解得:
    s = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) s = -[\nabla ^2 f(x_k)]^{-1} \nabla f(x_k) s=[2f(xk)]1f(xk) s s s 被称为牛顿方向,记为 d k d_k dk
    d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) d_k= -[\nabla ^2f(x_k)]⁻¹ \nabla f(x_k) dk=[2f(xk)]1f(xk)

  • 特点: \textbf{特点:} 特点:

    • 在最优解附近具有二次收敛速度,收敛速度非常快
    • 需要计算并存储Hessian矩阵及其逆矩阵,计算量大,特别是对于高维问题
    • Hessian矩阵可能非正定,此时牛顿方向可能不是下降方向,甚至算法可能失效。需要修正(如阻尼牛顿法、修正Cholesky分解等)
    • 对初始点要求较高,远离最优解时可能不收敛
    • 属于二阶方法(使用一阶和二阶导数)
  • 算法收敛特征: \textbf{算法收敛特征:} 算法收敛特征:

    • “一步最优”情况
      当目标函数为二次函数,其Hessian矩阵正定且可逆,牛顿法可以一步到达最优点,原因在于该二次函数的泰勒展开是精确的。 假设目标函数 f ( x ) f(x) f(x) 是一个二次函数
      f ( x ) = c + b T x + 1 2 x T A x f(x) = c + b^T x + \frac{1}{2}x^T Ax f(x)=c+bTx+21xTAx 其中 A A A 是对称的Hessian矩阵, b b b 是一个常向量, c c c 是一个常数,则 ∇ f ( x ) = b + A x \nabla f(x) = b + Ax f(x)=b+Ax ∇ 2 f ( x ) = A \nabla ^2 f(x) = A 2f(x)=A f ( x ) f(x) f(x)在任意点 x k x_k xk 的二阶泰勒展开为:
      q ( s ) = f ( x k ) + ∇ f ( x k ) T s + 1 2 s T ∇ 2 f ( x k ) s q(s) = f(x_k) + \nabla f(x_k)^Ts + \frac{1}{2}s^T\nabla^2 f(x_k)s q(s)=f(xk)+f(xk)Ts+21sT2f(xk)s 代入 ∇ f ( x k ) \nabla f(x_k) f(xk) ∇ 2 f ( x k ) \nabla ^2 f(x_k) 2f(xk) 后得到:
      q ( s ) = f ( x k ) + ( b + A x k ) T s + 1 2 s T A s q(s) = f(x_k) + (b + Ax_k)^Ts + \frac{1}{2}s^TAs q(s)=f(xk)+(b+Axk)Ts+21sTAs x = s + x k x = s + x_k x=s+xk代入 f ( x ) f(x) f(x)中,得到1
      f ( x k + s ) = c + b T ( x k + s ) + 1 2 ( x k + s ) T A ( x k + s ) = c + b T x k + b T s + 1 2 ( x k T A + s T A ) ( x k + s ) = c + b T x k + b T s + 1 2 ( x k T A x k + s T A x k + x k T A s + s T A s ) = c + b T x k + b T s + 1 2 x k T A x k + x k T A s + 1 2 s T A s = [ c + b T x k + 1 2 x k T A x k ] + ( b + A x k ) T s + 1 2 s T A s − − − ( f ( x k ) 代入 q ( s ) 后的形式 ) = q ( s ) \begin{align*} f(x_k+s) & = c + b^T (x_k+s) + \frac{1}{2}(x_k+s)^T A(x_k+s)\\ & = c+b^Tx_k +b^T s +\frac{1}{2}(x_k^T A+s^T A)(x_k+s)\\ & = c+b^Tx_k +b^T s +\frac{1}{2}(x_k^T A x_k +s^T A x_k + x_k^T A s + s^T A s)\\ & = c+b^Tx_k +b^T s +\frac{1}{2}x_k^T A x_k + x_k^T A s + \frac{1}{2}s^T A s\\ & = [c+b^Tx_k+\frac{1}{2}x_k^T A x_k] + (b+Ax_k)^T s+ \frac{1}{2}s^T A s ---(f(x_k)代入q(s)后的形式)\\ & = q(s) \end{align*} f(xk+s)=c+bT(xk+s)+21(xk+s)TA(xk+s)=c+bTxk+bTs+21(xkTA+sTA)(xk+s)=c+bTxk+bTs+21(xkTAxk+sTAxk+xkTAs+sTAs)=c+bTxk+bTs+21xkTAxk+xkTAs+21sTAs=[c+bTxk+21xkTAxk]+(b+Axk)Ts+21sTAs(f(xk)代入q(s)后的形式)=q(s) 因此函数 f ( x ) f(x) f(x)在任意点 x k x_k xk 处的二阶泰勒展开 q ( s ) q(s) q(s)完全等于其本身。那么 q ( s ) q(s) q(s) 的驻点(通过 ∇ q ( s k ) = 0 \nabla q(s_k) = 0 q(sk)=0)所对应的点 x k + 1 = x k + s k x_{k+1}=x_k + s_k xk+1=xk+sk,就是原始函数 f ( x ) f(x) f(x) 的驻点。当 f ( x ) f(x) f(x) 是一个凸二次函数( A A A 是正定的),它有唯一的全局最小值。其求解过程如下:
      ∇ q ( s k ) = ∇ f ( x k ) + A s k = 0 \nabla q(s_k) = \nabla f(x_k) + As_k = 0 q(sk)=f(xk)+Ask=0 得到 s k = − A − 1 ∇ f ( x k ) s_k = -A^{-1} \nabla f(x_k) sk=A1f(xk),则 x k + 1 = x k − A − 1 ∇ f ( x k ) x_{k+1}=x_k -A^{-1}\nabla f(x_k) xk+1=xkA1f(xk) f ( x ) f(x) f(x) x k + 1 x_{k+1} xk+1处的梯度为:
      ∇ f ( x k + 1 ) = b + A x k + 1 = b + A ( x k − A − 1 ∇ f ( x k ) ) = b + A x k − A A − 1 ∇ f ( x k ) = ∇ f ( x k ) − ∇ f ( x k ) = 0 \begin{align*} \nabla f(x_{k+1}) & = b + A x_{k+1}\\ & = b + A(x_k - A^{-1} \nabla f(x_k))\\ & = b + Ax_k - AA^{-1} \nabla f(x_k)\\ & = \nabla f(x_k) - \nabla f(x_k)\\ & = 0 \end{align*} f(xk+1)=b+Axk+1=b+A(xkA1f(xk))=b+AxkAA1f(xk)=f(xk)f(xk)=0 所以,在 x k + 1 x_{k+1} xk+1 处,梯度为零。这意味着 x k + 1 x_{k+1} xk+1 就是该二次函数的全局最小值点。其收敛特征如下图所示。
      在这里插入图片描述
    • 需要修正情况
      当牛顿法面临Hessian矩阵非正定、计算成本过高或全局收敛性差等问题时,需要对其进行修正。主要的修正方向包括:
      • 阻尼牛顿法
        问题:
        在远离最优解或Hessian非正定时,牛顿法的步长 α k = 1 α_k=1 αk=1 可能过大,将导致牛顿法难以收敛。
        修正思想:
        在计算出牛顿方向 d k = − [ ∇ 2 f ( x k ) ] − 1 ∇ f ( x k ) d_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k) dk=[2f(xk)]1f(xk) 后,不直接取单位步长,而是引入一个步长因子 α k > 0 α_k > 0 αk>0,并通过线搜索来确定 α k α_k αk 的值,即 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk,则 α k = a r g m i n α > 0 f ( x k + α d k ) α_k = argmin_{α>0} f(x_k + α d_k) αk=argminα>0f(xk+αdk) ( α k α_k αk通过满足Armijo或Wolfe条件的非精确线搜索寻找)。当迭代点接近最优解时, α k α_k αk 通常会自动趋向于1,从而保留牛顿法的快速局部收敛性。然而,由于Hessian矩阵非正定,牛顿方向 d k d_k dk 可能不是一个下降方向,即存在 ∇ f ( x k ) T d k ≥ 0 \nabla f(x_k)^T d_k \ge 0 f(xk)Tdk0,那么无论 α k > 0 α_k>0 αk>0取何值都不可能使目标函数值下降。因此,阻尼牛顿法通常需要先对Hessian矩阵进行修正,以确保 d k d_k dk 是一个下降方向。
        Hessian矩阵的修正方法
        ∇ 2 f ( x k ) ( 或 H k ) \nabla ^2 f(x_k)(或H_k) 2f(xk)(Hk) 非正定时,对其进行修正,得到一个正定矩阵 B k B_k Bk,然后用 B k B_k Bk 代替 H k H_k Hk 来计算搜索方向: d k = − B k ∇ f ( x k ) d_k = -B_k \nabla f(x_k) dk=Bkf(xk)
        1. 对角加载
          方法: 给Hessian矩阵的对角线元素加上一个正数 μ k \mu_k μk 使其得到的矩阵正定,即 B k = H k + μ k I B_k = H_k + \mu_k I Bk=Hk+μkI ,其中 I I I 为单位矩阵。
          特性:

          • μ k \mu_k μk 很小时, B k ≈ H k B_k \approx H_k BkHk,搜索方向接近牛顿方向。
          • μ k \mu_k μk 很大时, B k B_k Bk μ k I \mu_k I μkI 主导( B k ≈ μ k I B_k \approx \mu_k I BkμkI),则 d k ≈ − 1 μ k ∇ f ( x k ) d_k \approx -\frac{1}{\mu_k} \nabla f(x_k) dkμk1f(xk),此时搜索方向接近最速下降方向。
        2. 修正Cholesky分解
          方法: H k H_k Hk进行Cholesky分解 ,即 H k = L D L T H_k = LDL^T Hk=LDLT,其中 L L L 是单位下三角矩阵, D = ( d 1 , d 2 , ⋯   , d n ) D=(d_1,d_2,\cdots,d_n) D=(d1,d2,,dn) 是对角矩阵。在计算 d i d_i di时,如果 d i < 0 d_i<0 di<0,则令 d i ~ = d i + σ \tilde{d_i} = d_i+\sigma di~=di+σ,使得 d i ~ > 0 \tilde{d_i}>0 di~>0,同时,单位下三角矩阵 L L L中的元素发生相应的修改。对所有元素修正完成后得到的新的单位下三角矩阵 L ~ \tilde{L} L~和对角矩阵 D ~ \tilde{D} D~。最后, B k = L ~ D ~ L T ~ B_k=\tilde{L} \tilde{D} \tilde{L^T} Bk=L~D~LT~

拟牛顿法

  • 思想: 试图兼顾梯度下降法的低计算成本和牛顿法的快速收敛性。通过直接构造一个矩阵 H k H_k Hk 来近似Hessian矩阵的逆 ∇ 2 f ( x k ) − 1 \nabla ^2 f(x_k)^{-1} 2f(xk)1。降低了每一步的计算成本,同时保持了超线性收敛的特性。

  • 计算: d k = − H k ∇ f ( x k ) d_k = -H_k \nabla f(x_k) dk=Hkf(xk)。拟牛顿法的一般性迭代框架为:

    1. 给定初始点 x 0 x_0 x0,初始Hessian逆近似 H 0 H_0 H0 (通常为单位矩阵 I I I);
    2. 计算梯度 g k = ∇ f ( x k ) g_k = \nabla f(x_k) gk=f(xk);
    3. 计算搜索方向 d k = − H k g k d_k = -H_k g_k dk=Hkgk;
    4. 通过线搜索确定步长 α k α_k αk,使得 f ( x k + α k d k ) f(x_k + α_k d_k) f(xk+αkdk) 充分下降;
    5. 更新迭代点 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk;
    6. 计算 s k = x k + 1 − x k s_k = x_{k+1} - x_k sk=xk+1xk y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) y_k = \nabla f(x_{k+1}) - \nabla f(x_k) yk=f(xk+1)f(xk);
    7. 使用 s k s_k sk y k y_k yk 以及当前的 H k H_k Hk 来更新得到 H k + 1 H_{k+1} Hk+1;
    8. 检查收敛条件,若不满足则返回步骤 2 2 2
  • 常用的 H k H_k Hk更新公式:

    • DFP公式:
      H k + 1 = H k + s k s k T s k T y k − H k y k y k T H k y k T H k y k H_{k+1} = H_k + \frac{s_k s_k^T}{s_k^T y_k} - \frac{H_k y_k y_k^T H_k}{y_k^T H_k y_k} Hk+1=Hk+skTykskskTykTHkykHkykykTHk
    • BFGS公式:
      H k + 1 = ( I − s k y k T y k T s k ) H k ( I − y k s k T y k T s k ) + s k s k T y k T s k H_{k+1} = (I - \frac{s_k y_k^T}{y_k^T s_k}) H_k (I - \frac{y_k s_k^T}{y_k^T s_k}) + \frac{s_k s_k^T}{y_k^T s_k} Hk+1=(IykTskskykT)Hk(IykTskykskT)+ykTskskskT
    • L-BFGS公式:
      适用于大规模问题,不存储完整的 H k H_k Hk 矩阵,而是存储最近几次迭代的梯度和变量变化信息来隐式地计算 H k ∇ f ( x k ) H_k \nabla f(x_k) Hkf(xk)
  • 特点:

    • 具有超线性收敛速度,通常比最速下降快得多,比牛顿法在全局收敛性上更鲁棒
    • 适用于Hessian矩阵非正定或奇异
    • 适用于中等规模问题(BFGS)或大规模问题(L-BFGS)
    • 属于(近似的)二阶方法

基于BFGS方法的搜索结果如下图所示。
在这里插入图片描述

共轭梯度法

  • 线性共轭梯度法
    • 思想: \textbf{思想:} 思想: 对于二次规划问题 φ ( x ) = 1 2 x T A x − b T x \varphi(x) = \frac{1}{2} x^T A x - b^T x φ(x)=21xTAxbTx A A A对称正定),迭代地生成一组关于矩阵 A A A共轭的搜索方向 { d 0 , d 1 , . . . , d k } \{d_0, d_1, ..., d_{k}\} {d0,d1,...,dk} ( A A A-共轭),则对所有的 i ≠ j i \neq j i=j,满足 d i T A d j = 0 d_i^T A d_j = 0 diTAdj=0。从初始点 x 0 x_0 x0 开始,沿着这些 A A A-共轭方向 d k d_k dk 进行精确线搜索,通过迭代公式 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk逐步逼近最优解 x ∗ x^* x。理论上,对于 n n n 维问题,在没有舍入误差的情况下,线性共轭梯度法最多经过 n n n 步迭代就能找到精确解。
    • 计算流程: \textbf{计算流程:} 计算流程:
      1. 初始化迭代点 x 0 x_0 x0,初始残差 r 0 = b − A x 0 r_0 = b - Ax_0 r0=bAx0,初始搜索方向 d 0 = r 0 d_0 = r_0 d0=r0;设置容差 σ \sigma σ
      2. 计算迭代步长 α k = r k T r k d k T A d k \alpha_k=\frac{r_k^T r_k}{d_k^T A d_k} αk=dkTAdkrkTrk;
      3. 更新新的迭代点 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk;
      4. 更新残差 r k + 1 = r k − α k A d k r_{k+1} = r_k - α_k A d_k rk+1=rkαkAdk 2;
      5. 检查是否收敛,如果 ∣ ∣ r k + 1 ∣ ∣ < σ ||r_{k+1}||< \sigma ∣∣rk+1∣∣<σ ,则停止迭代;
      6. 计算下一个共轭方向 β k = r k + 1 T r k + 1 r k T r k β_k = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k} βk=rkTrkrk+1Trk+1;
      7. 更新搜索方向 d k + 1 = r k + 1 + β k d k d_{k+1} = r_{k+1} + β_k d_k dk+1=rk+1+βkdk d k + 1 d_{k+1} dk+1与之前的所有搜索方向 d j d_j dj共轭, d k + 1 A d j = 0 d_{k+1} A d_j = 0 dk+1Adj=0);
      8. 转到步骤 2 2 2.
    • 算法特点: \textbf{算法特点:} 算法特点:
      • 存储高效:只需要存储少数几个 n n n 维向量 ( x k , r k , d k , A d k x_k, r_k, d_k, Ad_k xk,rk,dk,Adk),存储复杂度为 O ( n ) O(n) O(n)
      • 有限步收敛(理论上):对于 n n n 维问题,在精确算术下,最多 n n n 步收敛。
      • 单调下降:目标函数 φ ( x ) \varphi(x) φ(x) 的值在迭代过程中是单调下降的。

线性共轭梯度法的搜索结果如下图所示。
在这里插入图片描述
其中,
A = ( 3 2 2 6 ) , b = ( 2 8 ) , X 0 = ( − 3 − 2 ) A= \begin{pmatrix} 3 & 2 \\ 2 & 6 \end{pmatrix}, b = \begin{pmatrix} 2 \\ 8 \end{pmatrix}, X_0= \begin{pmatrix} -3 \\ -2 \end{pmatrix} A=(3226),b=(28),X0=(32)

  • 非线性共轭梯度法

    • 思想: \textbf{思想:} 思想:
      非线性共轭梯度法将线性共轭梯度法的思想推广到非二次目标函数。由于Hessian矩阵 ∇ 2 f ( x ) \nabla^2 f(x) 2f(x) 不再是常数 A A A,且通常难以计算或存储,非线性共轭梯度法避免了直接使用Hessian矩阵。与线性共轭梯度法类似,非线性共轭梯度法也生成一系列搜索方向 d k d_k dk,使得新的方向 d k + 1 d_{k+1} dk+1 是当前负梯度 − ∇ f ( x k + 1 ) -∇f(x_{k+1}) f(xk+1) 和前一个搜索方向 d k d_k dk 的线性组合: d k + 1 = − ∇ f ( x k + 1 ) + β k + 1 d k d_{k+1} = -∇f(x_{k+1}) + β_{k+1} d_k dk+1=f(xk+1)+βk+1dk
    • 计算流程: \textbf{计算流程:} 计算流程:
      1. 初始化迭代点 x 0 x_0 x0,初始梯度 g 0 = ∇ f ( x 0 ) g_0 = \nabla f(x_0) g0=f(x0),初始搜索方向 d 0 = − g 0 d_0 = -g_0 d0=g0;设置容差 σ \sigma σ
      2. 计算迭代步长 α k \alpha_k αk(使用非精确线搜索,Armijo/Wolfe条件);
      3. 更新新的迭代点 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk;
      4. 计算新梯度 g k + 1 = ∇ f ( x k + 1 ) g_{k+1} = \nabla f(x_{k+1}) gk+1=f(xk+1);
      5. 检查是否收敛,如果 ∣ ∣ g k + 1 ∣ ∣ < σ ||g_{k+1}||< \sigma ∣∣gk+1∣∣<σ ,则停止迭代;
      6. 计算 β k + 1 β_{k+1} βk+1,几种不同的方法如下:
        • Fletcher-Reeves (FR): β k + 1 = g k + 1 T g k + 1 g k T g k β_{k+1} = \frac{g_{k+1}^T g_{k+1}}{g_k^T g_k} βk+1=gkTgkgk+1Tgk+1
        • Polak-Ribière (PR): β k + 1 = g k + 1 T ( g k + 1 − g k ) g k T g k β_{k+1} = \frac{g_{k+1}^T (g_{k+1} - g_k)}{g_k^T g_k} βk+1=gkTgkgk+1T(gk+1gk)
        • Hestenes-Stiefel (HS): β k + 1 = g k + 1 T ( g k + 1 − g k ) ( g k + 1 − g k ) T d k β_{k+1} = \frac{g_{k+1}^T (g_{k+1} - g_k)}{(g_{k+1} - g_k)^T d_k} βk+1=(gk+1gk)Tdkgk+1T(gk+1gk)
        • Dai-Yuan (DY): β k + 1 = g k + 1 T g k + 1 ( g k + 1 − g k ) T d k β_{k+1} = \frac{g_{k+1}^T g_{k+1}}{(g_{k+1} - g_k)^T d_k} βk+1=(gk+1gk)Tdkgk+1Tgk+1
        • Polak-Ribière Plus (PRP): β k + 1 = m a x ( 0 , β k + 1 P R ) β_{k+1}= max(0, β_{k+1}^{PR}) βk+1=max(0,βk+1PR)
      7. 更新搜索方向 d k + 1 = − g k + 1 + β k + 1 d k d_{k+1} = -g_{k+1} + β_{k+1} d_k dk+1=gk+1+βk+1dk;
      8. 转到步骤 2 2 2.
  • 算法特点: \textbf{算法特点:} 算法特点:

    • 不需要存储Hessian矩阵,计算量小,适用于大规模问题
    • 收敛速度通常介于最速下降和拟牛顿法之间。
    • 对于非二次函数,需要周期性地重置(例如,每 n n n次迭代后将方向设为负梯度方向)。
    • 属于一阶方法(主要依赖梯度信息,历史信息被用来构造共轭性)。

基于Fletcher-Reeves和Polak-Ribière方法的非线性牛顿法的求解Rosenbrock函数结果分别如下所示。
在这里插入图片描述在这里插入图片描述

搜索步长的选择

在非线性优化中,一旦确定了搜索方向 d k d_k dk,下一步就是确定沿着这个方向走多远,即选择一个合适的步长 α k > 0 α_k > 0 αk>0,使得目标函数 f ( x ) f(x) f(x) 的值在 x k + 1 = x k + α k d k x_{k+1} = x_k + α_k d_k xk+1=xk+αkdk 处有足够的下降。选择步长的方法主要分为线搜索方法信赖域方法两大类,其中线搜索方法包括精确线搜索非精确线搜索

精确线搜索

  • 思想: \textbf{思想:} 思想:
    给定一个搜索方向 d k d_k dk ,沿着这个方向找到一个最优解 α k \alpha_k αk,使得目标函数 f ( x ) f(x) f(x) 在这个方向上达到最小值。也就是说,给定当前迭代点 x k x_k xk 和搜索方向 d k d_k dk,定义一个关于迭代步长 α α α 的一维函数 φ ( α ) = f ( x k + α ∗ d k ) \varphi(α) = f(x_k + α * d_k) φ(α)=f(xk+αdk) 精确线搜索的目标就是求解该一维优化问题,即
    α k = arg min ⁡ α > 0 φ ( α ) α_k = \argmin_{α > 0} \varphi (α) αk=α>0argminφ(α)

  • 计算方法: \textbf{计算方法:} 计算方法:
    基于函数 φ ( α ) \varphi(α) φ(α) 的性质以及可用的信息,求解 φ ( α ) = f ( x k + α ∗ d k ) \varphi(α) = f(x_k + α * d_k) φ(α)=f(xk+αdk)的方法可以分为解析方法和数值方法。

    • 解析方法:
      如果函数 f ( x ) f(x) f(x) 的形式比较简单,使得 φ ( α ) \varphi(α) φ(α) 的导数 φ ′ ( α ) \varphi'(α) φ(α) 容易计算,并且方程 φ ′ ( α ) = 0 \varphi'(α)=0 φ(α)=0 可以解析地求解,则 α k \alpha_k αk 往往可以求解出解析式。假设目标函数 f ( x ) f(x) f(x) 是一个凸二次函数:
      f ( x ) = 1 2 x T A x − b T x + c f(x) = \frac{1}{2}x^T A x - b^T x + c f(x)=21xTAxbTx+c 其中 A A A是对称正定矩阵,则 ∇ f ( x ) = A x − b \nabla f(x) = Ax - b f(x)=Axb。对 φ ( α ) \varphi(α) φ(α)求导:
      φ ′ ( α ) = ∇ f ( x k + α ∗ d k ) T d k = ( A ( x k + α ∗ d k ) − b ) T d k = ( A x k − b + α A d k ) T d k = ( ∇ f ( x k ) + α A d k ) T d k \begin{align*} \varphi' (\alpha) & = \nabla f(x_k + \alpha * d_k)^Td_k \\ & = (A(x_k + \alpha * d_k)-b)^Td_k \\ & = (Ax_k -b +\alpha A d_k)^T d_k \\ & = (\nabla f(x_k) + \alpha A d_k)^T d_k \end{align*} φ(α)=f(xk+αdk)Tdk=(A(xk+αdk)b)Tdk=(Axkb+αAdk)Tdk=(f(xk)+αAdk)Tdk φ ′ ( α ) = 0 \varphi'(α)=0 φ(α)=0,则:
      ∇ f ( x k ) T d k + α d k T A d k = 0 \begin{align*} \nabla f(x_k)^T d_k + \alpha d_k^T A d_k = 0 \end{align*} f(xk)Tdk+αdkTAdk=0 解得:
      α k = − ∇ f ( x k ) T d k d k T A d k α_k = -\frac{\nabla f(x_k)^T d_k}{d_k^T A d_k} αk=dkTAdkf(xk)Tdk 这正是最速下降法或共轭梯度法在求解二次规划问题时使用的精确迭代步长公式。对于一般的非线性函数 f ( x ) f(x) f(x),方程 ∇ f ( x k + α d k ) T d k = 0 \nabla f(x_k + α d_k)^T d_k = 0 f(xk+αdk)Tdk=0 很难解析求解。

    • 数值方法
      当解析解不可行时,则需要使用数值优化技术来求解一维问题 min ⁡ α > 0 φ ( α ) \min_{α > 0} \varphi (α) minα>0φ(α),即对 φ ( α ) \varphi (α) φ(α) 这个单变量函数进行优化。常用的方法有:

      • 基于导数的方法 (需要计算 φ ′ ( α ) \varphi '(α) φ(α) φ ′ ′ ( α ) \varphi ''(α) φ′′(α))
        1. 牛顿法:
          α j + 1 = α j − φ ′ ( α j ) / φ ′ ′ ( α j ) \alpha_{j+1}= \alpha_{j} - \varphi '(α_{j}) / \varphi ''(\alpha_{j}) αj+1=αjφ(αj)/φ′′(αj) 迭代该过程直至收敛3,令 α k = α j + 1 \alpha_k=\alpha_{j+1} αk=αj+1
        2. 割线法:
          α j + 1 = α j − φ ′ ( α j ) ( α j − α j − 1 ) φ ′ ( α j ) − φ ′ ( α j − 1 ) \alpha_{j+1} = \alpha_j - \frac{\varphi'(\alpha_j) (\alpha_j - \alpha_{j-1})}{\varphi'(α_j) - \varphi'(\alpha_{j-1})} αj+1=αjφ(αj)φ(αj1)φ(αj)(αjαj1)
        3. 求根方法:
          使用二分法或区间收缩法在某个区间内寻找 φ ′ ( α ) = 0 \varphi'(\alpha) = 0 φ(α)=0 的根。需要先确定一个包含根的区间 [ a , b ] [a, b] [a,b] 使得 φ ′ ( a ) \varphi'(a) φ(a) φ ′ ( b ) \varphi '(b) φ(b) 异号。
      • 不依赖导数的方法 (只需要 φ ( α ) \varphi(\alpha) φ(α) 的函数值)
        1. 黄金分割法:
          如果已知 φ ( α ) φ(α) φ(α) 在某个区间 [ a , b ] [a, b] [a,b] 内是单峰的,黄金分割法可以通过不断缩小这个区间来逼近最小值点。
        2. 抛物线插值法:
          通过三个点 ( α 1 , φ ( α 1 ) ) (\alpha_1, \varphi(\alpha_1)) (α1,φ(α1)), ( α 2 , φ ( α 2 ) ) (\alpha_2, \varphi(\alpha_2)) (α2,φ(α2)), ( α 3 , φ ( α 3 ) ) (\alpha_3, \varphi(\alpha_3)) (α3,φ(α3)) 构造一个二次抛物线,然后取抛物线的极小值点作为新的近似。可以与黄金分割法结合以提高鲁棒性和效率。
  • 特点: \textbf{特点:} 特点:

    • 理论上的最优性: 在给定的搜索方向上,精确线搜索确保了目标函数的最大下降。
    • 计算成本高昂: 算法的每一次迭代都执行一次精确线搜索,总的计算量可能非常巨大。
    • 对噪声敏感: 如果梯度计算本身带有噪声,那么追求基于噪声方向的“精确”迭代步长意义不大。

非精确线搜索

  • 思想: \textbf{思想:} 思想:
    不追求找到一维最优解,而是寻找一个能给出“足够”函数值下降且满足某些简单条件的步长 α k α_k αk。目标是在保证算法收敛性的前提下,尽快找到一个可接受的步长,以节省计算时间。
  • 常用的准则和方法: \textbf{常用的准则和方法:} 常用的准则和方法:
    • Armijo 条件

      • 思想: 确保步长 α k α_k αk 能够带来目标函数值的显著下降。实际的函数下降量至少是基于当前点线性近似所预测的下降量的一个固定比例,其主要作用是防止步长 α k \alpha_k αk过大,避免函数值上升。
      • 数学表达: f ( x k + α k d k ) ≤ f ( x k ) + c 1 α k ∇ f ( x k ) T d k f(x_k + α_k d_k) \le f(x_k) + c_1 α_k \nabla f(x_k)^T d_k f(xk+αkdk)f(xk)+c1αkf(xk)Tdk 其中 c 1 c_1 c1 是一个小的正数 ( 0 < c 1 < 1 ) (0 < c_1 < 1) (0<c1<1) ∇ f ( x k ) T d k \nabla f(x_k)^T d_k f(xk)Tdk 是函数 f ( x ) f(x) f(x) x k x_k xk处沿 d k d_k dk 的方向导数( ∇ f ( x k ) T d k < 0 \nabla f(x_k)^T d_k <0 f(xk)Tdk<0
      • 推导过程:
        如下图所示,设 ϕ ( α ) = f ( x k + α d k ) \phi(\alpha)=f(x_k+\alpha d_k) ϕ(α)=f(xk+αdk) d k d_k dk方向上函数值的变化曲线,根据泰勒展开, ϕ ( α ) \phi(\alpha) ϕ(α) α = 0 \alpha = 0 α=0 附近的一阶近似为:
        ϕ ( α ) = ϕ ( 0 ) + α ϕ ′ ( 0 ) = ∇ f ( x k ) + α ∇ f ( x k ) T d k \begin{align*} \phi(\alpha) = \phi(0) + \alpha \phi'(0) = \nabla f(x_k)+ \alpha \nabla f(x_k)^T d_k \end{align*} ϕ(α)=ϕ(0)+αϕ(0)=f(xk)+αf(xk)Tdk l ( α ) = ∇ f ( x k ) + α ∇ f ( x k ) T d k l(\alpha) =\nabla f(x_k)+ \alpha \nabla f(x_k)^T d_k l(α)=f(xk)+αf(xk)Tdk l ( α ) l(\alpha) l(α)表示从点 ( 0 , f ( x k ) ) (0, f(x_k)) (0,f(xk)) 出发,斜率为 ∇ f ( x k ) T d k \nabla f(x_k)^T d_k f(xk)Tdk 的一条直线,即图中 l c 1 = 1 ( α ) l_{c_1=1}(\alpha) lc1=1(α)表示的直线。为了保证在步长 α k \alpha_k αk处函数有足够的下降量,即满足 ϕ ( α k ) ≤ l ( α k ) \phi(\alpha_k) \le l(\alpha_k) ϕ(αk)l(αk),Armijo条件对直线 l ( α ) l(\alpha) l(α)的斜率利用比例 c 1 c_1 c1进行调节,即直线 l ( α ) l(\alpha) l(α)表示为:
        l ( α ) = ∇ f ( x k ) + c 1 α ∇ f ( x k ) T d k l(\alpha) = \nabla f(x_k)+ c_1 \alpha \nabla f(x_k)^T d_k l(α)=f(xk)+c1αf(xk)Tdk 则直线 l ( α ) l(\alpha) l(α)在图中 l c 1 = 1 ( α ) l_{c_1=1}(\alpha) lc1=1(α) l c 1 = 0 ( α ) l_{c_1=0}(\alpha) lc1=0(α)之间变动。根据不同的 c 1 c_1 c1值,在 ϕ ( α k ) ≤ l ( α k ) \phi(\alpha_k) \le l(\alpha_k) ϕ(αk)l(αk)条件下,可以获得不同的可接受步长区间,则 α k \alpha_k αk在该区间内取值。
        在这里插入图片描述
    • 曲率条件

      • 思想: Armijo条件允许 α \alpha α 任意小,当 α \alpha α 趋近于0,其总能满足Armijo条件,但会导致算法几乎没有进展。基于该缺陷,曲率条件则通过约束确保步长 α \alpha α 不能太小。它要求在新的迭代点 x k + α k d k x_k + \alpha_k d_k xk+αkdk 处的方向导数(沿 d k d_k dk 方向)比在 x k x_k xk 处的方向导数要平缓。
      • 数学表达: ∇ f ( x k + α k d k ) T d k ≥ c 2 ∇ f ( x k ) T d k \nabla f(x_k + \alpha_k d_k)^T d_k \ge c_2 \nabla f(x_k)^T d_k f(xk+αkdk)Tdkc2f(xk)Tdk 其中 c 2 c_2 c2 是一个常数,满足 c 1 < c 2 < 1 c_1 < c_2 < 1 c1<c2<1。在梯度下降法中 c 2 c_2 c2 常取 0.9 0.9 0.9,在拟牛顿法中常取 0.1 0.1 0.1
      • 推导过程:
        ϕ ( α ) = f ( x k + α d k ) \phi(\alpha)=f(x_k+\alpha d_k) ϕ(α)=f(xk+αdk) d k d_k dk方向上函数值的变化曲线,由链式法则可知 ϕ ( α ) \phi(\alpha) ϕ(α) 在新迭代点 ( α > 0 \alpha > 0 α>0) 处的斜率为:
        ϕ ′ ( α ) = ∇ f ( x k + α d k ) T d k \phi'(\alpha) = \nabla f(x_k+\alpha d_k)^T d_k ϕ(α)=f(xk+αdk)Tdk
        ϕ ′ ( 0 ) = ∇ f ( x k ) T d k \phi'(0) = \nabla f(x_k)^T d_k ϕ(0)=f(xk)Tdk 表示曲线在起点( α = 0 \alpha=0 α=0)处,沿 d k d_k dk 方向的斜率,其值为负。曲率条件设计了公式: ϕ ′ ( α ) ≥ c 2 ϕ ′ ( 0 ) \phi'(\alpha) \ge c_2 \phi'(0) ϕ(α)c2ϕ(0) 其意味着新迭代点的斜率 ϕ ′ ( α ) \phi'(\alpha) ϕ(α) 必须比当前点的斜率 ϕ ′ ( 0 ) \phi'(0) ϕ(0) “更不负一些”(即更平坦或值为正),该条件保证了 α \alpha α值足够大。将 ϕ ′ ( α k ) \phi'(\alpha_k) ϕ(αk) ϕ ′ ( 0 ) \phi'(0) ϕ(0)的定义代入上述条件,就得到了最终的曲率条件公式: ∇ f ( x k + α k d k ) T d k ≥ c 2 ∇ f ( x k ) T d k \nabla f(x_k + \alpha_k d_k)^T d_k \ge c_2 \nabla f(x_k)^T d_k f(xk+αkdk)Tdkc2f(xk)Tdk
        \qquad 下图展示了曲率条件所约束的移动步长 α \alpha α的取值区间。其通俗解释为:将函数想象成一个具有不同坡度的山谷(负坡度为下降方向)。在当前位置,我们向下走的坡度较大( − 30 -30 30度),走了一步之后,在新的位置其坡度仍然较大( − 28 -28 28度),为了快速到达谷底,那么我们仍然往下走。一直走到某个位置,其坡度较小( − 5 -5 5度),此时我们认为当前地势比较平坦,再往下走所下降的高度比较小(被认为比较接近谷底),因此,我们选择该位置为分界点,超过该分界点的区域被认为是我们想要到达的地方,即图中可接受步长区间。
        在这里插入图片描述
    • Wolfe 条件

      • 思想: 一般而言, W o l f e 条件 = A r m i j o 条件 + 曲率条件 Wolfe条件=Armijo条件+曲率条件 Wolfe条件=Armijo条件+曲率条件 为了避免 Armijo 条件可能导致步长 α k \alpha_k αk过小的问题,Wolfe 条件在Armijo条件的基础上引入了曲率条件,从而保证迭代步长既不过大,也不过小。
      • 数学表达:
        f ( x k + α k d k ) ≤ f ( x k ) + c 1 α k ∇ f ( x k ) T d k ( A r m i j o 条件 ) ∇ f ( x k + α k d k ) T d k ≥ c 2 ∇ f ( x k ) T d k ( 曲率条件 ) \begin{align*} & f(x_k + \alpha_k d_k) \le f(x_k) + c_1 \alpha_k \nabla f(x_k)^T d_k (Armijo条件) \\ & \nabla f(x_k + \alpha_k d_k)^T d_k ≥ c_2 \nabla f(x_k)^T d_k (曲率条件) \end{align*} f(xk+αkdk)f(xk)+c1αkf(xk)Tdk(Armijo条件)f(xk+αkdk)Tdkc2f(xk)Tdk(曲率条件) 其中 0 < c 1 < c 2 < 1 0<c_1< c_2<1 0<c1<c2<1。Wolfe条件所刻画的取值区间如下图所示。
        在这里插入图片描述
    • 强 Wolfe 条件

      • 思想: 标准曲率条件是单侧约束,其给斜率设定了一个“下限”,只约束了步长的最小值,当Armijo条件约束过大时, ϕ ′ ( α ) \phi'(\alpha) ϕ(α)可能取值较大。强Wolfe条件中的曲率条件则采用了双侧约束的强曲率条件,通过同时给新斜率 ϕ ′ ( α ) \phi'(\alpha) ϕ(α) 设定了下限和上限,使得新斜率的绝对值足够小,即 ϕ ′ ( α ) \phi'(\alpha) ϕ(α)被限制在一个靠近0的区间内。
      • 数学表达:
        f ( x k + α k d k ) ≤ f ( x k ) + c 1 α k ∇ f ( x k ) T d k ( A r m i j o 条件 ) ∣ ∇ f ( x k + α k d k ) T d k ∣ ≤ c 2 ∣ ∇ f ( x k ) T d k ∣ ( 强曲率条件 ) \begin{align*} & f(x_k + \alpha_k d_k) \le f(x_k) + c_1 \alpha_k \nabla f(x_k)^T d_k (Armijo条件) \\ & |\nabla f(x_k + \alpha_k d_k)^T d_k| \le c_2 |\nabla f(x_k)^T d_k| (强曲率条件) \\ \end{align*} f(xk+αkdk)f(xk)+c1αkf(xk)Tdk(Armijo条件)∣∇f(xk+αkdk)Tdkc2∣∇f(xk)Tdk(强曲率条件) 由于 ∇ f ( x k ) T d k < 0 \nabla f(x_k)^T d_k<0 f(xk)Tdk<0,则上面的强曲率条件等价于:
        c 2 ∇ f ( x k ) T d k ≤ ∇ f ( x k + α k d k ) T d k ≤ − c 2 ∇ f ( x k ) T d k c_2 \nabla f(x_k)^T d_k \le \nabla f(x_k + \alpha_k d_k)^T d_k \le -c_2 \nabla f(x_k)^T d_k c2f(xk)Tdkf(xk+αkdk)Tdkc2f(xk)Tdk
    • Goldstein 条件

      • 思想: 与Wolfe设计思想类似,通过构造一个双边不等式给步长 α \alpha α限定一个取值区间。
      • 数学表达:
        f ( x k ) + ( 1 − c ) α k ∇ f ( x k ) T d k ≤ f ( x k + α k d k ) ≤ f ( x k ) + c α k ∇ f ( x k ) T d k f(x_k) + (1-c) \alpha_k \nabla f(x_k)^T d_k \le f(x_k + \alpha_k d_k) \le f(x_k) + c \alpha_k \nabla f(x_k)^T d_k f(xk)+(1c)αkf(xk)Tdkf(xk+αkdk)f(xk)+cαkf(xk)Tdk 其中 0 < c < 0.5 0 < c < 0.5 0<c<0.5。令 ϕ ( α ) = f ( x k + α d k ) \phi(\alpha)=f(x_k+\alpha d_k) ϕ(α)=f(xk+αdk),则上式等价为:
        ϕ ( 0 ) + ( 1 − c ) α k ϕ ′ ( 0 ) ≤ ϕ ( α k ) ≤ ϕ ( 0 ) + c α k ϕ ′ ( 0 ) \phi(0)+(1-c) \alpha_k \phi'(0) \le \phi(\alpha_k) \le \phi(0)+c \alpha_k \phi'(0) ϕ(0)+(1c)αkϕ(0)ϕ(αk)ϕ(0)+cαkϕ(0) 该式的右边为Armijo条件,保证了函数值是下降的。如下图所示,从几何的角度该式可理解为,从当前迭代点 α = 0 \alpha=0 α=0处引出两条不同斜率的直线,则可接受的移动步长位于两条直线所形成区间内。在这里插入图片描述
        图中 L 1 L_1 L1为上界线,斜率为 c ϕ ′ ( 0 ) c \phi'(0) cϕ(0),要求 ( α , ϕ ( α ) ) (\alpha,\phi(\alpha)) (α,ϕ(α))处的点位于该线的下方;
        L 2 L_2 L2为下界线,斜率为 ( 1 − c ) ϕ ′ ( 0 ) (1-c) \phi'(0) (1c)ϕ(0),要求 ( α , ϕ ( α ) ) (\alpha,\phi(\alpha)) (α,ϕ(α))处的点位于该线的上方。如图中所示,GoldStein条件可能会将极值点排除在 α \alpha α取值区间外。

信赖域法

  • 思想: \textbf{思想:} 思想: 线搜索方法的核心思想是先确定搜索方向,再寻找搜索步长( x k + 1 = x k + α k d k x_{k+1}=x_k+\alpha_k d_k xk+1=xk+αkdk),而当Hessian矩阵 H k H_k Hk非正定时,此时搜索方向可能不是下降方向,因此线搜索方法需要额外设计复杂的修正策略来处理这种情况。为了避免该问题,信赖域方法则遵循先确定范围,再确定步长( x k + 1 = x k + p k x_{k+1}=x_k+p_k xk+1=xk+pk) 的步骤来寻找优化问题。其大致思路为:在当前迭代点 x k x_k xk 附近定义一个“信赖域”(通常是球形或椭球形),在该区域内用一个简单模型 m k m_k mk(如二次模型)近似目标函数。然后求解这个模型在信赖域内的最优解,得到迭代位移 p k p_k pk

  • 数学表达: \textbf{数学表达:} 数学表达:
    在当前迭代点 x k x_k xk处,构建一个二次模型 m k ( p ) m_k(p) mk(p)来近似 f ( x k + p ) f(x_k+p) f(xk+p)
    m k ( p ) = f ( x k ) + ∇ f ( x k ) T p + 1 2 p T B k p m_k(p)=f(x_k)+\nabla f(x_k)^T p + \frac{1}{2} p^T B_k p mk(p)=f(xk)+f(xk)Tp+21pTBkp 其中, B k B_k Bk 是一个近似Hessian矩阵的对称矩阵。则该优化问题将转换为一个如下形式的子优化问题:
    min ⁡ m k ( p ) s . t . ∣ ∣ p ∣ ∣ ≤ Δ k \begin{align*} \min \qquad &m_k(p) \\ s.t. \qquad & ||p||\le \Delta_k \end{align*} mins.t.mk(p)∣∣p∣∣Δk 其中, Δ k \Delta_k Δk为信赖域半径。该子优化问题的解 p k p_k pk即为本次迭代的候选解。

  • 算法框架: \textbf{算法框架:} 算法框架:
    1. 初始化迭代点 x 0 x_0 x0,初始信赖域半径 Δ 0 \Delta_0 Δ0,设置接受/拒绝阈值 η \eta η η ∈ ( 0 , 1 / 4 ) \eta \in(0,1/4) η(0,1/4),最大迭代次数,最大信赖域半径 Δ ^ \hat{\Delta} Δ^
    2. 构建二次模型 m k ( p ) = f ( x k ) + ∇ f ( x k ) T p + 1 2 p T B k p m_k(p) = f(x_k) + \nabla f(x_k)^T p + \frac{1}{2} p^T B_k p mk(p)=f(xk)+f(xk)Tp+21pTBkp;
    3. 近似求解 min ⁡ ∥ p ∥ ≤ Δ k m k ( p ) \min_{\|p\| \le \Delta_k} m_k(p) minpΔkmk(p),得到候选步长 p k p_k pk;
    4. 计算 ρ k \rho_k ρk ρ k = f ( x k ) − f ( x k + p k ) m k ( 0 ) − m k ( p k ) \rho_k = \frac{f(x_k) - f(x_k + p_k)}{m_k(0) - m_k(p_k)} ρk=mk(0)mk(pk)f(xk)f(xk+pk);
    5. 更新下一迭代点 x k + 1 x_{k+1} xk+1:
    \qquad 如果 ρ k > η \rho_k > \eta ρk>η:
    \qquad\qquad x k + 1 = x k + p k x_{k+1} = x_k + p_k xk+1=xk+pk (接受步长)
    \qquad 否则:
    \qquad\qquad x k + 1 = x k x_{k+1} = x_k xk+1=xk (拒绝步长,停在原地)
    6. 更新下一信赖域半径 Δ k + 1 \Delta_{k+1} Δk+1
    \qquad 如果 ρ k < 1 / 4 \rho_k < 1/4 ρk<1/4:
    \qquad\qquad Δ k + 1 = ( 1 / 4 ) Δ k \Delta_{k+1} = (1/4) \Delta_k Δk+1=(1/4)Δk (模型差,缩小半径)
    \qquad 如果 ρ k > 3 / 4 \rho_k > 3/4 ρk>3/4 并且 ∥ p k ∥ = Δ k \|p_k\| = \Delta_k pk=Δk (在边界上取得好效果):
    \qquad\qquad Δ k + 1 = min ⁡ ( 2 Δ k , Δ ^ ) \Delta_{k+1} = \min(2 \Delta_k, \hat{\Delta}) Δk+1=min(2Δk,Δ^) (模型好,扩大半径)
    \qquad 否则:
    \qquad\qquad Δ k + 1 = Δ k \Delta_{k+1} = \Delta_k Δk+1=Δk (模型一般,保持半径不变)
    7. 检查是否达到终止准则,是,停止迭代;否,转到步骤 2 2 2.

  • 特点: \textbf{特点:} 特点:

    • 具有较好的全局收敛性
    • 能处理Hessian矩阵非正定的情况

  1. 运用到的转置的性质:

    • ( M + N ) T = M T + N T (M+N)^T = M^T + N^T (M+N)T=MT+NT
    • ( M N ) T = N T M T (MN)^T = N^T M^T (MN)T=NTMT
    • ( M N P ) T = P T N T M T (MNP)^T = P^T N^T M^T (MNP)T=PTNTMT
    ↩︎
  2. 几个前提:

    • r k + 1 r_{k+1} rk+1推导: r k + 1 = b − A x k + 1 = b − A ( x k + α k p k ) = ( b − A x k ) − α k A p k = r k − α k A p k r_{k+1} = b - Ax_{k+1} = b - A(x_k + α_k p_k) = (b - Ax_k) - α_k A p_k = r_k - α_k A p_k rk+1=bAxk+1=bA(xk+αkpk)=(bAxk)αkApk=rkαkApk)
    • 残差的正交性: r i T r j = 0 , ∀ i ≠ j r_i^T r_j = 0, \forall i ≠ j riTrj=0,i=j,即迭代产生的残差序列是相互正交的。
    • 残差与先前搜索方向的关系: r k T p j = 0 , ∀ j < k r_k^T p_j = 0, \forall j < k rkTpj=0,j<k
    ↩︎
  3. 收敛性准则:

    • 一阶导数足够小: ∣ φ ′ ( α j + 1 ) ∣ < ε 1 , ε 1 = 10 − 5 |\varphi '(\alpha_{j+1})| < \varepsilon_1,\varepsilon_1=10^{-5} φ(αj+1)<ε1ε1=105
    • 步长变化足够小: ∣ α j + 1 − α j ∣ < ε 2 |\alpha_{j+1} - \alpha_j| < \varepsilon_2 αj+1αj<ε2
    • 函数值变化足够小 : ∣ φ ( α j + 1 ) − φ ( α j ) ∣ < ε 3 |\varphi(\alpha_{j+1}) - \varphi(\alpha_j)| < \varepsilon_3 φ(αj+1)φ(αj)<ε3
    • 达到最大迭代次数: j + 1 ≥ T m a x j + 1 \ge T_{max} j+1Tmax
    ↩︎
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值