非线性规划中的拟牛顿法（BFGS）基础题
题目描述：
考虑一个无约束非线性规划问题：
\[\min f(x) = (x_1 - 1)^2 + 2(x_2 - 2)^2 + 3(x_3 - 3)^2 \]
其中 \(x = (x_1, x_2, x_3)^T \in \mathbb{R}^3\)。该目标函数是二次可微的正定二次函数，具有唯一全局最小值点 \(x^* = (1, 2, 3)^T\)，对应最优值 \(f(x^*) = 0\)。要求使用拟牛顿法中的BFGS算法求解该问题，从初始点 \(x^{(0)} = (0, 0, 0)^T\) 开始，初始近似Hessian矩阵取为单位矩阵 \(B_0 = I\)，并详细展示BFGS迭代的完整过程（包括梯度计算、步长选择、矩阵更新等），直到满足梯度范数小于 \(10^{-3}\) 的收敛条件。
解题过程循序渐进讲解：
1. 背景与算法思想
拟牛顿法（Quasi-Newton Method）是求解无约束优化问题的一类迭代算法，旨在模拟牛顿法的快速收敛性，同时避免计算精确Hessian矩阵的高成本。BFGS（Broyden–Fletcher–Goldfarb–Shanno）算法是其中最著名的一种，它通过迭代更新一个近似Hessian矩阵（或其逆矩阵），使得该矩阵满足拟牛顿条件（也称为割线方程）。
核心思想：在第 \(k\) 次迭代，用当前近似Hessian矩阵 \(B_k\) 代替牛顿法中的精确Hessian，求解线性方程组 \(B_k p_k = -\nabla f(x^{(k)})\) 得到搜索方向 \(p_k\)，然后沿 \(p_k\) 进行一维线搜索确定步长 \(\alpha_k\)，更新迭代点 \(x^{(k+1)} = x^{(k)} + \alpha_k p_k\)，最后利用新旧迭代点的梯度差和位移差更新 \(B_k\) 到 \(B_{k+1}\)。
2. 问题初始化
- 目标函数：\(f(x) = (x_1 - 1)^2 + 2(x_2 - 2)^2 + 3(x_3 - 3)^2\)
- 梯度（解析形式）：
\[ \nabla f(x) = \begin{bmatrix} 2(x_1 - 1) \\ 4(x_2 - 2) \\ 6(x_3 - 3) \end{bmatrix} \]
- 初始点：\(x^{(0)} = (0, 0, 0)^T\)
- 初始近似Hessian矩阵：\(B_0 = I_3\)（3阶单位矩阵）
- 收敛阈值：\(\epsilon = 10^{-3}\)（当梯度范数 \(\| \nabla f(x^{(k)}) \| < \epsilon\) 时停止）
3. 第0次迭代（k=0）
(1) 计算梯度：
\[\nabla f(x^{(0)}) = \begin{bmatrix} 2(0-1) \\ 4(0-2) \\ 6(0-3) \end{bmatrix} = \begin{bmatrix} -2 \\ -8 \\ -18 \end{bmatrix}, \quad \| \nabla f(x^{(0)}) \| = \sqrt{(-2)^2 + (-8)^2 + (-18)^2} = \sqrt{4+64+324} = \sqrt{392} \approx 19.799 > \epsilon \]
不收敛，继续迭代。
(2) 计算搜索方向：解方程 \(B_0 p_0 = -\nabla f(x^{(0)})\)，由于 \(B_0 = I\)，直接得：
\[p_0 = -\nabla f(x^{(0)}) = \begin{bmatrix} 2 \\ 8 \\ 18 \end{bmatrix} \]
该方向实际上就是最速下降方向（因为 \(B_0\) 是单位阵）。
(3) 线搜索确定步长 \(\alpha_0\)：
线搜索目标：\(\min_{\alpha > 0} \phi(\alpha) = f(x^{(0)} + \alpha p_0)\)。
代入表达式：
\[x^{(0)} + \alpha p_0 = (2\alpha,\; 8\alpha,\; 18\alpha)^T \]
\[\phi(\alpha) = (2\alpha - 1)^2 + 2(8\alpha - 2)^2 + 3(18\alpha - 3)^2 \]
展开并合并同类项（也可通过求导解最优步长）：
\[\phi'(\alpha) = 2(2)(2\alpha-1) + 2\cdot 2(8)(8\alpha-2) + 3\cdot 2(18)(18\alpha-3) = 4(2\alpha-1) + 32(8\alpha-2) + 108(18\alpha-3) \]
整理系数：
\(4(2\alpha-1) = 8\alpha - 4\)
\(32(8\alpha-2) = 256\alpha - 64\)
\(108(18\alpha-3) = 1944\alpha - 324\)
合计：\((8+256+1944)\alpha + (-4-64-324) = 2208\alpha - 392\)
令 \(\phi'(\alpha) = 0\) 得：
\[\alpha_0 = \frac{392}{2208} = \frac{49}{276} \approx 0.177536 \]
这就是精确线搜索步长（因为目标函数是凸二次函数，线搜索有解析解）。
(4) 更新迭代点：
\[x^{(1)} = x^{(0)} + \alpha_0 p_0 = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix} + 0.177536 \cdot \begin{bmatrix} 2 \\ 8 \\ 18 \end{bmatrix} = \begin{bmatrix} 0.355072 \\ 1.420288 \\ 3.195648 \end{bmatrix} \]
(5) 计算新梯度：
\[\nabla f(x^{(1)}) = \begin{bmatrix} 2(0.355072-1) \\ 4(1.420288-2) \\ 6(3.195648-3) \end{bmatrix} = \begin{bmatrix} -1.289856 \\ -2.318848 \\ 1.173888 \end{bmatrix} \]
(6) 准备BFGS更新所需量：
位移 \(s_0 = x^{(1)} - x^{(0)} = \begin{bmatrix} 0.355072 \\ 1.420288 \\ 3.195648 \end{bmatrix}\)
梯度差 \(y_0 = \nabla f(x^{(1)}) - \nabla f(x^{(0)}) = \begin{bmatrix} -1.289856 - (-2) \\ -2.318848 - (-8) \\ 1.173888 - (-18) \end{bmatrix} = \begin{bmatrix} 0.710144 \\ 5.681152 \\ 19.173888 \end{bmatrix}\)
检查曲率条件：\(s_0^T y_0\) 必须为正，以保证 \(B_{k+1}\) 正定（对于凸问题通常成立）。计算：
\[s_0^T y_0 = 0.355072*0.710144 + 1.420288*5.681152 + 3.195648*19.173888 \approx 0.252 + 8.068 + 61.287 = 69.607 > 0 \]
满足条件。
(7) BFGS公式更新近似Hessian矩阵 \(B_k\)（这里更新的是Hessian近似，不是其逆；若更新逆矩阵 \(H_k = B_k^{-1}\) 则用对偶公式，步骤类似）：
BFGS更新公式为：
\[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} \]
计算各项（为节省篇幅，下面用近似数值，并展示计算逻辑）：
首先计算 \(B_0 s_0 = I s_0 = s_0\)（因为 \(B_0 = I\)）
分母1：\(s_0^T B_0 s_0 = s_0^T s_0 = 0.355072^2 + 1.420288^2 + 3.195648^2 \approx 0.126 + 2.017 + 10.212 = 12.355\)
分子1：外积矩阵 \(B_0 s_0 s_0^T B_0^T = s_0 s_0^T\) 是一个3×3矩阵，其第(i,j)元素为 \(s_{0,i} s_{0,j}\)。
第一项修正：\(-\frac{s_0 s_0^T}{12.355}\)（相当于每个元素除以12.355后从 \(B_0\) 中减去）。
分母2：\(s_0^T y_0 \approx 69.607\)（已算）
分子2：外积 \(y_0 y_0^T\) 是3×3矩阵。
第二项修正：\(+\frac{y_0 y_0^T}{69.607}\)。
因此：
\[B_1 = I - \frac{s_0 s_0^T}{12.355} + \frac{y_0 y_0^T}{69.607} \]
具体计算一个元素示例（比如(1,1)位置）：
\(s_0 s_0^T\) 的(1,1)元素 = \(0.355072^2 \approx 0.126\)
\(y_0 y_0^T\) 的(1,1)元素 = \(0.710144^2 \approx 0.504\)
所以：
\(B_1[1,1] = 1 - 0.126/12.355 + 0.504/69.607 \approx 1 - 0.0102 + 0.00724 = 0.99704\)
类似可算其他元素。实际上，由于目标函数是二次函数，BFGS更新应在有限步内达到精确Hessian矩阵 \(\nabla^2 f = \text{diag}(2, 4, 6)\)。经过一次更新，\(B_1\) 会逼近该对角阵。
4. 第1次迭代（k=1）
(1) 检查梯度范数：\(\| \nabla f(x^{(1)}) \| = \sqrt{(-1.289856)^2 + (-2.318848)^2 + (1.173888)^2} \approx \sqrt{1.664 + 5.377 + 1.378} = \sqrt{8.419} \approx 2.902 > \epsilon\)，继续。
(2) 解 \(B_1 p_1 = -\nabla f(x^{(1)})\) 得搜索方向 \(p_1\)。由于 \(B_1\) 已不是单位阵，需解3×3线性方程组。可验证 \(B_1\) 已接近对角阵，因此 \(p_1\) 的各分量近似为 \(-\nabla f(x^{(1)})\) 对应分量除以 \(B_1\) 对角元。
(3) 线搜索得步长 \(\alpha_1\)。
(4) 更新点 \(x^{(2)} = x^{(1)} + \alpha_1 p_1\)。
(5) 计算新梯度。
(6) 更新 \(B_2\) 用BFGS公式。
5. 后续迭代与收敛
对于这个凸二次函数，BFGS具有二次终止性质：从任意初始点出发，在n步内（n=3）达到精确解（在精确算术下）。实际上，由于本例是正定二次函数，且采用精确线搜索，BFGS会在至多3次迭代后使 \(B_k\) 收敛到精确Hessian，且 \(x^{(k)}\) 收敛到 \(x^*\)。
可以继续迭代计算，通常第2或第3次迭代后梯度范数就会小于 \(10^{-3}\)。
6. 关键点总结
- BFGS通过迭代更新 \(B_k\) 逐步逼近Hessian，避免了二阶导数计算。
- 线搜索需满足Wolfe条件以保证全局收敛；本例用精确线搜索演示。
- 每次迭代需存储和更新一个n×n矩阵（n为变量维数），适用于中小规模问题。
- 对于大规模问题，可采用有限内存L-BFGS，不显式存储矩阵。
通过以上步骤，你可清晰看到BFGS如何利用梯度信息和位移信息修正Hessian近似，从而产生比最速下降法更快的搜索方向。