【学习笔记】最优化方法

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

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

$$f(x_k + \alpha_k d_k) \leq f(x_k) + \rho \alpha_k \nabla f(x_k)^T d_k$$

参数说明：

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

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

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

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

---

参数说明：

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

算法注释：

1. 先尝试 $\alpha_k=1$，若满足下降条件（2.6），则直接采用单位步长
2. 若不满足，则从较大步长 $\beta$ 开始，通过循环逐步缩小步长（$\alpha_k \leftarrow \rho\alpha_k$）
3. 适用于梯度下降法、拟牛顿法等优化算法，保证迭代过程中目标函数值 $f(x_k)$ 充分下降

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(x_k + \alpha_k d_k) \leq f(x_k) + \rho \alpha_k \nabla f(x_k)^T d_k$$

2. 曲率条件：

$$\nabla f(x_k + \alpha_k d_k)^T d_k \geq \sigma \nabla f(x_k)^T d_k$$

参数约束：

$$0 < \rho < \sigma < 1$$

---

参数说明：

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

---

公式作用：

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

---

补充说明：

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

---

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

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

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

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

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

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

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

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

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

---

参数说明：

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

---

解题思路与算法逻辑：

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

---

补充说明：

1. 条件 (2.7) 应包含两个不等式（推测为类似 Wolfe-Powell 条件的 Armijo 条件和曲率条件）
2. 该算法比标准 Armijo 搜索更复杂，适用于需要同时满足多个条件的优化问题
3. 离散化集合 $\{ \beta \rho^j \}$ 和 $\{ \alpha_k^{(i)} + \rho_1^j (\beta_k^{(i)} - \alpha_k^{(i)}) \}$ 的设计平衡了计算效率与精度

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  # 继续搜索

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

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

1. 初始化：

   • 计算初始内点：

     $$\lambda_1 = a_1 + 0.382(b_1 - a_1)$$

     $$\mu_1 = a_1 + 0.618(b_1 - a_1)$$

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

2. 停止条件：

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

3. 区间更新：

   • 情况1：若 $f(\lambda_k) > f(\mu_k)$

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

   • 情况2：若 $f(\lambda_k) \leq f(\mu_k)$

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

4. 迭代循环：

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

---

参数说明：

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

---

算法原理：

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

---

解题思路：

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

2. 关键步骤验证：

      # 伪代码示例：判断如何更新区间
      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  # 返回最终区间

信赖域算法

修正后的规范算法步骤：

输入：初始点 $x_0$，最大信赖域半径 $\overline{\Delta}$，精度阈值 $\varepsilon \geq 0$，阈值参数 $0 < \eta_1 \leq \eta_2 < 1$
输出：优化解 $x_k$

1. 初始化：

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

2. 停止条件：

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

3. 子问题求解：

   • 在信赖域 $\|s\| \leq \Delta_k$ 内，近似求解：

     $$\min_{s} m_k(s) = f(x_k) + g_k^T s + \frac{1}{2}s^T B_k s$$

     得到试探步 $s_k$，其中 $B_k$ 是Hessian矩阵或其近似

4. 接受步长：

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

     $$r_k = \frac{\delta f}{\delta m}$$

   • 更新迭代点：
     $$x_{k+1} = \begin{cases} x_k + s_k, & \text{if } r_k > 0 \\ x_k, & \text{otherwise} \end{cases}$$

5. 调整信赖域半径：
   $$\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+1$，返回步骤2

---

参数说明：

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

---

算法原理：

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

---

解题思路：

1. 关键公式验证：

      # 伪代码示例：计算接受率
      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. 方向计算：

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

     $$d_k = -\nabla f(x_k)$$

4. 步长计算：

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

     $$\alpha_k = \arg\min_{\alpha > 0} f(x_k + \alpha d_k)$$

5. 迭代更新：

   • 更新迭代点：

     $$x_{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. 梯度计算：

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

2. 方向更新：

   $$d_k = - (Qx_k + b)$$

3. 精确步长解析解：

   $$\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条件降低计算成本

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. 方向计算：

   • 解线性方程组：

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

   • 得到搜索方向 $d_k$

4. 迭代更新：

   • 更新迭代点：

     $$x_{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连续且正定，则

     $$\|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：

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

2. 牛顿方向计算：

   $$Q d_k + (Qx_k + b) = 0 \Rightarrow d_k = -x_k - Q^{-1}b$$

3. 迭代更新：

   $$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}$（精确解）

---

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. 方向计算：

   • 解线性方程组：

     $$\nabla^2 f(x_k) d_k + \nabla f(x_k) = 0$$

   • 得到搜索方向 $d_k$

4. 步长计算：

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

5. 迭代更新：

   • 更新迭代点：

     $$x_{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：

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

  2. 牛顿方向：

     $$Q d_k + (Qx_k + b) = 0 \Rightarrow d_k = -x_k - Q^{-1}b$$

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

     $$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}$（精确解）

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：给定初始点 $x_0 \in \mathbb{R}^n$，精度 $\epsilon > 0$。令 $k = 0$；

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

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

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

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

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

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$，算法终止；否则，解线性方程组

$$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$ 的上界，即要求

   $$v_k \leq C \|\nabla f(x_k)\|$$

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

修正形式的收敛性

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

鞍点情况的特殊处理

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

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

   $$d^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 求解线性方程组：

$$B_k d + \nabla f(x_k) = 0$$

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

$$d_k = -H_k g_k \quad (\text{其中 } g_k = \nabla f(x_k))$$

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

$$x_{k+1} = x_k + \alpha_k d_k$$

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

$$B_{k+1} = B_k + \Delta B \quad \text{或} \quad H_{k+1} = H_k + \Delta H$$

使其满足拟Newton方程 $B_{k+1}s_k$ 或 $H_{k+1}y_k=s_k$ 。令 $k := k + 1$，返回步骤2。

修正矩阵

对称秩1

SR修正公式形式1

$$B_{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

$$H_{k+1}=H_k+ \dfrac {(s_k-H_ky_k)(s_k-H_ky_k)^T} {(s_k-H_ky_k)^Ty_k}$$

对称秩2

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

$$B_{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}$$

$BFGS$ 的逆修正公式

$$H_{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}$$

$DFP$ 修正公式

$$H_{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}$$

$DFP$ 的逆修正公式

$$B_{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：给定初始点 $x_0 \in \mathbb{R}^n$，初始对称正定矩阵 $B_0$，精度 $\epsilon > 0$。令 $k = 0$。

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

步骤3：解线性方程组

$$B_k d + \nabla f(x_k) = 0 \qquad$$

得解 $d_k$。

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

步骤5：令 $x_{k+1} = x_k + \alpha_k d_k$，若 $\| \nabla f(x_{k+1}) \| \leq \epsilon$，则得解 $x_{k+1}$，算法终止。否则由公式 $B_{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}$ 计算 $B_{k+1}$。

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

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：给定初始点 $x_0 \in \mathbb{R}^n$，初始对称正定矩阵 $H_0$，精度 $\epsilon > 0$。令 $k = 0$。

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

步骤3：计算搜索方向：

$$d_k = -H_k \nabla f(x_k)$$

（或通过解线性方程组 $B_k d + \nabla f(x_k) = 0$，其中 $B_k$ 由公式 $B_{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：通过线性搜索确定步长 $\alpha_k$。

步骤5：更新迭代点：

$$x_{k+1} = x_k + \alpha_k d_k$$

若 $\| \nabla f(x_{k+1}) \| \leq \epsilon$，输出解 $x_{k+1}$ 并终止；否则按公式 $H_{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}$ 计算 $H_{k+1}$。

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

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族

$$B_{k+1}^{\theta_k}=\theta_kB_{k+1}^{BFGS}+(1-\theta_k)B_{k+1}^{DFP}$$

$$H_{k+1}^{\phi_k}=\phi_kH_{k+1}^{BFGS}+(1-\phi_k)H_{k+1}^{DFP}$$

其中 $\theta_k=1-\dfrac {1-\phi_k}{1-\phi_k+\phi_k \mu_k}$ ， $\mu_k=\dfrac{y_k^TH_ky_ks_k^TB_ks_k}{s_k^Ty_k}$

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

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

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

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

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

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

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

一个 $m-n_1-n_2-q$ 前馈神经网络，表示该网络有 $m$ 个源节点， 第一个隐藏层有 $n_1$ 个神经元，第二个隐藏层有 $n2$ 个神经元，输出层有 $q$ 个神经元

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

考虑一阶常微分方程

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

其初始条件为

$$y(0)=y_0$$

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

$$y_t(x,p)$$

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

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

$$y_t(x,p)=y_0+xN(x,p)$$

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

$$N(x,p)= \sum_{j=1}^{n} v_jφ(z_j),z_j=w_jx-\theta_j$$

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

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

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

$$\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$$

来确定 $p$ 的值

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

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

$$\mathbf{d}_0 = -\mathbf{g}_0, \quad k := 0$$

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

$$\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$，使得：

$$\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

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：计算

$$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。

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。

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公式

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

3.FR公式

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

4.CD公式

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

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。

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 矩阵

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$: 步长质量指标（评估步长接受程度）

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：以 $x_{k-1}$ 为初始点求解无约束问题：

$$\min F_{\mu_k}(x), \quad \text{得解} \ x_k$$

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

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

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：选定初始点 $x_0 \in \mathbb{R}^n$；初始罚因子 $\mu_1 > 0$，放大系数 $\sigma > 1$，精度 $\varepsilon > 0$。置 $k=1$。

步1：构造增广目标函数

$$F_{\mu_k}(x) = f(x) + \mu_k^{-1} S(x)$$

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

$$\min F_{\mu_k}(x), \quad \text{得解} \ x_k$$

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

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

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：选定初始点 $x_0 \in \mathbb{R}^n$，初始乘子估计 $\lambda_1 \in \mathbb{R}^m$；初始罚因子 $\mu_1 > 0$，常数 $\sigma > 1$，$\beta \in (0,1)$；精度 $\varepsilon > 0$。置 $k := 1$。

步1：构造增广目标函数

$$L_{\mu_k}(x, \lambda_k) = L(x, \lambda_k) + \frac{1}{2} \mu_k S(x)$$

其中 $L(x, \lambda_k) = f(x) - \lambda_k^T h(x)$。

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

$$\min_x L_{\mu_k}(x, \lambda_k)$$

设其解为 $x_k$。

步3：若 $\|h(x_k)\| \leq \varepsilon$，则得解 $x_k$，STOP；否则转步4。

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

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

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：初始化参数

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

步1：构造增广Lagrange函数

$$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：求解无约束子问题
以 $x_{k-1}$ 为初始点，求解：

$$\min_x L_{\mu_k}(x, \lambda^{(k)})$$

得到解 $x_k$

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

步4：自适应罚因子更新

$$\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$$

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

步5：乘子更新

$$\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+1$，返回步1

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