【学习笔记】最优化方法

Armijo型非精确线性搜索求步长

Armijo 型线性搜索的数学条件定义为：

f(xk+αkdk)≤f(xk)+ραk∇f(xk)Tdkf(x_k + \alpha_k d_k) \leq f(x_k) + \rho \alpha_k \nabla f(x_k)^T d_k

参数说明：

  • ρ\rho 是预定义常数，且 ρ∈(0,1)\rho \in (0,1)
  • αk\alpha_k 表示第 kk 次迭代的步长
  • dkd_k 为搜索方向
  • ∇f(xk)\nabla f(x_k) 是目标函数在 xkx_k 处的梯度

公式作用：
该条件通过限制步长 αk\alpha_k 的取值，确保每次迭代中函数值 f(x)f(x) 充分下降，常用于梯度下降法等优化算法的步长选择。

算法2.3（Armijo型线性搜索）步骤：

步0：若 αk=1\alpha_k=1 满足上述数学条件，则取 αk=1\alpha_k=1；否则转步1
步1：给定常数 β>0\beta>0，ρ∈(0,1)\rho \in (0,1)，令 αk=β\alpha_k=\beta
步2：若 αk\alpha_k 满足条件（2.6），则得到步长 αk\alpha_k，终止计算；否则转步3
步3：令 αk=ραk\alpha_k=\rho\alpha_k，转步2

参数说明：

  • β\beta：初始试探步长（通常取 β=1\beta=1）
  • ρ\rho：步长缩减因子（控制步长收缩速度）
  • αk\alpha_k：第 kk 次迭代的候选步长

算法注释：

  1. 先尝试 αk=1\alpha_k=1，若满足下降条件（2.6），则直接采用单位步长
  2. 若不满足，则从较大步长 β\beta 开始，通过循环逐步缩小步长（αk←ραk\alpha_k \leftarrow \rho\alpha_k）
  3. 适用于梯度下降法、拟牛顿法等优化算法，保证迭代过程中目标函数值 f(xk)f(x_k) 充分下降
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
import numpy as np

def armijo_line_search(f, grad_f, x_k, d_k, beta=1, rho=0.5):
    """
    严格按照算法2.3实现的Armijo型线性搜索
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param x_k: 当前迭代点
    :param d_k: 当前搜索方向
    :param beta: 初始试探步长（默认值 1）
    :param rho: 步长缩减因子（默认值 0.5，范围 (0,1)）
    :return: 选定的步长 alpha_k
    """
    # 计算当前点的梯度
    grad_fk = grad_f(x_k)
    
    # 步0：先尝试 α_k = 1
    alpha_k = 1
    if f(x_k + alpha_k * d_k) <= f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        return alpha_k  # 满足 Armijo 条件，直接返回 1
    
    # 步1：否则使用 beta 作为初始步长
    alpha_k = beta
    
    # 步2 & 步3：逐步缩小步长，直到满足 Armijo 条件
    while f(x_k + alpha_k * d_k) > f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        alpha_k *= rho  # 步长缩小
    
    return alpha_k

# 示例目标函数和梯度
def f(x):
    return np.sum(x ** 2)  # 二次函数 f(x) = x^T x

def grad_f(x):
    return 2 * x  # 其梯度为 ∇f(x) = 2x

# 测试算法
x_k = np.array([1.0, 2.0])  # 初始点
d_k = -grad_f(x_k)  # 负梯度方向（梯度下降法）
alpha = armijo_line_search(f, grad_f, x_k, d_k)
print(f"选定步长: {alpha}")

Wolfe-Powell非精确线性搜索求步长

Wolfe-Powell 线性搜索的数学条件定义为：

  1. Armijo 条件（充分下降条件）：

f(xk+αkdk)≤f(xk)+ραk∇f(xk)Tdkf(x_k + \alpha_k d_k) \leq f(x_k) + \rho \alpha_k \nabla f(x_k)^T d_k

  1. 曲率条件：

∇f(xk+αkdk)Tdk≥σ∇f(xk)Tdk\nabla f(x_k + \alpha_k d_k)^T d_k \geq \sigma \nabla f(x_k)^T d_k

参数约束：

0<ρ<σ<10 < \rho < \sigma < 1

参数说明：

  • ρ\rho：控制下降量的常数（通常取 ρ=10−4\rho=10^{-4}）
  • σ\sigma：曲率条件的松弛系数（通常取 σ=0.1\sigma=0.1 或 0.90.9）
  • αk\alpha_k：第 kk 次迭代的候选步长
  • dkd_k：搜索方向（如梯度下降方向 dk=−∇f(xk)d_k = -\nabla f(x_k)）
  • ∇f(xk)\nabla f(x_k)：目标函数在 xkx_k 处的梯度

公式作用：

  1. Armijo 条件：确保步长 αk\alpha_k 使函数值 f(xk)f(x_k) 充分下降
  2. 曲率条件：保证步长 αk\alpha_k 不会过小（避免算法停滞），同时允许梯度模长适当增大
  3. 参数约束：ρ<σ\rho < \sigma 保证存在满足两条件的步长

补充说明：

  • 与 Armijo 准则相比，Wolfe-Powell 增加了曲率条件，避免步长选取过于保守
  • 广泛应用于共轭梯度法、拟牛顿法等需要保证“充分下降且合理步长”的优化算法
  • 条件中梯度内积 ∇f(xk)Tdk\nabla f(x_k)^T d_k 通常为负数（当 dkd_k 是下降方向时）

算法步骤（修正后规范表达）：

  1. 步0：若 αk=1\alpha_k=1 满足条件 (2.7)，则取 αk=1\alpha_k=1；否则转步1

    • 作用：优先尝试单位步长，若满足条件则直接采用

  2. 步1：给定参数 β>0\beta>0，ρ,ρ1∈(0,1)\rho,\rho_1 \in (0,1)，初始化 αk(0)\alpha_k^{(0)} 为集合 {βρj ∣ j∈Z}\{ \beta \rho^j \ | \ j \in \mathbb{Z} \} 中满足第一个不等式的最大候选步长，令 i=0i=0

    • 数学意义：在离散化网格 βρj\beta \rho^j 中搜索满足条件的最大步长

  3. 步2：若 αk(i)\alpha_k^{(i)} 满足第二个条件，则取 αk=αk(i)\alpha_k = \alpha_k^{(i)}；否则令 βk(i)=ρ−1αk(i)\beta_k^{(i)} = \rho^{-1} \alpha_k^{(i)}

    • 逻辑说明：若不满足条件，则放大搜索范围上限（βk(i)\beta_k^{(i)} 是原步长的 ρ−1\rho^{-1} 倍）

  4. 步3：在区间 [αk(i),βk(i)][\alpha_k^{(i)}, \beta_k^{(i)}] 内，通过集合 {αk(i)+ρ1j(βk(i)−αk(i)) ∣ j=0,1,2,...}\{ \alpha_k^{(i)} + \rho_1^j (\beta_k^{(i)} - \alpha_k^{(i)}) \ | \ j=0,1,2,... \} 重新搜索满足第一个不等式的最大步长 αk(i+1)\alpha_k^{(i+1)}，令 i=i+1i=i+1，转步2

    • 几何解释：在调整后的步长区间内，以 ρ1\rho_1 为细分因子进行二次搜索

参数说明：

  • β\beta：初始试探步长基准值（控制搜索范围）
  • ρ\rho：步长网格收缩因子（用于生成离散候选步长 βρj\beta \rho^j）
  • ρ1\rho_1：区间细分因子（控制二次搜索的精细程度）
  • αk(i)\alpha_k^{(i)}：第 ii 次迭代时的候选步长

解题思路与算法逻辑：

  1. 核心目标：通过动态调整步长区间，高效找到同时满足两个条件 (2.7) 的 αk\alpha_k
  2. 创新点：
    • 双重网格搜索：先用 ρ\rho 生成粗粒度候选步长，再用 ρ1\rho_1 在调整后的区间内细粒度搜索
    • 区间动态扩展：当候选步长不满足条件时，通过 βk(i)=ρ−1αk(i)\beta_k^{(i)} = \rho^{-1}\alpha_k^{(i)} 扩大搜索范围
  3. 收敛保证：ρ,ρ1<1\rho, \rho_1 < 1 确保步长在有限次迭代内缩小到满足条件

补充说明：

  1. 条件 (2.7) 应包含两个不等式（推测为类似 Wolfe-Powell 条件的 Armijo 条件和曲率条件）
  2. 该算法比标准 Armijo 搜索更复杂，适用于需要同时满足多个条件的优化问题
  3. 离散化集合 {βρj}\{ \beta \rho^j \} 和 {αk(i)+ρ1j(βk(i)−αk(i))}\{ \alpha_k^{(i)} + \rho_1^j (\beta_k^{(i)} - \alpha_k^{(i)}) \} 的设计平衡了计算效率与精度
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
import numpy as np

def wolfe_powell_search(f, grad_f, x_k, d_k, beta=1.0, rho=1e-4, sigma=0.1, rho1=0.5):
    """
    Wolfe-Powell 线性搜索算法
    :param f: 目标函数 f(x)
    :param grad_f: 目标函数的梯度函数 ∇f(x)
    :param x_k: 当前点 x_k
    :param d_k: 搜索方向 d_k
    :param beta: 初始试探步长
    :param rho: Armijo 条件控制参数（0 < rho < 1）
    :param sigma: 曲率条件控制参数（rho < sigma < 1）
    :param rho1: 区间细分因子（0 < rho1 < 1）
    :return: 满足 Wolfe-Powell 条件的步长 α_k
    """
    
    # 计算初始梯度内积 g_k_d = ∇f(x_k)^T d_k
    g_k = grad_f(x_k)  # 计算梯度
    g_k_d = np.dot(g_k, d_k)  # 梯度方向上的导数
    
    if g_k_d >= 0:
        raise ValueError("搜索方向 d_k 不是下降方向，应满足 ∇f(x_k)^T d_k < 0")
    
    # 先尝试 α_k = 1
    alpha_k = 1.0
    x_new = x_k + alpha_k * d_k  # 计算新点
    
    # 判断是否满足 Armijo 条件
    if f(x_new) <= f(x_k) + rho * alpha_k * g_k_d:
        # 判断是否满足曲率条件
        if np.dot(grad_f(x_new), d_k) >= sigma * g_k_d:
            return alpha_k  # 满足条件，直接返回
    
    # 步1：初始化 α_k^(0) 在集合 {βρ^j} 中满足 Armijo 条件的最大值
    alpha_k = beta
    while f(x_k + alpha_k * d_k) > f(x_k) + rho * alpha_k * g_k_d:
        alpha_k *= rho  # 缩小步长
    
    # 步2：如果 α_k^(0) 也满足曲率条件，则返回
    if np.dot(grad_f(x_k + alpha_k * d_k), d_k) >= sigma * g_k_d:
        return alpha_k
    
    # 步3：设定 β_k^(0) 并进行区间搜索
    beta_k = alpha_k / rho  # 设定上界
    i = 0  # 迭代计数
    while True:
        # 在 [α_k^(i), β_k^(i)] 之间搜索满足 Armijo 条件的最大步长
        alpha_new = alpha_k + rho1 ** i * (beta_k - alpha_k)
        if f(x_k + alpha_new * d_k) <= f(x_k) + rho * alpha_new * g_k_d:
            alpha_k = alpha_new
        else:
            beta_k = alpha_new
        
        # 检查是否满足曲率条件
        if np.dot(grad_f(x_k + alpha_k * d_k), d_k) >= sigma * g_k_d:
            return alpha_k
        
        i += 1  # 继续搜索

黄金分割法精确线性搜索求步长

输入：初始区间 [a1,b1][a_1, b_1]，精度要求 ε>0\varepsilon > 0
输出：极小值点所在区间 [ak,bk][a_k, b_k]

  1. 初始化：

    • 计算初始内点：

      λ1=a1+0.382(b1−a1)\lambda_1 = a_1 + 0.382(b_1 - a_1)

      μ1=a1+0.618(b1−a1)\mu_1 = a_1 + 0.618(b_1 - a_1)

    • 计算函数值 f(λ1)f(\lambda_1) 和 f(μ1)f(\mu_1)，令 k=1k=1

  2. 停止条件：

    • 若 bk−ak<εb_k - a_k < \varepsilon，终止计算，输出区间 [ak,bk][a_k, b_k]

  3. 区间更新：

    • 情况1：若 f(λk)>f(μk)f(\lambda_k) > f(\mu_k)

      • 更新左端点：ak+1=λka_{k+1} = \lambda_k，bk+1=bkb_{k+1} = b_k
      • 新内点：μk+1=ak+1+0.618(bk+1−ak+1)\mu_{k+1} = a_{k+1} + 0.618(b_{k+1} - a_{k+1})
      • 保留原右内点：λk+1=μk\lambda_{k+1} = \mu_k
      • 仅需计算 f(μk+1)f(\mu_{k+1})

    • 情况2：若 f(λk)≤f(μk)f(\lambda_k) \leq f(\mu_k)

      • 更新右端点：ak+1=aka_{k+1} = a_k，bk+1=μkb_{k+1} = \mu_k
      • 新内点：λk+1=ak+1+0.382(bk+1−ak+1)\lambda_{k+1} = a_{k+1} + 0.382(b_{k+1} - a_{k+1})
      • 保留原左内点：μk+1=λk\mu_{k+1} = \lambda_k
      • 仅需计算 f(λk+1)f(\lambda_{k+1})

  4. 迭代循环：

    • 令 k=k+1k = k+1，返回步骤2

参数说明：

  • λk,μk\lambda_k, \mu_k：第 kk 次迭代的黄金分割点（比例系数 0.382 和 0.618）
  • ε\varepsilon：区间长度精度阈值（控制算法终止条件）
  • [ak,bk][a_k, b_k]：第 kk 次迭代时的候选区间

算法原理：

  1. 黄金分割比：0.618 是 5−12\frac{\sqrt{5}-1}{2} 的近似值，确保每次迭代区间长度按固定比例 0.6180.618 收缩
  2. 单点复用：在区间更新时保留一个旧内点，只需计算一个新点函数值，减少计算量
  3. 单调收敛性：每次迭代排除不可能包含极小值的子区间，保证 bk−ak→0b_k - a_k \to 0

解题思路：

  1. 初始区间选择：需满足 f(x)f(x) 在 [a1,b1][a_1, b_1] 内为单峰函数

  2. 关键步骤验证：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
       # 伪代码示例：判断如何更新区间
       if f(lambda_k) > f(mu_k):
           a_new = lambda_k      # 舍弃左区间
           b_new = b_k
       else:
           a_new = a_k
           b_new = mu_k          # 舍弃右区间

    ```python
    import numpy as np

    def golden_section_search(f, a, b, epsilon=1e-5):
        """
        黄金分割法求解极小值所在区间
        :param f: 目标函数
        :param a: 初始区间左端点
        :param b: 初始区间右端点
        :param epsilon: 精度要求，控制停止条件
        :return: 极小值所在的区间 [a_k, b_k]
        """
        phi = (np.sqrt(5) - 1) / 2  # 黄金比例 0.618
        
        # 计算初始内点
        lambda_k = a + (1 - phi) * (b - a)
        mu_k = a + phi * (b - a)
        f_lambda = f(lambda_k)
        f_mu = f(mu_k)
        
        while (b - a) > epsilon:
            if f_lambda > f_mu:
                # 情况1：f(λ_k) > f(μ_k)，收缩左侧区间
                a = lambda_k
                lambda_k = mu_k  # 保留原右内点
                f_lambda = f_mu  # 复用已计算值
                mu_k = a + phi * (b - a)  # 计算新点
                f_mu = f(mu_k)
            else:
                # 情况2：f(λ_k) <= f(μ_k)，收缩右侧区间
                b = mu_k
                mu_k = lambda_k  # 保留原左内点
                f_mu = f_lambda  # 复用已计算值
                lambda_k = a + (1 - phi) * (b - a)  # 计算新点
                f_lambda = f(lambda_k)
        
        return a, b  # 返回最终区间

信赖域算法

修正后的规范算法步骤：

输入：初始点 x0x_0，最大信赖域半径 Δ‾\overline{\Delta}，精度阈值 ε≥0\varepsilon \geq 0，阈值参数 0<η1≤η2<10 < \eta_1 \leq \eta_2 < 1
输出：优化解 xkx_k

  1. 初始化：

    • 设定初始信赖域半径 Δ0∈(0,Δ‾)\Delta_0 \in (0, \overline{\Delta})，迭代计数器 k=0k=0

  2. 停止条件：

    • 计算梯度 gk=∇f(xk)g_k = \nabla f(x_k)
    • 若 ∥gk∥≤ε\|g_k\| \leq \varepsilon，终止算法，输出 xkx_k

  3. 子问题求解：

    • 在信赖域 ∥s∥≤Δk\|s\| \leq \Delta_k 内，近似求解：

      min⁡smk(s)=f(xk)+gkTs+12sTBks\min_{s} m_k(s) = f(x_k) + g_k^T s + \frac{1}{2}s^T B_k s

      得到试探步 sks_k，其中 BkB_k 是Hessian矩阵或其近似

  4. 接受步长：

    • 计算实际下降量 δf=f(xk)−f(xk+sk)\delta f = f(x_k) - f(x_k + s_k)
    • 计算预测下降量 δm=mk(0)−mk(sk)\delta m = m_k(0) - m_k(s_k)
    • 计算比值：

      rk=δfδmr_k = \frac{\delta f}{\delta m}

    • 更新迭代点：

      xk+1={xk+sk,if rk>0xk,otherwise x_{k+1} = \begin{cases} x_k + s_k, & \text{if } r_k > 0 \\ x_k, & \text{otherwise} \end{cases}

  5. 调整信赖域半径：

    Δk+1={min⁡(2Δk,Δ‾),if rk>η2Δk/4,if rk<η1Δk,otherwise \Delta_{k+1} = \begin{cases} \min(2\Delta_k, \overline{\Delta}), & \text{if } r_k > \eta_2 \\ \Delta_k/4, & \text{if } r_k < \eta_1 \\ \Delta_k, & \text{otherwise} \end{cases}

  6. 循环迭代：

    • 令 k=k+1k = k+1，返回步骤2

参数说明：

  • Δ‾\overline{\Delta}：信赖域半径的最大允许值（防止过度扩张）
  • η1,η2\eta_1, \eta_2：阈值参数，控制半径调整的敏感性（典型值：η1=0.1\eta_1=0.1, η2=0.75\eta_2=0.75）
  • BkB_k：二阶模型中的Hessian近似矩阵（如BFGS更新、SR1更新等）
  • rkr_k：实际下降与预测下降的比值（衡量模型可靠性）

算法原理：

  1. 信赖域框架：通过局部二次模型 mk(s)m_k(s) 近似目标函数，在半径 Δk\Delta_k 内寻找最优步长
  2. 半径动态调整：
    • rkr_k 接近 1：扩大信赖域（模型可信度高）
    • rkr_k 接近 0：缩小信赖域（模型不可靠）
  3. 子问题求解：常用Steihaug-CG方法或狗腿法(Dogleg)近似求解
  4. 全局收敛性：在适当条件下，算法满足 lim inf⁡k→∞∥gk∥=0\liminf_{k\to\infty} \|g_k\| = 0

解题思路：

  1. 关键公式验证：

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
       # 伪代码示例：计算接受率
       actual_reduction = f(xk) - f(xk + sk)
       predicted_reduction = - (g_k.T @ sk + 0.5 * sk.T @ B_k @ sk)
       rk = actual_reduction / predicted_reduction

    ```python
    import numpy as np

    def trust_region_method(f, grad_f, hess_f, x0, delta_max, epsilon=1e-5, eta1=0.1, eta2=0.75):
        """
        信赖域优化算法
        :param f: 目标函数
        :param grad_f: 目标函数梯度函数 ∇f(x)
        :param hess_f: 目标函数 Hessian 矩阵函数 B_k
        :param x0: 初始点 x_0
        :param delta_max: 最大信赖域半径 \overline{\Delta}
        :param epsilon: 终止条件，梯度范数阈值
        :param eta1: 阈值参数 η_1
        :param eta2: 阈值参数 η_2
        :return: 优化解 x_k
        """
        x_k = np.array(x0, dtype=float)  # 设定初始点，确保数据类型兼容计算
        delta_k = delta_max / 2  # 设定初始信赖域半径
        
        while True:
            g_k = grad_f(x_k)  # 计算梯度
            if np.linalg.norm(g_k) <= epsilon:
                break  # 终止条件
            
            B_k = hess_f(x_k)  # 计算 Hessian 矩阵
            
            # 处理 Hessian 矩阵奇异性，确保可逆
            try:
                s_k = -np.linalg.solve(B_k, g_k)  # 计算 Newton 步长
            except np.linalg.LinAlgError:
                s_k = -g_k  # 退化为梯度下降
            
            # 限制步长在信赖域半径内
            s_norm = np.linalg.norm(s_k)
            if s_norm > delta_k:
                s_k = (delta_k / s_norm) * s_k
            
            # 计算实际下降量
            delta_f = f(x_k) - f(x_k + s_k)
            
            # 计算预测下降量
            delta_m = g_k.T @ s_k + 0.5 * s_k.T @ B_k @ s_k
            
            # 计算 r_k
            r_k = delta_f / delta_m if delta_m != 0 else 0
            
            # 更新 x_k
            if r_k > 0:
                x_k = x_k + s_k
            
            # 更新信赖域半径
            if r_k > eta2:
                delta_k = min(2 * delta_k, delta_max)
            elif r_k < eta1:
                delta_k = delta_k / 4
            
        return x_k

最速下降法

输入：初始点 $ x_0 \in \mathbb{R}^n $，精度阈值 $ \epsilon > 0 $
输出：优化解 $ x_k $

  1. 初始化：

    • 设定 $ k = 0 $，给定初始点 $ x_0 $

  2. 停止条件：

    • 计算梯度 $ \nabla f(x_k) $
    • 若 $ |\nabla f(x_k)| \leq \epsilon $，终止算法，输出 $ x_k $

  3. 方向计算：

    • 取搜索方向为负梯度方向：

      dk=−∇f(xk)d_k = -\nabla f(x_k)

  4. 步长计算：

    • 通过一维精确线搜索求解：

      αk=arg⁡min⁡α>0f(xk+αdk)\alpha_k = \arg\min_{\alpha > 0} f(x_k + \alpha d_k)

  5. 迭代更新：

    • 更新迭代点：

      xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k

    • 令 $ k = k + 1 $，返回步骤2

关键点解析：

  1. 核心思想：

    • 沿目标函数下降最快的方向（负梯度方向）更新迭代点
    • 通过精确线搜索保证每步下降量最大

  2. 方向选择：

    • 方向 $ d_k = -\nabla f(x_k) $ 是当前点的最速下降方向
    • 梯度方向是函数值局部上升最快的方向，负梯度则为下降最快方向

  3. 步长计算：

    • 精确线搜索需解一维优化问题，对简单函数可解析求导（如二次函数）
    • 一般函数可能需要黄金分割法、抛物线插值法等数值方法

收敛性分析：

  1. 收敛条件：

    • 在目标函数 $ f(x) $ 为凸且光滑的条件下，算法全局收敛
    • 线性收敛速度：相邻迭代点误差满足 $ |x_{k+1} - x^| \leq \rho |x_k - x^| $，其中 $ 0 < \rho < 1 $

  2. 局限性：

    • 高维问题中可能出现“锯齿现象”，收敛速度显著下降
    • 对非凸函数可能收敛到局部极小或鞍点

解题思路（以二次函数为例）：

设目标函数为 $ f(x) = \frac{1}{2}x^T Q x + b^T x $，其中 $ Q $ 对称正定

  1. 梯度计算：

    ∇f(x)=Qx+b\nabla f(x) = Qx + b

  2. 方向更新：

    dk=−(Qxk+b)d_k = - (Qx_k + b)

  3. 精确步长解析解：

    αk=∇f(xk)T∇f(xk)∇f(xk)TQ∇f(xk)\alpha_k = \frac{\nabla f(x_k)^T \nabla f(x_k)}{\nabla f(x_k)^T Q \nabla f(x_k)}

示例（二维情况）：

  • 令 $ Q = \begin{bmatrix} 2 & 0 \ 0 & 1 \end{bmatrix} $, $ b = \begin{bmatrix} -1 \ 0 \end{bmatrix} $, $ x_0 = \begin{bmatrix} 0 \ 0 \end{bmatrix} $
  • 第1次迭代：
    • $ \nabla f(x_0) = \begin{bmatrix} -1 \ 0 \end{bmatrix} $
    • $ d_0 = \begin{bmatrix} 1 \ 0 \end{bmatrix} $
    • $ \alpha_0 = \frac{1}{2} $
    • $ x_1 = \begin{bmatrix} 0.5 \ 0 \end{bmatrix} $

对比与扩展：

  1. 与信赖域算法区别：

    • 最速下降法仅依赖梯度，而信赖域法需构建二阶模型
    • 最速下降法步长由线搜索确定，信赖域法步长受半径约束

  2. 改进方向：

    • 预处理技术：通过坐标变换减少条件数，缓解锯齿现象
    • 非精确线搜索：使用Armijo准则、Wolfe条件降低计算成本
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
import numpy as np

def steepest_descent(f, grad_f, x0, epsilon=1e-5, max_iter=1000):
    """
    最速下降法优化算法
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数 ∇f(x)
    :param x0: 初始点 x_0
    :param epsilon: 精度阈值
    :param max_iter: 最大迭代次数
    :return: 近似最优解 x_k
    """
    x_k = np.array(x0, dtype=float)  # 确保初始点为 NumPy 数组
    k = 0  # 迭代计数
    
    while k < max_iter:
        g_k = grad_f(x_k)  # 计算梯度
        if np.linalg.norm(g_k) <= epsilon:
            break  # 满足停止条件，返回结果
        
        # 计算步长 α_k 通过精确线搜索（以二次函数为例）
        d_k = -g_k  # 负梯度方向
        alpha_k = (g_k.T @ g_k) / (g_k.T @ hess_f(x_k) @ g_k)  # 解析步长计算
        
        # 更新 x_k
        x_k = x_k + alpha_k * d_k
        k += 1
    
    return x_k

基本Newton法

输入：初始点 $ x_0 \in \mathbb{R}^n $，精度阈值 $ \varepsilon > 0 $
输出：优化解 $ x_k $

  1. 初始化：

    • 设定 $ k = 0 $，给定初始点 $ x_0 $

  2. 停止条件：

    • 计算梯度 $ \nabla f(x_k) $
    • 若 $ |\nabla f(x_k)| \leq \varepsilon $，终止算法，输出 $ x_k $

  3. 方向计算：

    • 解线性方程组：

      ∇2f(xk)dk+∇f(xk)=0(公式3.3)\nabla^2 f(x_k) d_k + \nabla f(x_k) = 0 \quad \text{(公式3.3)}

    • 得到搜索方向 $ d_k $

  4. 迭代更新：

    • 更新迭代点：

      xk+1=xk+dkx_{k+1} = x_k + d_k

    • 令 $ k = k + 1 $，返回步骤2

关键点解析：

  1. 核心思想：

    • 利用目标函数的二阶导数（Hessian矩阵）构造搜索方向，实现局部快速收敛
    • 通过牛顿方程（公式3.3）直接求解优化方向，默认步长 $ \alpha_k = 1 $

  2. 方向计算：

    • 方向 $ d_k $ 是牛顿方向，满足 $ d_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k) $
    • 当Hessian矩阵正定时，方向为下降方向

  3. 步长特性：

    • 基本Newton法固定步长为1，适用于靠近极值点的情况
    • 改进方法（如阻尼Newton法）会引入线搜索调整步长

收敛性分析：

  1. 收敛速度：

    • 局部二次收敛：若初始点 $ x_0 $ 充分接近极小点且Hessian连续且正定，则

      ∥xk+1−x∗∥≤C∥xk−x∗∥2\|x_{k+1} - x^*\| \leq C \|x_k - x^*\|^2

    • 相比最速下降法（线性收敛），收敛速度显著提升

  2. 局限性：

    • Hessian矩阵可能非正定或奇异，导致方向非下降或无法求解
    • 初始点选择敏感，远离极值点时可能发散

解题示例（二次函数优化）：

设目标函数为 $ f(x) = \frac{1}{2}x^T Q x + b^T x $，其中 $ Q $ 对称正定

  1. 梯度与Hessian：

    ∇f(x)=Qx+b,∇2f(x)=Q\nabla f(x) = Qx + b, \quad \nabla^2 f(x) = Q

  2. 牛顿方向计算：

    Qdk+(Qxk+b)=0⇒dk=−xk−Q−1bQ d_k + (Qx_k + b) = 0 \Rightarrow d_k = -x_k - Q^{-1}b

  3. 迭代更新：

    xk+1=xk+dk=−Q−1b（一步收敛到极小点）x_{k+1} = x_k + d_k = -Q^{-1}b \quad \text{（一步收敛到极小点）}

数值算例：

  • 令 $ Q = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix} $, $ b = \begin{bmatrix} -3 \ -3 \end{bmatrix} $, $ x_0 = \begin{bmatrix} 0 \ 0 \end{bmatrix} $
  • 第1次迭代：
    • $ \nabla f(x_0) = \begin{bmatrix} -3 \ -3 \end{bmatrix} $
    • 解方程 $ Q d_0 = -\nabla f(x_0) \Rightarrow d_0 = \begin{bmatrix} 1 \ 1 \end{bmatrix} $
    • $ x_1 = x_0 + d_0 = \begin{bmatrix} 1 \ 1 \end{bmatrix} $（精确解）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
import numpy as np

def newton_method(f_grad, f_hessian, x0, epsilon):
    """
    使用牛顿法（Newton's Method）进行优化求解。
    :param f_grad: 目标函数的梯度函数，输入 x 返回梯度向量。
    :param f_hessian: 目标函数的 Hessian 矩阵函数，输入 x 返回 Hessian 矩阵。
    :param x0: 初始点，numpy 数组。
    :param epsilon: 终止精度阈值。
    :return: 近似最优解 x_k。
    """
    k = 0  # 迭代计数
    x_k = x0  # 初始化 x_k
    
    while True:
        grad = f_grad(x_k)  # 计算梯度 ∇f(x_k)
        if np.linalg.norm(grad) <= epsilon:  # 停止条件
            break
        
        hessian = f_hessian(x_k)  # 计算 Hessian 矩阵 ∇²f(x_k)
        
        # 求解线性方程组 ∇²f(x_k) d_k = -∇f(x_k)
        d_k = np.linalg.solve(hessian, -grad)
        
        # 更新迭代点
        x_k = x_k + d_k
        k += 1
    
    return x_k

阻尼Newton法

算法（阻尼Newton法）
输入：初始点 $ x_0 \in \mathbb{R}^n $，精度阈值 $ \varepsilon > 0 $
输出：优化解 $ x_k $

  1. 初始化：

    • 设定 $ k = 0 $，给定初始点 $ x_0 $

  2. 停止条件：

    • 计算梯度 $ \nabla f(x_k) $
    • 若 $ |\nabla f(x_k)| \leq \varepsilon $，终止算法，输出 $ x_k $

  3. 方向计算：

    • 解线性方程组：

      ∇2f(xk)dk+∇f(xk)=0\nabla^2 f(x_k) d_k + \nabla f(x_k) = 0

    • 得到搜索方向 $ d_k $

  4. 步长计算：

    • 通过线搜索（如Armijo准则、Wolfe条件）确定步长 $ \alpha_k $，确保目标函数值下降

  5. 迭代更新：

    • 更新迭代点：

      xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k

    • 令 $ k = k + 1 $，返回步骤2

关键点解析：

  1. 核心改进：

    • 相比基本Newton法（固定步长 $ \alpha_k = 1 $），阻尼法引入线搜索动态调整步长，提升全局收敛性。
    • 适用于初始点远离极值点或Hessian矩阵非正定的情况。

  2. 方向与步长的协同：

    • 方向：仍为牛顿方向 $ d_k = -[\nabla^2 f(x_k)]^{-1} \nabla f(x_k) $，需Hessian正定以保证下降性。
    • 步长：通过线搜索找到 $ \alpha_k \in (0, 1] $，使 $ f(x_k + \alpha_k d_k) < f(x_k) $。

  3. 适用场景：

    • 目标函数局部凸性较弱时，步长调整可避免迭代发散。
    • Hessian矩阵条件数较差时（如接近奇异），阻尼法更稳定。

对比基本Newton法：

特性 基本Newton法 阻尼Newton法
步长策略 固定步长 $ \alpha_k = 1 $ 动态步长 $ \alpha_k \in (0, 1] $
收敛性 局部二次收敛 全局收敛（需合理线搜索）
计算成本 低（无需步长搜索） 较高（需多次函数值计算）
初始点敏感性 高（依赖初始点接近极值点） 低（适应性更强）

收敛性分析：

  1. 全局收敛：
    • 在适当线搜索条件下（如Armijo准则），算法对任意初始点 $ x_0 $ 均收敛到局部极小点。
  2. 局部收敛速度：
    • 若初始点充分接近极小点且Hessian连续正定，仍保持局部二次收敛速度。

  • 解题示例（二次函数优化）：

    设目标函数 $ f(x) = \frac{1}{2}x^T Q x + b^T x $，其中 $ Q $ 对称正定

    1. 梯度与Hessian：

      ∇f(x)=Qx+b,∇2f(x)=Q\nabla f(x) = Qx + b, \quad \nabla^2 f(x) = Q

    2. 牛顿方向：

      Qdk+(Qxk+b)=0⇒dk=−xk−Q−1bQ d_k + (Qx_k + b) = 0 \Rightarrow d_k = -x_k - Q^{-1}b

    3. 步长计算：
      • 二次函数下，线搜索直接得 $ \alpha_k = 1 $（一步收敛）
    4. 迭代更新：

      xk+1=xk+dk=−Q−1b（精确极小点）x_{k+1} = x_k + d_k = -Q^{-1}b \quad \text{（精确极小点）}

    数值算例（同基本Newton法）：

    • $ Q = \begin{bmatrix} 2 & 1 \ 1 & 2 \end{bmatrix} $, $ b = \begin{bmatrix} -3 \ -3 \end{bmatrix} $, $ x_0 = \begin{bmatrix} 0 \ 0 \end{bmatrix} $
    • 第1次迭代：
      • $ \nabla f(x_0) = \begin{bmatrix} -3 \ -3 \end{bmatrix} $
      • 解方程 $ Q d_0 = -\nabla f(x_0) \Rightarrow d_0 = \begin{bmatrix} 1 \ 1 \end{bmatrix} $
      • 线搜索得 $ \alpha_0 = 1 $，故 $ x_1 = x_0 + d_0 = \begin{bmatrix} 1 \ 1 \end{bmatrix} $（精确解）
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import numpy as np

def newton_method(f, grad_f, hess_f, x0, epsilon=1e-6, max_iter=100):
    """
    牛顿法优化算法
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param hess_f: 目标函数的 Hessian 矩阵函数
    :param x0: 初始点 (n 维向量)
    :param epsilon: 收敛阈值
    :param max_iter: 最大迭代次数
    :return: 近似优化解 x_k
    """
    x_k = np.array(x0, dtype=np.float64)  # 确保是浮点类型
    k = 0  # 初始化迭代次数
    
    while k < max_iter:
        grad = grad_f(x_k)  # 计算梯度 ∇f(x_k)
        
        # 检查终止条件 ||∇f(x_k)|| ≤ ε
        if np.linalg.norm(grad) <= epsilon:
            break
        
        hess = hess_f(x_k)  # 计算 Hessian 矩阵 ∇²f(x_k)
        
        # 解线性方程组 ∇²f(x_k) * d_k = -∇f(x_k)
        try:
            d_k = np.linalg.solve(hess, -grad)
        except np.linalg.LinAlgError:
            print("Hessian 矩阵不可逆，终止优化。")
            break
        
        # 线搜索确定步长 α_k（使用 Armijo 准则）
        alpha = armijo_line_search(f, grad_f, x_k, d_k)
        
        # 迭代更新
        x_k = x_k + alpha * d_k
        k += 1
    
    return x_k

def armijo_line_search(f, grad_f, x_k, d_k, beta=1, rho=0.5):
    """
    严格按照算法2.3实现的Armijo型线性搜索
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param x_k: 当前迭代点
    :param d_k: 当前搜索方向
    :param beta: 初始试探步长（默认值 1）
    :param rho: 步长缩减因子（默认值 0.5，范围 (0,1)）
    :return: 选定的步长 alpha_k
    """
    # 计算当前点的梯度
    grad_fk = grad_f(x_k)
    
    # 步0：先尝试 α_k = 1
    alpha_k = 1
    if f(x_k + alpha_k * d_k) <= f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        return alpha_k  # 满足 Armijo 条件，直接返回 1
    
    # 步1：否则使用 beta 作为初始步长
    alpha_k = beta
    
    # 步2 & 步3：逐步缩小步长，直到满足 Armijo 条件
    while f(x_k + alpha_k * d_k) > f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        alpha_k *= rho  # 步长缩小
    
    return alpha_k

Newton-最速下降混合法

步骤1：给定初始点 x0∈Rnx_0 \in \mathbb{R}^n，精度 ϵ>0\epsilon > 0。令 k=0k = 0；

步骤2：若 ∥∇f(xk)∥≤ϵ\| \nabla f(x_k) \| \leq \epsilon，则得解 xkx_k，算法终止。否则，解线性方程组

∇2f(xk)d+∇f(xk)=0\nabla^2 f(x_k)d + \nabla f(x_k) = 0

若有解 dkd_k 且满足 ∇f(xk)⊤dk<0\nabla f(x_k)^\top d_k < 0，转步骤3；否则取 dk=−∇f(xk)d_k = -\nabla f(x_k)；

步骤3：由线性搜索计算步长 αk\alpha_k；

步骤4：令 xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k，k:=k+1k := k + 1，转步骤2。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import numpy as np

def newton_method(f, grad_f, hess_f, x0, epsilon=1e-6, max_iter=100):
    """
    牛顿法优化算法
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param hess_f: 目标函数的 Hessian 矩阵函数
    :param x0: 初始点 (n 维向量)
    :param epsilon: 收敛阈值
    :param max_iter: 最大迭代次数
    :return: 近似优化解 x_k
    """
    x_k = np.array(x0, dtype=np.float64)  # 确保是浮点类型
    k = 0  # 初始化迭代次数
    
    while k < max_iter:
        grad = grad_f(x_k)  # 计算梯度 ∇f(x_k)
        
        # 检查终止条件 ||∇f(x_k)|| ≤ ε
        if np.linalg.norm(grad) <= epsilon:
            break
        
        hess = hess_f(x_k)  # 计算 Hessian 矩阵 ∇²f(x_k)
        
        # 解线性方程组 ∇²f(x_k) * d_k = -∇f(x_k)
        try:
            d_k = np.linalg.solve(hess, -grad)
            if np.dot(grad, d_k) >= 0:
                d_k = -grad  # 若解不满足下降方向，则取 d_k = -∇f(x_k)
        except np.linalg.LinAlgError:
            d_k = -grad  # 若 Hessian 不可逆，使用梯度下降方向
        
        # 线搜索确定步长 α_k（使用 Armijo 准则）
        alpha = armijo_line_search(f, grad_f, x_k, d_k)
        
        # 迭代更新
        x_k = x_k + alpha * d_k
        k += 1
    
    return x_k

def armijo_line_search(f, grad_f, x_k, d_k, beta=1, rho=0.5):
    """
    严格按照算法2.3实现的Armijo型线性搜索
    :param f: 目标函数
    :param grad_f: 目标函数的梯度函数
    :param x_k: 当前迭代点
    :param d_k: 当前搜索方向
    :param beta: 初始试探步长（默认值 1）
    :param rho: 步长缩减因子（默认值 0.5，范围 (0,1)）
    :return: 选定的步长 alpha_k
    """
    # 计算当前点的梯度
    grad_fk = grad_f(x_k)
    
    # 步0：先尝试 α_k = 1
    alpha_k = 1
    if f(x_k + alpha_k * d_k) <= f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        return alpha_k  # 满足 Armijo 条件，直接返回 1
    
    # 步1：否则使用 beta 作为初始步长
    alpha_k = beta
    
    # 步2 & 步3：逐步缩小步长，直到满足 Armijo 条件
    while f(x_k + alpha_k * d_k) > f(x_k) + rho * alpha_k * np.dot(grad_fk, d_k):
        alpha_k *= rho  # 步长缩小
    
    return alpha_k

修正Newton-LM法

步骤1 给定初始点 $ x_0 \in \mathbb{R}^n $，精度 $ \varepsilon > 0 $。令 $ k = 0 $。

步骤2 若 $ |\nabla f(x_k)| \leq \varepsilon $，则得解 $ x_k $，算法终止；否则，解线性方程组

Akd+∇f(xk)=0(3.3)A_k d + \nabla f(x_k) = 0 \quad (3.3)

得解 $ d_k $，其中 $ A_k = \nabla^2 f(x_k) + v_k I ，， v_k > 0 $ 使得 $ A_k $ 正定。

步骤3 通过线性搜索计算步长 $ \alpha_k $。

步骤4 令 $ x_{k+1} = x_k + \alpha_k d_k ，， k := k + 1 $，转至步骤2。

参数说明

在修正Newton法中：

  1. 为确保 $ A_k $ 的正定性，需选择足够大的 $ v_k > 0 $；
  2. 为确保算法收敛性，需限制 $ v_k $ 的上界，即要求

    vk≤C∥∇f(xk)∥v_k \leq C \|\nabla f(x_k)\|

    其中 $ C $ 为常数，但其具体确定较为困难。

修正形式的收敛性

  • Newton法的两种修正形式在较弱条件下具有：
    • 超线性收敛性
    • 二次收敛性
  • 存在多种其他修正形式

鞍点情况的特殊处理

当 $ x $ 为鞍点时（即 $ \nabla f(x) = 0 $ 但 $ \nabla^2 f(x) $ 不定）：

  1. 所有修正方法均失效
  2. 可取负曲率方向 $ d $，需满足：

    dT∇2f(x)d<0d^T \nabla^2 f(x) d < 0

    沿此方向搜索时，目标函数值必然下降

方法特性对比

优点 缺点
收敛速度快 对初始点要求高
- 计算量大

算法改进方向

通过修改Newton法：

  1. 保留优势：利用快速收敛特性
  2. 克服缺陷：降低对初始点的敏感性，减少计算量
  3. 已衍生出多种高效新算法

拟Newton法

步骤1 给定初始点 $ x_0 \in \mathbb{R}^n $，初始化矩阵 $ B_0 $ 或 $ H_0 $，精度阈值 $ \varepsilon > 0 $。令迭代计数器 $ k = 0 $。

步骤2 若满足终止准则 $ |\nabla f(x_k)| \leq \varepsilon $，输出解 $ x_k $ 并终止算法。

步骤3 求解线性方程组：

Bkd+∇f(xk)=0B_k d + \nabla f(x_k) = 0

得搜索方向 $ d_k $，或通过逆矩阵公式直接计算：

dk=−Hkgk(其中 gk=∇f(xk))d_k = -H_k g_k \quad (\text{其中 } g_k = \nabla f(x_k))

步骤4 通过线搜索计算步长 $ \alpha_k $，更新迭代点：

xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k

步骤5 按拟Newton条件（DFP/SRP1等）修正矩阵：

Bk+1=Bk+ΔB或Hk+1=Hk+ΔHB_{k+1} = B_k + \Delta B \quad \text{或} \quad H_{k+1} = H_k + \Delta H

使其满足拟Newton方程 Bk+1skB_{k+1}s_k 或 Hk+1yk=skH_{k+1}y_k=s_k 。令 $ k := k + 1 $，返回步骤2。

修正矩阵

对称秩1

SR修正公式形式1

Bk+1=Bk+(yk−Bksk)(yk−Bksk)T(yk−Bksk)TskB_{k+1}=B_k+\dfrac {(y_k-B_ks_k)(y_k-B_ks_k)^T} {(y_k-B_ks_k)^Ts_k}

SR修正公式形式2

Hk+1=Hk+(sk−Hkyk)(sk−Hkyk)T(sk−Hkyk)TykH_{k+1}=H_k+ \dfrac {(s_k-H_ky_k)(s_k-H_ky_k)^T} {(s_k-H_ky_k)^Ty_k}

对称秩2

BFGSBFGS 修正公式(迄今最好的拟牛顿修正公式)

Bk+1=Bk−BkskskTBkskTBksk+ykykTykTskB_{k+1}=B_k-\dfrac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\dfrac {y_ky_k^T}{y_k^Ts_k}

BFGSBFGS 的逆修正公式

Hk+1=(I−skykTykTsk)Hk(I−skykTykTsk)T+skskTykTskH_{k+1}=(I-\dfrac {s_ky_k^T}{y_k^Ts_k})H_k(I-\dfrac {s_ky_k^T}{y_k^Ts_k})^T+\dfrac {s_ks_k^T}{y_k^Ts_k}

DFPDFP 修正公式

Hk+1=Hk−HkykykTHkykTHkyk+skskTskTykH_{k+1}=H_k-\dfrac {H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\dfrac {s_ks_k^T}{s_k^Ty_k}

DFPDFP 的逆修正公式

Bk+1=(I−ykskTskTyk)Bk(I−ykskTskTyk)+ykykTskTykB_{k+1}=(I-\dfrac {y_ks_k^T}{s_k^Ty_k})B_k(I-\dfrac {y_ks_k^T}{s_k^Ty_k})+\dfrac {y_ky_k^T}{s_k^Ty_k}

BFGS算法

步骤1：给定初始点 x0∈Rnx_0 \in \mathbb{R}^n，初始对称正定矩阵 B0B_0，精度 ϵ>0\epsilon > 0。令 k=0k = 0。

步骤2：若 ∥∇f(xk)∥≤ϵ\| \nabla f(x_k) \| \leq \epsilon，则得解 xkx_k，算法终止。否则转下一步。

步骤3：解线性方程组

Bkd+∇f(xk)=0B_k d + \nabla f(x_k) = 0 \qquad

得解 dkd_k。

步骤4：由线性搜索计算步长 αk\alpha_k。

步骤5：令 xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k，若 ∥∇f(xk+1)∥≤ϵ\| \nabla f(x_{k+1}) \| \leq \epsilon，则得解 xk+1x_{k+1}，算法终止。否则由公式 Bk+1=Bk−BkskskTBkskTBksk+ykykTykTskB_{k+1}=B_k-\dfrac{B_ks_ks_k^TB_k}{s_k^TB_ks_k}+\dfrac {y_ky_k^T}{y_k^Ts_k} 计算 Bk+1B_{k+1}。

步骤6：令 k:=k+1k := k + 1，转步骤3。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
import numpy as np

def function(x):
    """ 目标函数 f(x) """
    return np.sum(x**2)  # 这里使用 f(x) = x^T x 作为示例

def gradient(x):
    """ 计算目标函数的梯度 \nabla f(x) """
    return 2 * x  # 目标函数 f(x) = x^T x 的梯度为 2x

#这里也可以用Wolfe-Powell搜索
def line_search(f, grad, xk, dk, alpha=1, rho=0.8, c=1e-4):
    """
    线性搜索计算步长 \alpha_k
    使用 Armijo 条件进行搜索，确保下降足够快。
    """
    while f(xk + alpha * dk) > f(xk) + c * alpha * np.dot(grad(xk).T, dk):
        alpha *= rho  # 按比例缩小步长
    return alpha

def bfgs(x0, B0, epsilon=1e-6, max_iter=100):
    """ BFGS 算法实现 """
    xk = x0  # 初始点
    Bk = B0  # 初始对称正定矩阵
    k = 0
    
    while np.linalg.norm(gradient(xk)) > epsilon and k < max_iter:
        # 计算搜索方向 d_k
        dk = -np.linalg.solve(Bk, gradient(xk))
        
        # 线性搜索确定步长 alpha_k
        alpha_k = line_search(function, gradient, xk, dk)
        
        # 更新 x
        xk_new = xk + alpha_k * dk
        
        # 计算梯度变化 y_k 和步长向量 s_k
        sk = xk_new - xk
        yk = gradient(xk_new) - gradient(xk)
        
        # 终止条件检查
        if np.linalg.norm(gradient(xk_new)) <= epsilon:
            return xk_new
        
        # 计算 BFGS 公式中的矩阵更新项
        Bk_sk = Bk @ sk
        sk_Bk_sk = sk.T @ Bk_sk
        yk_ykT = np.outer(yk, yk)
        yk_sk = yk.T @ sk
        
        # BFGS 更新公式
        Bk = Bk - np.outer(Bk_sk, Bk_sk) / sk_Bk_sk + yk_ykT / yk_sk
        
        # 进入下一轮迭代
        xk = xk_new
        k += 1
    
    return xk

# 设置初始值
n = 2  # 变量维数
x0 = np.array([1.0, 1.5])  # 初始点
B0 = np.eye(n)  # 初始正定矩阵，单位矩阵

# 运行 BFGS 算法
optimal_x = bfgs(x0, B0)
print("最优解:", optimal_x)

DFP算法

步骤1：给定初始点 x0∈Rnx_0 \in \mathbb{R}^n，初始对称正定矩阵 H0H_0，精度 ϵ>0\epsilon > 0。令 k=0k = 0。

步骤2：若 ∥∇f(xk)∥≤ϵ\| \nabla f(x_k) \| \leq \epsilon，则输出解 xkx_k，算法终止；否则继续。

步骤3：计算搜索方向：

dk=−Hk∇f(xk)d_k = -H_k \nabla f(x_k)

（或通过解线性方程组 Bkd+∇f(xk)=0B_k d + \nabla f(x_k) = 0，其中 BkB_k 由公式 Bk+1=(I−ykskTskTyk)Bk(I−ykskTskTyk)+ykykTskTykB_{k+1}=(I-\dfrac {y_ks_k^T}{s_k^Ty_k})B_k(I-\dfrac {y_ks_k^T}{s_k^Ty_k})+\dfrac {y_ky_k^T}{s_k^Ty_k} 计算）

步骤4：通过线性搜索确定步长 αk\alpha_k。

步骤5：更新迭代点：

xk+1=xk+αkdkx_{k+1} = x_k + \alpha_k d_k

若 ∥∇f(xk+1)∥≤ϵ\| \nabla f(x_{k+1}) \| \leq \epsilon，输出解 xk+1x_{k+1} 并终止；否则按公式 Hk+1=Hk−HkykykTHkykTHkyk+skskTskTykH_{k+1}=H_k-\dfrac {H_ky_ky_k^TH_k}{y_k^TH_ky_k}+\dfrac {s_ks_k^T}{s_k^Ty_k} 计算 Hk+1H_{k+1}。

步骤6：令 k:=k+1k := k + 1，转步骤3。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
import numpy as np

def line_search(f, x, d, grad_x, alpha_init=1.0, rho=0.8, c=1e-4, max_backtrack=20):
    """
    回溯线搜索算法寻找满足Armijo条件的步长
    
    参数:
    f -- 目标函数
    x -- 当前点
    d -- 搜索方向
    grad_x -- 当前点的梯度
    alpha_init -- 初始步长
    rho -- 步长衰减因子
    c -- Armijo条件常数
    max_backtrack -- 最大回溯次数
    
    返回:
    alpha -- 选择的步长
    """
    alpha = alpha_init
    f_x = f(x)
    grad_dot_d = np.dot(grad_x, d)
    
    # 确保方向是下降方向
    if grad_dot_d >= 0:
        return 0.0
    
    for _ in range(max_backtrack):
        x_new = x + alpha * d
        f_new = f(x_new)
        if f_new <= f_x + c * alpha * grad_dot_d:
            return alpha
        alpha *= rho
    return alpha  # 返回最小找到的alpha

def dfp_algorithm(f, grad_f, x0, H0, epsilon, max_iter=1000, line_search_params=None):
    """
    DFP拟牛顿法实现
    
    参数:
    f -- 目标函数
    grad_f -- 梯度函数
    x0 -- 初始点 (numpy数组)
    H0 -- 初始正定矩阵 (numpy数组)
    epsilon -- 收敛精度
    max_iter -- 最大迭代次数
    line_search_params -- 线搜索参数字典
    
    返回:
    x_k -- 近似最优解
    """
    # 步骤1: 初始化参数
    x_k = x0.copy()
    H_k = H0.copy()
    grad_k = grad_f(x_k)
    n = x_k.shape[0]
    
    # 线搜索参数设置
    if line_search_params is None:
        line_search_params = {'alpha_init': 1.0, 'rho': 0.5, 'c': 1e-4}
    alpha_init = line_search_params['alpha_init']
    rho = line_search_params['rho']
    c = line_search_params['c']
    
    for k in range(max_iter):
        # 步骤2: 检查梯度收敛条件
        if np.linalg.norm(grad_k) <= epsilon:
            print(f"在迭代 {k} 次后收敛")
            return x_k
        
        # 步骤3: 计算搜索方向
        d_k = -H_k @ grad_k  # 计算搜索方向
        
        # 步骤4: 执行线搜索
        alpha_k = line_search(f, x_k, d_k, grad_k, alpha_init, rho, c)
        if alpha_k <= 1e-12:
            print("步长过小，终止计算")
            return x_k
        
        # 步骤5: 更新迭代点
        x_next = x_k + alpha_k * d_k
        grad_next = grad_f(x_next)
        
        # 检查新点的收敛条件
        if np.linalg.norm(grad_next) <= epsilon:
            print(f"在迭代 {k+1} 次后收敛")
            return x_next
        
        # 计算向量s_k和y_k
        s_k = x_next - x_k
        y_k = grad_next - grad_k
        sy_dot = np.dot(s_k, y_k)  # 计算s_k^T y_k
        
        # 步骤5续: 更新H矩阵
        if abs(sy_dot) < 1e-8:  # 防止数值不稳定
            H_next = H_k.copy()
        else:
            Hy = H_k @ y_k
            yHy = y_k @ Hy
            if abs(yHy) < 1e-8:
                term1 = 0.0
            else:
                term1 = np.outer(Hy, Hy) / yHy  # 第一修正项
            
            term2 = np.outer(s_k, s_k) / sy_dot  # 第二修正项
            H_next = H_k - term1 + term2  # 更新H矩阵
        
        # 步骤6: 迭代变量更新
        x_k = x_next
        grad_k = grad_next
        H_k = H_next.copy()
    
    print(f"达到最大迭代次数 {max_iter} 次未收敛")
    return x_k

# 示例使用（二次函数优化）
def quadratic(x):
    return x[0]**2 + 2*x[1]**2 + x[0]*x[1]

def quadratic_grad(x):
    return np.array([2*x[0] + x[1], 4*x[1] + x[0]])

# 初始参数设置
x0 = np.array([2.0, 1.0])
H0 = np.eye(2)
epsilon = 1e-6

# 执行DFP算法
solution = dfp_algorithm(quadratic, quadratic_grad, x0, H0, epsilon)
print("优化解:", solution)
print("最终梯度:", quadratic_grad(solution))

Broyden族

对互为对偶的BFGS公式和DFP公式进行加权组合得到一类重要的拟Newton矩阵修正公式：Broden族

Bk+1θk=θkBk+1BFGS+(1−θk)Bk+1DFPB_{k+1}^{\theta_k}=\theta_kB_{k+1}^{BFGS}+(1-\theta_k)B_{k+1}^{DFP}

Hk+1ϕk=ϕkHk+1BFGS+(1−ϕk)Hk+1DFPH_{k+1}^{\phi_k}=\phi_kH_{k+1}^{BFGS}+(1-\phi_k)H_{k+1}^{DFP}

其中 θk=1−1−ϕk1−ϕk+ϕkμk\theta_k=1-\dfrac {1-\phi_k}{1-\phi_k+\phi_k \mu_k} ， μk=ykTHkykskTBkskskTyk\mu_k=\dfrac{y_k^TH_ky_ks_k^TB_ks_k}{s_k^Ty_k}

由人工神经网络方法解微分方程导出的最优化问题

思维学一般把人类的大脑思维分为抽象（逻辑）思维，形象（直观）思维和灵感（顿悟）思维三种基本方式，人工神经网络就是模拟人类思维的形象思维方式去解决问题的。神经网络的基础是神经元，神经元是以生物神经系统的神经细胞为基础的生物模型。人们对生物神经系统进行研究以探讨人工智能的机制时，把神经元数学化，从而产生了神经元数学模型，大量相同形式的神经元连接在一起就组成了神经网络

神经网络是一个高度非线性动力学系统，虽然每个神经元的结构和功能都不复杂，但是神经网络的动态行为却是十分复杂的，因此，用神经网络可以表达实际物理世界的各种现象。

我们可以用有向图表示神经网络，这里有向图是节点与节点之间的有向连接。一般来说，神经网络有两种基本结构：前馈神经网络和递归神经网络。这里我们仅考虑前馈神经网络，它包括：

输入层：源节点构成输入层，它向网络提供输入信号

隐藏层：前馈神经网络有一层或多层隐藏节点，相应的节点称为隐藏神经元

输出层：该层给出相对于源节点的激活模式的网络输出

一个 m−n1−n2−qm-n_1-n_2-q 前馈神经网络，表示该网络有 mm 个源节点， 第一个隐藏层有 n1n_1 个神经元，第二个隐藏层有 n2n2 个神经元，输出层有 qq 个神经元

求解微分方程的方法是多种多样的，如下有一种用神经网络的方法求解微分方程的方法

考虑一阶常微分方程

dy(x)dx=f(x,y),x∈[0,1]\dfrac {dy(x)}{dx}=f(x,y),x\in [0,1]

其初始条件为

y(0)=y0y(0)=y_0

利用神经网络的方法我们可以得到一个近似解

yt(x,p)y_t(x,p)

可以用它去近似方程的精准解 y(x)y(x) ，其中实验解中的参数 pp 为神经网络方法所产生的参数，下面我们讨论神经网络实验解的表示

考虑一个 1−n−11-n-1 的前馈神经网络，一阶微分方程的神经网络实验解可以表示为

yt(x,p)=y0+xN(x,p)y_t(x,p)=y_0+xN(x,p)

其中 N(x,p)N(x,p) 为前馈神经网络的输出， xx 是单一输入， pp 是神经网络参数向量。前馈神经网络的输出为

N(x,p)=∑j=1nvjφ(zj),zj=wjx−θjN(x,p)= \sum_{j=1}^{n} v_jφ(z_j),z_j=w_jx-\theta_j

其中 wjw_j 是从源节点到隐藏节点 jj 的权重 ， vjv_j 是从隐藏节点 jj 到输出节点的权重， θj\theta_j 是隐藏节点 jj 的阈值， p=(w,v,θ)Tp=(w,v,\theta)^T 是一个 3n3n 维向量 ， φ(z)φ(z) 是 Sigmoid型激活函数

φ(z)=11+e−zφ(z)=\dfrac {1}{1+e^{-z}}

为了使得实验解能尽可能好地近似方程的解，我们要用最优化的方法去确定实验解中的参数 pp 。对于一阶微分方程，我们要将 xx 在 [0,1][0,1] 离散化，得到 mm 个值 $0 \leq x_1<x_2……\leq x_m=1 $ 以及 dyt(x,p)dx∣x=xi\left. \dfrac {dy_t(x,p)}{dx}\right|_{x=x_i} ，f(xi,yt(xi,p))f(x_i,y_t(x_i,p)) .其中 x1x_1 到 xmx_m 称为训练集元素.然后求解最优化问题

min⁡∑i=1m{dyt(x,p)dx∣x=xi−f(xi,yt(xi,p))}2\min\sum_{i=1}^m \{\left.\dfrac{dy_t(x,p)}{dx}\right|_{x=x_i}-f(x_i,y_t(x_i,p))\}^2

来确定 pp 的值

对正定二次型的共轭梯度法

步0（初始化）:
给定线性无关向量组 $ \mathbf{g} \in \mathbb{R}^n $，令

d0=−g0,k:=0\mathbf{d}_0 = -\mathbf{g}_0, \quad k := 0

步1（方向更新）:
计算搜索方向：

dk=−gk+∑j=0k−1(djTGgkdjTGdj)dj\mathbf{d}_k = -\mathbf{g}_k + \sum_{j=0}^{k-1} \left( \frac{\mathbf{d}_j^T \mathbf{G} \mathbf{g}_k}{\mathbf{d}_j^T \mathbf{G} \mathbf{d}_j} \right) \mathbf{d}_j

步2（终止条件）:
若 $ k = n-1 $，则停止；否则继续。

步3（线搜索）:
计算步长 $ \alpha_k $，使得：

αk=arg⁡min⁡αf(xk+αdk)\alpha_k = \arg\min_{\alpha} f(\mathbf{x}_k + \alpha \mathbf{d}_k)

步4（迭代更新）:

  1. 计算梯度 $ \nabla f(\mathbf{x}_k + \alpha_k \mathbf{d}_k) $

  2. 令 $ k := k + 1 $，返回步1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import numpy as np

def conjugate_direction_method(f, grad_f, G, x0):
    """
    共轭方向法求解无约束优化问题。
    
    参数：
    f : function - 目标函数 f(x)
    grad_f : function - 目标函数的梯度 ∇f(x)
    G : ndarray -  正定矩阵，表示二次型 f(x) 的 Hessian 矩阵
    x0 : ndarray - 初始点
    
    返回：
    x_opt : ndarray - 最优解
    """
    n = len(x0)  # 变量维度
    g = grad_f(x0)  # 计算初始梯度
    d = [-g]  # 初始化搜索方向 d_0 = -g_0
    x = x0  # 当前点
    
    for k in range(n):
        if k > 0:
            # 计算共轭方向
            d_k = -g
            for j in range(k):
                beta = (d[j].T @ G @ g) / (d[j].T @ G @ d[j])
                d_k += beta * d[j]
            d.append(d_k)
        
        # 线搜索，找到最优步长 α_k
        alpha_k = - (g.T @ d[k]) / (d[k].T @ G @ d[k])
        
        # 迭代更新
        x = x + alpha_k * d[k]
        g = grad_f(x)  # 计算新的梯度
        
        if k == n - 1:
            break  # 终止条件
    
    return x

# 示例使用
if __name__ == "__main__":
    # 目标函数 f(x) = 0.5 * x^T G x - b^T x
    G = np.array([[4, 1], [1, 3]])  # 正定矩阵
    b = np.array([1, 2])
    
    def f(x):
        return 0.5 * x.T @ G @ x - b.T @ x
    
    def grad_f(x):
        return G @ x - b
    
    x0 = np.array([2, 1])  # 初始点
    x_opt = conjugate_direction_method(f, grad_f, G, x0)
    print("最优解:", x_opt)

（FR）对正定二次型的共轭梯度法

算法步骤：

步0：给定线性无关向量组：$ g \in \mathbb{R}^n $，令 $ d_0 = -g_0 ，， k = 0 $。
步1：计算

dk=−gk+dk−1TGgkdk−1TGdk−1dk−1.d_k = -g_k + \frac{d_{k-1}^T G g_k}{d_{k-1}^T G d_{k-1}} d_{k-1}.

步2：若 $ k = n - 1 $，则停止。
步3：计算线搜索 $ f(x_k + \alpha d_k) $ 的最小值点 $ \alpha_k $。
步4：求梯度 $ \nabla f(x_k + \alpha_k d_k) $，令 $ k = k + 1 $，转步1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
import numpy as np

def conjugate_gradient_method(f, grad_f, G, x0):
    """
    共轭梯度法求解无约束优化问题。
    
    参数：
    f : function - 目标函数 f(x)
    grad_f : function - 目标函数的梯度 ∇f(x)
    G : ndarray -  正定矩阵，表示二次型 f(x) 的 Hessian 矩阵
    x0 : ndarray - 初始点
    
    返回：
    x_opt : ndarray - 最优解
    """
    n = len(x0)  # 变量维度
    g = grad_f(x0)  # 计算初始梯度
    d = -g  # 初始化搜索方向 d_0 = -g_0
    x = x0  # 当前点
    
    for k in range(n):
        # 线搜索，找到最优步长 α_k
        alpha_k = - (g.T @ d) / (d.T @ G @ d)
        
        # 迭代更新
        x = x + alpha_k * d
        g_new = grad_f(x)  # 计算新的梯度
        
        if k == n - 1:
            break  # 终止条件
        
        # 计算新的搜索方向
        beta_k = (d.T @ G @ g_new) / (d.T @ G @ d)
        d = -g_new + beta_k * d
        g = g_new  # 更新梯度
    
    return x

# 示例使用
if __name__ == "__main__":
    # 目标函数 f(x) = 0.5 * x^T G x - b^T x
    G = np.array([[4, 1], [1, 3]])  # 正定矩阵
    b = np.array([1, 2])
    
    def f(x):
        return 0.5 * x.T @ G @ x - b.T @ x
    
    def grad_f(x):
        return G @ x - b
    
    x0 = np.array([2, 1])  # 初始点
    x_opt = conjugate_gradient_method(f, grad_f, G, x0)
    print("最优解:", x_opt)

HS共轭梯度法

步骤1：
给定初始点 $ x_0 \in \mathbb{R}^n $，精度 $ \varepsilon > 0 $。
令 $ d_0 = -\nabla f(x_0) $，令 $ k := 0 $。

步骤2：
若 $ |\nabla f(x_k)| \leq \varepsilon $，则得解 $ x_k $，算法终止。否则转步骤3。

步骤3：
由线性搜索计算步长 $ \alpha_k $。

步骤4：
令 $ x_{k+1} = x_k + \alpha_k d_k $。

步骤5：
由式 $$
d_k = \begin{cases}
-\nabla f(x_0), & k=0 \
-\nabla f(x_k) + \beta_k d_{k-1}, & k \geq 1
\end{cases}

确定 $ d_{k+1} $，其中 $ \beta_{k+1} = \beta_{k+1}^{HS} $ 其中 $ \beta_{k}^{\text{HS}} = \frac{ \nabla f(x_k)^T [ \nabla f(x_k) - \nabla f(x_{k-1})] }{d_{k-1}^T [ \nabla f(x_k) - \nabla f(x_{k-1})] }, \quad k=1, 2, \ldots $ 令 $ k := k + 1 $，转步骤2。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import numpy as np

def hs_conjugate_gradient(f, grad_f, x0, epsilon=1e-6):
    """
    HS 共轭梯度法求解无约束优化问题。
    
    参数：
    f : function - 目标函数 f(x)
    grad_f : function - 目标函数的梯度 ∇f(x)
    x0 : ndarray - 初始点
    epsilon : float - 终止精度
    
    返回：
    x_opt : ndarray - 最优解
    """
    x = x0
    g = grad_f(x)  # 计算初始梯度
    d = -g  # 初始化搜索方向
    k = 0  # 迭代步数
    
    while np.linalg.norm(g) > epsilon:
        # 线搜索计算步长 α_k
        alpha_k = - (g.T @ d) / (d.T @ grad_f(x + d))
        
        # 更新 x
        x_new = x + alpha_k * d
        g_new = grad_f(x_new)  # 计算新的梯度
        
        # 计算 β_{k+1}^{HS} (Hestenes-Stiefel 公式)
        beta_hs = (g_new.T @ (g_new - g)) / (d.T @ (g_new - g))
        
        # 更新搜索方向
        d = -g_new + beta_hs * d
        
        # 更新 x 和梯度
        x, g = x_new, g_new
        k += 1
    
    return x

# 示例使用
if __name__ == "__main__":
    # 目标函数 f(x) = 0.5 * x^T G x - b^T x
    G = np.array([[4, 1], [1, 3]])  # 正定矩阵
    b = np.array([1, 2])
    
    def f(x):
        return 0.5 * x.T @ G @ x - b.T @ x
    
    def grad_f(x):
        return G @ x - b
    
    x0 = np.array([2, 1])  # 初始点
    x_opt = hs_conjugate_gradient(f, grad_f, x0)
    print("最优解:", x_opt)
1.HS公式

\beta_k^{HS} = \dfrac {[∇f(x_k)^T (∇f(x_k) - ∇f(x_{k-1}))]} {[d_{k-1}^T (∇f(x_k) - ∇f(x_{k-1}))]}

2.PRP公式2.PRP公式

\beta_k^{\text{PRP}} = \frac{
∇ f(x_k)^T [
∇ f(x_k) -
∇ f(x_{k-1})]}{|
∇ f(x_{k-1})|^2}

3.FR公式3.FR公式

\beta_k^{FR}=\dfrac {|∇f(x_k)|^2}{|∇f(x_{k-1})|^2}

4.CD公式4.CD公式

\beta_k^{CD}=-\dfrac {|∇f(x_k)|^2}{d_{k-1}^T∇f(x_{k-1})}

5.DY公式5.DY公式

\beta_k^{DY}=\dfrac {|∇f(x_{k})|^2}{d_{k-1}^T[∇f(x_k)-∇f(x_{k-1})]}

### 非线性共轭梯度法 1. **步1** 给定初始点 $x_0 \in \mathbb{R}^n$，精度 $\varepsilon > 0$。令 $d_0 = -\nabla f(x_0)$，令 $k := 0$。 2. **步2** 若 $\|\nabla f(x_k)\| \leq \varepsilon$，则得解 $x_k$，算法终止。否则转步3。 3. **步3** 由线性搜索计算步长 $\alpha_k$。 4. **步4** 令 $x_{k+1} = x_k + \alpha_k d_k$。 5. **步5** 确定 $d_{k+1}$ 和 $\beta_k$，其中 $d_{k+1} = -\nabla f(x_{k+1}) + \beta_k d_k$。令 $k := k + 1$，转步2。 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
import numpy as np
from typing import Callable, Tuple

def gradient_descent(
    f: Callable[[np.ndarray], float],          # 目标函数，输入x，返回标量函数值f(x)
    grad_f: Callable[[np.ndarray], np.ndarray], # 目标函数的梯度函数，输入x，返回梯度向量∇f(x)
    x0: np.ndarray,                            # 初始点x_0∈R^n
    epsilon: float = 1e-6,                     # 精度ε>0，用于判断梯度是否足够小以终止算法
    max_iter: int = 1000,                      # 最大迭代次数，防止无限循环
    line_search: Callable = None               # 线性搜索函数，用于计算步长α_k
) -> Tuple[np.ndarray, int]:
    """
    最速下降法（梯度下降法）实现
    
    参数:
        f: 目标函数
        grad_f: 目标函数的梯度
        x0: 初始点
        epsilon: 收敛精度
        max_iter: 最大迭代次数
        line_search: 步长搜索函数
        
    返回:
        x_k: 找到的解
        k: 迭代次数
    """
    # 步1：初始化
    x_k = x0.copy()                           # 初始点x_0
    d_k = -grad_f(x_k)                        # 初始搜索方向d_0 = -∇f(x_0)
    k = 0                                     # 迭代次数k初始化为0
    
    # 步2：迭代直到满足终止条件
    while k < max_iter:
        grad_norm = np.linalg.norm(grad_f(x_k))  # 计算当前梯度的范数||∇f(x_k)||
        
        # 检查是否满足终止条件：梯度范数是否小于等于ε
        if grad_norm <= epsilon:
            print(f"算法在 {k} 次迭代后收敛")
            return x_k, k
        
        # 步3：由线性搜索计算步长α_k
        # 如果没有提供line_search函数，则使用固定步长（这里简化为0.01，实际应用中需要调整）
        if line_search is None:
            alpha_k = 0.01  # 固定步长，实际应用应使用线搜索或其他自适应方法
        else:
            alpha_k = line_search(f, grad_f, x_k, d_k)  # 使用线搜索计算最优步长
        
        # 步4：更新x_{k+1} = x_k + α_k * d_k
        x_k_plus_1 = x_k + alpha_k * d_k
        
        # 步5：计算新的搜索方向d_{k+1} = -∇f(x_{k+1}) + β_k * d_k
        # 在最速下降法中，β_k = 0，即每次都是负梯度方向
        # 这里实现的是最速下降法（β_k=0），但保留了β_k的接口以便扩展为共轭梯度法
        beta_k = 0.0  # 对于最速下降法，β_k=0
        d_k_plus_1 = -grad_f(x_k_plus_1) + beta_k * d_k
        
        # 更新变量，准备下一次迭代
        x_k = x_k_plus_1
        d_k = d_k_plus_1
        k += 1
    
    print(f"达到最大迭代次数 {max_iter}，未收敛到所需精度")
    return x_k, k

# 示例使用：优化二次函数 f(x) = (x-1)^2 + (y-2)^2
if __name__ == "__main__":
    # 定义目标函数和梯度函数
    def f(x):
        return (x[0]-1)**2 + (x[1]-2)**2
    
    def grad_f(x):
        return np.array([2*(x[0]-1), 2*(x[1]-2)])
    
    # 初始点
    x0 = np.array([0.0, 0.0])
    # 调用梯度下降法
    solution, iterations = gradient_descent(f, grad_f, x0, epsilon=1e-6)
    
    print(f"找到的解: {solution}")
    print(f"实际最小值点: [1.0, 2.0]")
    print(f"迭代次数: {iterations}")
### Gauss-Newton法 **目的:**解决最小二乘法的问题 **步骤：** 1. 给定解的初始估计 $ x_0 $，$ \varepsilon > 0 $，$ k := 0 $； 2. 如果 $ x^{(k)} $ 满足精度要求，停止迭代； 3. 解方程组 $ J_k^T J_k d_k = -J_k^T r_k $； 4. 更新解： $ x_k := x_k + \alpha_k d_k $， 其中 $ \alpha_k $ 为一维线搜索结果，$ k := k + 1 $，转步2。 **注释：** - **基本Gauss-Newton法**：步长固定为 $ \alpha_k = 1 $； - **阻尼Gauss-Newton法**：步长 $ \alpha_k $ 通过一维线搜索动态确定。 **Gauss-Newton法的优缺点** **优点** 1. **零残量问题** 当 $ S^* = 0 $ 时，具有**局部二阶收敛速度**； 2. **小残量问题** 当 $ S^* $ 较小或问题接近线性时，具有**快速的局部收敛速度**； 3. **线性最小二乘问题** 可通过一步迭代达到极小点。 **缺点** 1. **大残量问题** 对于非严重的大残量问题，收敛速度较慢； 2. **强非线性或极大残量问题** 当残量极大或 $ f(x) $ 非线性程度过高时，算法可能**不收敛**； 3. **Jacobian矩阵不满秩** 若 $ J_k $ 不满秩，方法无法定义； 4. **全局收敛性** 算法**不一定保证总体收敛**。 **公式说明：** - $ S^* $: 残量平方和的最小值 - $ J_k $: 第 $ k $ 次迭代的 Jacobian 矩阵 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
import numpy as np

def gauss_newton(residual, jacobian, x0, epsilon, max_iter, damped=False):
    """Gauss-Newton法实现（支持阻尼选项）

    参数:
    residual -- 残差函数，输入参数向量返回残差向量
    jacobian -- 雅可比矩阵函数，输入参数向量返回雅可比矩阵
    x0 -- 初始估计值向量
    epsilon -- 收敛精度（残差范数阈值）
    max_iter -- 最大迭代次数
    damped -- 是否启用阻尼（默认False）

    返回:
    x -- 最终迭代得到的最优解向量
    """
    
    x = np.array(x0, dtype=float)  # 将初始估计值转换为numpy数组
    
    for k in range(max_iter):
        # Step 2: 检查当前解是否满足精度要求
        r = residual(x)  # 计算当前参数的残差向量
        residual_norm = np.linalg.norm(r)  # 计算残差的欧氏范数
        if residual_norm < epsilon:  # 残差足够小则停止迭代
            break
        
        # Step 3: 解线性方程组 J^T J d = -J^T r
        J = jacobian(x)  # 计算当前参数的雅可比矩阵
        JtJ = J.T @ J  # 计算J的转置与J的乘积
        Jtr = J.T @ r  # 计算J的转置与残差向量的乘积
        
        try:
            # 尝试直接解线性方程组
            d = np.linalg.solve(JtJ, -Jtr)
        except np.linalg.LinAlgError:  # 若矩阵不可逆则使用最小二乘解
            d, *_ = np.linalg.lstsq(JtJ, -Jtr, rcond=None)
        
        # Step 4: 计算步长alpha（阻尼时需进行一维搜索）
        if damped:  # 阻尼Gauss-Newton模式
            # 定义目标函数（残差平方和的一半）
            def cost(alpha):
                return 0.5 * np.linalg.norm(residual(x + alpha * d))**2
            
            alpha = 1.0  # 初始步长设为1
            c = 0.1  # Armijo条件中的充分下降系数
            rho = 0.5  # 步长缩减因子
            grad = J.T @ r  # 计算梯度∇f(x) = J^T r
            grad_dot_d = grad.dot(d)  # 梯度与方向d的点积
            
            # 检查Armijo条件：cost(alpha) ≤ cost(0) + c*alpha*grad_dot_d
            sufficient_decrease = cost(alpha) <= cost(0) + c * alpha * grad_dot_d
            while not sufficient_decrease and alpha > 1e-6:  # 当步长过小时停止搜索
                alpha *= rho  # 步长减半
                sufficient_decrease = cost(alpha) <= cost(0) + c * alpha * grad_dot_d
        else:  # 基本Gauss-Newton模式
            alpha = 1.0  # 固定步长
            
        # 更新当前解向量
        x += alpha * d  # 按步长alpha进行更新
    
    return x  # 返回最终迭代结果
### LMF方法 **步骤：** 1. **初始化** 给定 $ x_0 \in \mathbb{R}^n $，$ v_0 > 0 $，$ \varepsilon > 0 $，$ k := 0 $； 2. **终止条件** 若满足终止条件，输出结果并停止迭代； 3. **求解方程** 计算 $ (J_k^T J_k + v_k I) d_k = -J_k^T r_k $，得到步长 $ d_k $； 4. **计算步长质量** 计算 $ \rho_k = \dfrac {\triangle f_k}{\triangle q_k} $； 5. **更新阻尼因子** - 若 $ \rho_k < 0.25 $，则 $ v_{k+1} = 4v_k $； - 若 $ \rho_k > 0.75 $，则 $ v_{k+1} = v_k / 2 $； - 否则 $ v_{k+1} = v_k $； 6. **迭代更新** - 若 $ \rho_k \leq 0 $，则 $ x_{k+1} = x_k $； - 否则 $ x_{k+1} = x_k + d_k $； - 令 $ k := k + 1 $，转步骤2。 **公式说明**： - $ J_k $: 第 $ k $ 次迭代的 Jacobian 矩阵 - $ r_k $: 第 $ k $ 次迭代的残差向量 - $ v_k $: 阻尼因子（控制步长可靠性） - $ \rho_k $: 步长质量指标（评估步长接受程度） 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
import numpy as np

def levenberg_marquardt(residual, jacobian, x0, epsilon, max_iter, v0=0.01):
    """Levenberg-Marquardt算法实现（非线性最小二乘问题）

    参数:
    residual -- 残差函数，输入参数向量返回残差向量
    jacobian -- 雅可比矩阵函数，输入参数向量返回雅可比矩阵
    x0 -- 初始估计值向量
    epsilon -- 收敛精度（残差范数阈值）
    max_iter -- 最大迭代次数
    v0 -- 初始阻尼因子（默认0.01）

    返回:
    x -- 最终迭代得到的最优解向量
    """
    
    x_prev = np.array(x0, dtype=float)  # 初始参数估计
    v = v0  # 初始阻尼因子
    k = 0  # 迭代计数器
    
    while k < max_iter:
        # Step 2: 检查终止条件
        r = residual(x_prev)  # 计算当前参数的残差向量
        residual_norm = np.linalg.norm(r)  # 计算残差的欧氏范数
        if residual_norm < epsilon:  # 残差足够小则停止迭代
            break
        
        # Step3: 解方程 (J^T J + vI) d = -J^T r
        J = jacobian(x_prev)  # 当前参数的雅可比矩阵
        JtJ = J.T @ J  # 计算J的转置与J的乘积
        Jtr = J.T @ r  # 计算J的转置与残差的乘积
        
        # 构建矩阵A = JtJ + vI
        A = JtJ + v * np.eye(len(x_prev))
        
        try:
            # 尝试直接解线性方程组
            d = np.linalg.solve(A, -Jtr)
        except np.linalg.LinAlgError:  # 若矩阵奇异则使用最小二乘解
            d, *_ = np.linalg.lstsq(A, -Jtr, rcond=None)
        
        # 计算候选解x_candidate和新残差
        x_candidate = x_prev + d
        r_new = residual(x_candidate)
        
        # Step4: 计算步长质量指标ρ
        # 计算实际减少量Δf = 0.5*(||r||^2 - ||r_new||^2)
        delta_f = 0.5 * (np.linalg.norm(r)**2 - np.linalg.norm(r_new)**2)
        
        # 计算预测减少量Δq = d^T (JtJ d) + v*(d^T d)
        JtJd = JtJ @ d
        delta_q = d.dot(JtJd) + v * (d.dot(d))
        
        # 当预测减少量为负时强制拒绝步长
        if delta_q <= 0:
            rho = -np.inf  # 强制ρ为负无穷
        else:
            rho = delta_f / delta_q  # 计算步长质量
        
        # Step5: 更新阻尼因子v
        if rho < 0.25:  # 步长质量差，增大阻尼因子
            v_new = v * 4.0
        elif rho > 0.75:  # 步长质量好，减小阻尼因子
            v_new = v / 2.0
        else:  # 保持当前阻尼因子
            v_new = v
        
        # Step6: 根据ρ值决定是否接受步长
        if rho <= 0:  # 步长导致残差增大，不接受更新
            x_next = x_prev.copy()  # 保持当前参数
            v = v_new  # 更新阻尼因子
        else:  # 接受步长，更新参数和阻尼因子
            x_next = x_candidate
            v = v_new
        
        # 更新迭代变量
        x_prev = x_next
        k += 1
    
    return x_prev  # 返回最终解
### 外点罚函数法 **步0**：选定初始点 $x_0 \in \mathbb{R}^n$；初始罚因子 $\mu_1 > 0$，放大系数 $\sigma > 1$，精度 $\varepsilon > 0$。置 $k=1$。 **步1**：构造增广目标函数 $$ F_{\mu_k}(x) = f(x) + \dfrac 1 2\mu_k P(x)

步2：以 xk−1x_{k-1} 为初始点求解无约束问题：

min⁡Fμk(x),得解 xk\min F_{\mu_k}(x), \quad \text{得解} \ x_k

步3：若 μkP(xk)<ε\mu_k P(x_k) < \varepsilon，停止迭代，输出 xkx_k；否则转步4。

步4：更新罚因子 μk+1=σμk\mu_{k+1} = \sigma \mu_k，置 k=k+1k = k + 1，转步1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np

# 外点罚函数法
def exterior_penalty_method(f, P, x0, mu1, sigma, epsilon, solver):
    """
    参数说明：
    f: 原目标函数 f(x)
    P: 罚函数 P(x)
    x0: 初始点，numpy数组
    mu1: 初始罚因子，正数
    sigma: 放大系数，sigma > 1
    epsilon: 精度要求
    solver: 无约束优化子程序，输入一个目标函数和初始点，输出局部极小点
    """

    # 步0：初始化
    xk = np.array(x0)
    muk = mu1
    k = 1

    while True:
        # 步1：构造增广目标函数
        def F_mu(x):
            x = np.array(x)
            return f(x) + (1 / (2 * muk)) * P(x)

        # 步2：以 x_{k-1} 为初始点求解无约束问题
        xk = solver(F_mu, xk)

        # 步3：检查停止条件
        if muk * P(xk) < epsilon:
            return xk

        # 步4：更新罚因子和迭代次数
        muk = sigma * muk
        k += 1

内点罚函数法

步0：选定初始点 x0∈Rnx_0 \in \mathbb{R}^n；初始罚因子 μ1>0\mu_1 > 0，放大系数 σ>1\sigma > 1，精度 ε>0\varepsilon > 0。置 k=1k=1。

步1：构造增广目标函数

Fμk(x)=f(x)+μk−1S(x)F_{\mu_k}(x) = f(x) + \mu_k^{-1} S(x)

步2：以 xk−1x_{k-1} 为初始点求解无约束问题：

min⁡Fμk(x),得解 xk\min F_{\mu_k}(x), \quad \text{得解} \ x_k

步3：若 μkS(xk)<ε\mu_k S(x_k) < \varepsilon，停止迭代，输出 xkx_k；否则转步4。

步4：更新罚因子 μk+1=σμk\mu_{k+1} = \sigma \mu_k，置 k=k+1k = k + 1，转步1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import numpy as np

# 内点罚函数法
def interior_penalty_method(f, S, x0, mu1, sigma, epsilon, solver):
    """
    参数说明：
    f: 原目标函数 f(x)
    S: 罚函数 S(x)（比如 -log形式）
    x0: 初始点，numpy数组
    mu1: 初始罚因子，正数
    sigma: 放大系数，sigma > 1
    epsilon: 精度要求
    solver: 无约束优化子程序，输入一个目标函数和初始点，输出局部极小点
    """

    # 步0：初始化
    xk = np.array(x0)
    muk = mu1
    k = 1

    while True:
        # 步1：构造增广目标函数
        def F_mu(x):
            x = np.array(x)
            return f(x) + (1 / muk) * S(x)

        # 步2：以 x_{k-1} 为初始点求解无约束问题
        xk = solver(F_mu, xk)

        # 步3：检查停止条件
        if muk * S(xk) < epsilon:
            return xk

        # 步4：更新罚因子和迭代次数
        muk = sigma * muk
        k += 1

等式约束问题的乘子法

步0：选定初始点 x0∈Rnx_0 \in \mathbb{R}^n，初始乘子估计 λ1∈Rm\lambda_1 \in \mathbb{R}^m；初始罚因子 μ1>0\mu_1 > 0，常数 σ>1\sigma > 1，β∈(0,1)\beta \in (0,1)；精度 ε>0\varepsilon > 0。置 k:=1k := 1。

步1：构造增广目标函数

Lμk(x,λk)=L(x,λk)+12μkS(x)L_{\mu_k}(x, \lambda_k) = L(x, \lambda_k) + \frac{1}{2} \mu_k S(x)

其中 L(x,λk)=f(x)−λkTh(x)L(x, \lambda_k) = f(x) - \lambda_k^T h(x)。

步2：以 xk−1x_{k-1} 为初始点求解无约束问题：

min⁡xLμk(x,λk)\min_x L_{\mu_k}(x, \lambda_k)

设其解为 xkx_k。

步3：若 ∥h(xk)∥≤ε\|h(x_k)\| \leq \varepsilon，则得解 xkx_k，STOP；否则转步4。

步4：若 ∥h(xk)∥∥h(xk−1)∥≥β\frac{\|h(x_k)\|}{\|h(x_{k-1})\|} \geq \beta 成立，则令 μk+1=μk\mu_{k+1} = \mu_k；否则 μk+1=σμk\mu_{k+1} = \sigma \mu_k。

步5：令 λk+1=λk−μkh(xk)\lambda_{k+1} = \lambda_k - \mu_k h(x_k)，置 k:=k+1k := k+1，转步1。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import numpy as np

# 等式约束问题的乘子法
def multiplier_method(f, h, x0, lambda0, mu1, sigma, beta, epsilon, solver):
    """
    参数说明：
    f: 目标函数 f(x)
    h: 等式约束函数 h(x)，返回一个向量
    x0: 初始点，numpy数组
    lambda0: 初始乘子估计，numpy数组
    mu1: 初始罚因子，正数
    sigma: 放大系数，sigma > 1
    beta: 收敛判别系数，0 < beta < 1
    epsilon: 精度要求
    solver: 无约束优化子程序，输入一个目标函数和初始点，输出局部极小点
    """

    # 步0：初始化
    xk = np.array(x0)
    lambdak = np.array(lambda0)
    muk = mu1
    k = 1

    while True:
        # 步1：构造增广拉格朗日函数
        def L_mu(x):
            x = np.array(x)
            L = f(x) - np.dot(lambdak, h(x))
            penalty = (1/2) * muk * np.linalg.norm(h(x))**2
            return L + penalty

        # 步2：以 x_{k-1} 为初始点求解无约束极小化问题
        x_prev = xk.copy()
        xk = solver(L_mu, xk)

        # 步3：检查停止条件
        if np.linalg.norm(h(xk)) <= epsilon:
            return xk

        # 步4：更新罚因子
        if np.linalg.norm(h(xk)) / np.linalg.norm(h(x_prev)) >= beta:
            muk = muk  # 保持不变
        else:
            muk = sigma * muk

        # 步5：更新拉格朗日乘子并迭代
        lambdak = lambdak - muk * h(xk)
        k += 1

一般约束问题的乘子法

步0：初始化参数

  • 初始点 x0∈Rnx_0 \in \mathbb{R}^n
  • 乘子估计 λ(1)∈Rm\lambda^{(1)} \in \mathbb{R}^m（需满足 λi(1)≥0, ∀i∈I\lambda_i^{(1)} \geq 0,\ \forall i \in I）
  • 初始罚因子 μ1>0\mu_1 > 0，放大系数 σ>1\sigma > 1，收缩率 β∈(0,1)\beta \in (0,1)
  • 精度阈值 ε>0\varepsilon > 0
  • 置迭代计数器 k=1k = 1

步1：构造增广Lagrange函数

Lμ(x,λ)=f(x)+12μ∑i∈I(max⁡2{0,λi−μgi(x)}−λi2)−∑j∈Eλjhj(x)+μ2∑j∈Ehj2(x)L_{\mu}(x,\lambda) = f(x) + \frac{1}{2\mu} \sum_{i \in I} \left( \max^2 \{0, \lambda_i - \mu g_i(x)\} - \lambda_i^2 \right) - \sum_{j \in E} \lambda_j h_j(x) + \frac{\mu}{2} \sum_{j \in E} h_j^2(x)

步2：求解无约束子问题
以 xk−1x_{k-1} 为初始点，求解：

min⁡xLμk(x,λ(k))\min_x L_{\mu_k}(x, \lambda^{(k)})

得到解 xkx_k

步3：终止判定
若 ∥h(xk)∥+∥min⁡(g(xk),0)∥≤ε\|h(x_k)\| + \|\min(g(x_k), 0)\| \leq \varepsilon，输出 xkx_k 并终止；否则转步4

步4：自适应罚因子更新

若 ∥h(xk)∥+∥min⁡(g(xk),0)∥∥h(xk−1)∥+∥min⁡(g(xk−1),0)∥≥β, 则 μk+1=μk\text{若}\ \frac{\|h(x_k)\| + \|\min(g(x_k),0)\|}{\|h(x_{k-1})\| + \|\min(g(x_{k-1}),0)\|} \geq \beta,\ \text{则}\ \mu_{k+1} = \mu_k

否则 μk+1=σμk\mu_{k+1} = \sigma \mu_k

步5：乘子更新

{λi(k+1)=max⁡{0, λi(k)−μgi(xk)},∀i∈Iλj(k+1)=λj(k)−μhj(xk),∀j∈E\begin{cases} \lambda_i^{(k+1)} = \max \left\{ 0,\ \lambda_i^{(k)} - \mu g_i(x_k) \right\}, & \forall i \in I \\ \lambda_j^{(k+1)} = \lambda_j^{(k)} - \mu h_j(x_k), & \forall j \in E \end{cases}

置 k:=k+1k := k+1，返回步1

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
import numpy as np

# 一般约束问题的乘子法
def general_multiplier_method(f, g, h, x0, lambda0, mu1, sigma, beta, epsilon, solver):
    """
    参数说明：
    f: 目标函数 f(x)
    g: 不等式约束函数组 g(x)，返回一个numpy数组，g_i(x) <= 0
    h: 等式约束函数组 h(x)，返回一个numpy数组，h_j(x) = 0
    x0: 初始点，numpy数组
    lambda0: 初始乘子估计（先 g 后 h 排列），numpy数组
    mu1: 初始罚因子，正数
    sigma: 放大系数，sigma > 1
    beta: 收缩率，0 < beta < 1
    epsilon: 精度阈值
    solver: 无约束优化子程序，输入一个目标函数和初始点，输出局部极小点
    """

    # 步0：初始化
    xk = np.array(x0)
    lambdak = np.array(lambda0)
    muk = mu1
    k = 1

    # 记录不等式与等式的个数
    m_g = len(g(xk))  # 不等式约束数目
    m_h = len(h(xk))  # 等式约束数目

    while True:
        # 步1：构造增广Lagrange函数
        def L_mu(x):
            x = np.array(x)
            gx = g(x)
            hx = h(x)
            term1 = f(x)
            term2 = (1 / (2 * muk)) * np.sum((np.maximum(0, lambdak[:m_g] - muk * gx))**2 - lambdak[:m_g]**2)
            term3 = - np.dot(lambdak[m_g:], hx)
            term4 = (muk / 2) * np.sum(hx**2)
            return term1 + term2 + term3 + term4

        # 步2：求解无约束子问题
        x_prev = xk.copy()
        xk = solver(L_mu, xk)

        # 步3：终止判定
        violation = np.linalg.norm(h(xk)) + np.linalg.norm(np.minimum(g(xk), 0))
        if violation <= epsilon:
            return xk

        # 步4：自适应罚因子更新
        prev_violation = np.linalg.norm(h(x_prev)) + np.linalg.norm(np.minimum(g(x_prev), 0))
        if violation / prev_violation >= beta:
            muk = muk  # 保持罚因子不变
        else:
            muk = sigma * muk

        # 步5：乘子更新
        gx = g(xk)
        hx = h(xk)
        lambda_g_new = np.maximum(0, lambdak[:m_g] - muk * gx)
        lambda_h_new = lambdak[m_g:] - muk * hx
        lambdak = np.concatenate([lambda_g_new, lambda_h_new])

        # 迭代次数加一
        k += 1