基于近似Hessian的序列二次规划算法（SQP with Approximate Hessian）基础题
题目描述：
考虑一个带有非线性等式约束的优化问题：
\[\min_{x \in \mathbb{R}^n} f(x) \quad \text{s.t.} \quad c(x) = 0 \]
其中，目标函数 \(f: \mathbb{R}^n \to \mathbb{R}\) 是二阶连续可微的，约束函数 \(c: \mathbb{R}^n \to \mathbb{R}^m\) 也是二阶连续可微的，且 \(m \leq n\)。假设计算精确的Hessian矩阵（即拉格朗日函数的二阶导数）计算代价很高，我们希望在序列二次规划（SQP）框架下，使用一个近似矩阵 \(B_k\) 来代替精确的Hessian，构建二次规划子问题，并迭代求解。请详细解释如何构造这样的近似Hessian SQP算法，包括子问题的形成、迭代格式、以及确保算法收敛的关键技巧（如BFGS更新）。
解题过程循序渐进讲解：
1. 问题背景与SQP基本思想回顾
序列二次规划（SQP）是求解非线性规划问题的一类重要算法。其核心思想是：在当前迭代点 \(x_k\) 处，利用目标函数和约束的一阶、二阶信息，构造一个二次规划（QP）子问题，求解该子问题得到搜索方向 \(p_k\)，然后沿该方向进行线性搜索或信赖域步长调整，得到下一个迭代点 \(x_{k+1}\)。
对于等式约束问题，在点 \(x_k\) 处，精确的SQP子问题为：
\[\begin{aligned} \min_{p \in \mathbb{R}^n} & \quad \nabla f(x_k)^T p + \frac{1}{2} p^T \nabla^2_{xx} L(x_k, \lambda_k) p \\ \text{s.t.} & \quad c(x_k) + A(x_k)^T p = 0 \end{aligned} \]
其中，\(L(x, \lambda) = f(x) + \lambda^T c(x)\) 是拉格朗日函数，\(A(x) = \nabla c(x) \in \mathbb{R}^{n \times m}\) 是约束的雅可比矩阵（每列是一个约束的梯度），\(\nabla^2_{xx} L(x_k, \lambda_k)\) 是拉格朗日函数关于 \(x\) 的Hessian矩阵。
难点：精确计算 \(\nabla^2_{xx} L(x_k, \lambda_k)\) 可能很昂贵（例如，需要二阶导数，或问题维度很高）。因此，我们想用一个近似矩阵 \(B_k \approx \nabla^2_{xx} L(x_k, \lambda_k)\) 来代替它。
2. 近似Hessian SQP的子问题构造
在迭代点 \(x_k\)，给定当前的拉格朗日乘子估计 \(\lambda_k\) 和一个对称矩阵 \(B_k\)，我们构造如下QP子问题：
\[\begin{aligned} \min_{p \in \mathbb{R}^n} & \quad \nabla f(x_k)^T p + \frac{1}{2} p^T B_k p \\ \text{s.t.} & \quad c(x_k) + A(x_k)^T p = 0 \end{aligned} \tag{1} \]
这里，我们直接用 \(B_k\) 取代了精确的Hessian。子问题(1)是一个等式约束的二次规划，其解 \(p_k\) 可以作为搜索方向。
为什么这样近似是合理的？
- 在最优解 \(x^*\) 处，如果 \(B_k\) 收敛到 \(\nabla^2_{xx} L(x^*, \lambda^*)\)，那么SQP迭代的局部收敛速度可以是超线性的（例如，采用拟牛顿更新）。
- 即使 \(B_k\) 只是粗略近似，只要它保持一定的正定性（或至少子空间曲率正确），算法仍可能全局收敛。
3. 求解QP子问题与迭代格式
子问题(1)的KKT条件为：
\[\begin{cases} B_k p + \nabla f(x_k) + A(x_k) \lambda_{QP} = 0 \\ c(x_k) + A(x_k)^T p = 0 \end{cases} \tag{2} \]
其中 \(\lambda_{QP}\) 是该QP子问题的乘子。我们可以将(2)写成线性系统：
\[\begin{bmatrix} B_k & A(x_k) \\ A(x_k)^T & 0 \end{bmatrix} \begin{bmatrix} p \\ \lambda_{QP} \end{bmatrix} = - \begin{bmatrix} \nabla f(x_k) \\ c(x_k) \end{bmatrix} \tag{3} \]
只要矩阵非奇异，就可解出 \(p_k\) 和 \(\lambda_{QP}\)。
然后，SQP的迭代格式为：
- 更新迭代点：\(x_{k+1} = x_k + \alpha_k p_k\)，其中 \(\alpha_k \in (0, 1]\) 是步长，通常通过线搜索（如满足某些价值函数下降）确定。
- 更新乘子：\(\lambda_{k+1} = \lambda_{QP}\)（即直接采用QP子问题的乘子作为新乘子估计）。
- 更新近似Hessian：\(B_{k+1} = \text{Update}(B_k, s_k, y_k)\)，其中 \(s_k = x_{k+1} - x_k\)，\(y_k\) 是与拉格朗日函数梯度差相关的向量，具体见下文。
4. 近似Hessian的更新：BFGS公式
为了不计算精确Hessian，我们采用拟牛顿法更新 \(B_k\)。对于无约束问题，BFGS更新利用梯度信息来近似Hessian。对于约束问题，我们需要近似的是拉格朗日函数的Hessian \(\nabla^2_{xx} L(x, \lambda)\)。
定义：
\[s_k = x_{k+1} - x_k, \quad y_k = \nabla_x L(x_{k+1}, \lambda_{k+1}) - \nabla_x L(x_k, \lambda_{k+1}) \]
其中，\(\nabla_x L(x, \lambda) = \nabla f(x) + A(x) \lambda\)。注意这里第二个梯度计算时使用新的乘子 \(\lambda_{k+1}\)，这是为了保证拟牛顿方程 \(B_{k+1} s_k = y_k\) 在拉格朗日函数梯度上近似成立。
BFGS更新公式为（保持 \(B_k\) 对称正定）：
\[B_{k+1} = B_k - \frac{B_k s_k s_k^T B_k}{s_k^T B_k s_k} + \frac{y_k y_k^T}{s_k^T y_k} \]
这个更新要求 \(s_k^T y_k > 0\)（曲率条件），在约束问题中，即使目标函数非凸，只要 \(B_k\) 初始正定且线搜索合适，该条件通常可满足。
为什么用 \(\lambda_{k+1}\) 而不是 \(\lambda_k\) 计算 \(y_k\)？
这样可以确保在解附近，\(y_k\) 近似等于 \(\nabla^2_{xx} L(\xi, \lambda_{k+1}) s_k\)（对某个 \(\xi\)），从而 \(B_{k+1}\) 更好地近似当前乘子下的Hessian，提高收敛速度。
5. 全局化策略：价值函数与线搜索
纯牛顿型迭代（步长 \(\alpha_k=1\)）可能只在初始点靠近解时收敛。为了全局收敛，我们引入价值函数（merit function）来衡量进展，并据此选择步长 \(\alpha_k\)。
常用的价值函数是 \(l_1\) 精确罚函数：
\[\Phi_\sigma (x) = f(x) + \sigma \| c(x) \|_1 \]
其中 \(\sigma > 0\) 是罚参数。我们要求搜索方向 \(p_k\) 是 \(\Phi_\sigma\) 的下降方向，即 \(D(p_k) := \nabla f(x_k)^T p_k - \sigma \| c(x_k) \|_1 < 0\)。通过调整 \(\sigma\) 足够大，可以确保这一点。
线搜索（如Armijo条件）：选择 \(\alpha_k\) 为 \(1, \beta, \beta^2, \dots\) 中最大的那个，使得
\[\Phi_\sigma (x_k + \alpha p_k) \leq \Phi_\sigma (x_k) + \eta \alpha D(p_k) \]
其中 \(\beta \in (0,1)\)，\(\eta \in (0, 0.5)\)。
6. 算法完整步骤
- 初始化：选择初始点 \(x_0\)，初始乘子 \(\lambda_0\)，初始正定矩阵 \(B_0\)（例如单位矩阵），罚参数 \(\sigma_0 > 0\)，参数 \(\beta \in (0,1)\)，\(\eta \in (0,0.5)\)，收敛容差 \(\epsilon > 0\)。设 \(k=0\)。
- 收敛检查：如果 \(\| \nabla_x L(x_k, \lambda_k) \| \leq \epsilon\) 且 \(\| c(x_k) \| \leq \epsilon\)，则停止。
- 求解QP子问题：求解(1)得到 \(p_k\) 和 \(\lambda_{QP}\)。
- 价值函数下降测试：计算 \(D(p_k) = \nabla f(x_k)^T p_k - \sigma_k \| c(x_k) \|_1\)。如果 \(D(p_k) \ge 0\)，则增大 \(\sigma_k\) 直到 \(D(p_k) < 0\)。
- 线搜索：用Armijo规则在价值函数 \(\Phi_{\sigma_k}\) 上搜索步长 \(\alpha_k\)。
- 更新迭代点：\(x_{k+1} = x_k + \alpha_k p_k\)。
- 更新乘子：\(\lambda_{k+1} = \lambda_{QP}\)。
- 更新近似Hessian：计算 \(s_k = x_{k+1} - x_k\)，\(y_k = \nabla_x L(x_{k+1}, \lambda_{k+1}) - \nabla_x L(x_k, \lambda_{k+1})\)。如果 \(s_k^T y_k > 0\)，用BFGS公式更新 \(B_{k+1}\)；否则跳过更新（或采用阻尼BFGS）。
- 更新罚参数：根据需要调整 \(\sigma_{k+1}\)（例如，若约束违反未充分下降，则增大）。
- 设 \(k = k+1\)，返回步骤2。
7. 关键点与注意事项
- 近似Hessian的正定性：BFGS更新保持正定，这有助于QP子问题有唯一解。如果 \(s_k^T y_k\) 不满足正曲率条件，可采用阻尼BFGS等技术。
- 局部超线性收敛：在适当条件下（如 \(B_k\) 收敛到精确Hessian，且步长最终取1），算法局部超线性收敛。
- 与精确Hessian SQP的关系：当 \(B_k = \nabla^2_{xx} L(x_k, \lambda_k)\) 时，就是标准SQP；这里用拟牛顿近似避免了二阶导数计算，适合大规模或黑箱函数问题。
- 实现细节：QP子问题求解可用积极集法或内点法；矩阵 \(B_k\) 通常以紧凑形式（如L-BFGS）存储以节省内存。
总结：基于近似Hessian的SQP算法通过用拟牛顿矩阵替代精确二阶导数，在保持快速局部收敛的同时，大幅降低了每步计算量。其核心是BFGS更新、价值函数线搜索的全局化，以及乘子的自然更新。该框架可扩展至不等式约束（通过引入松弛变量或积极集策略），是实际中应用广泛的非线性规划方法。