当前位置：首页 > article >正文

L-BFGS 方法实现

article 2025/1/31 6:19:27

L-BFGS 方法实现详细笔记

1. L-BFGS 算法简介

L-BFGS (Limited-memory Broyden-Fletcher-Goldfarb-Shanno) 是一种拟牛顿优化算法，适用于处理高维、无约束的优化问题，特别是可能非凸或非光滑的函数。

算法流程

初始化 $(x^0, g^0 \leftarrow \nabla f(x^0), B^0 \leftarrow I, k \leftarrow 0 )$ 。
进入循环：
- 计算搜索方向 $\leftarrow -B^k g^k )$ 。
- 使用 Lewis-Overton 线搜索找到步长 $(t)$ ，满足弱 Wolfe 条件。
- 更新变量 $(x^{k+1} \leftarrow x^k + t \cdot d )$ 。
- 更新梯度 $(g^{k+1} \leftarrow \nabla f(x^{k+1}) )$ 。
- 使用 Cautious-Limited-Memory-BFGS 更新近似 Hessian 信息 $B^{k+1} )$ 。
- 检查是否满足收敛条件或停止准则。
返回最优解 $x^* )$ 和对应的目标函数值 $f(x^*) )$ 。

2. 代码主要逻辑对照公式

1. 初始化阶段

公式：
$x^0, \, g^0 \leftarrow \nabla f(x^0), \, B^0 \leftarrow I, \, k \leftarrow 0$
对应代码：

Eigen::VectorXd xp(n);  // 当前点
Eigen::VectorXd g(n);   // 当前点的梯度
Eigen::VectorXd gp(n);  // 前一次的梯度
Eigen::VectorXd d(n);   // 搜索方向

// 计算初始目标函数值和梯度
fx = cd.proc_evaluate(cd.instance, x, g);
pf(0) = fx; // 存储目标函数值
d = -g;     // 初始搜索方向，负梯度方向

目标：初始化变量 $x^0 )$ 、梯度 $g^0 )$ 、目标函数值 $f(x^0) )$ ，并设置搜索方向为负梯度。
隐含假设：初始近似 Hessian 矩阵为单位矩阵 $B^0 = I )$ 。

2. 主循环条件

公式：
$\text{while} \, ||g^k|| > \delta$
对应代码：

gnorm_inf = g.cwiseAbs().maxCoeff();  // 梯度无穷范数
xnorm_inf = x.cwiseAbs().maxCoeff();  // 变量无穷范数

if (gnorm_inf / std::max(1.0, xnorm_inf) < param.g_epsilon) {
    ret = LBFGS_CONVERGENCE; // 收敛
    break;
}

目标：当梯度的无穷范数小于阈值 $(\delta )$ 时，认为优化已收敛。
代码：比较梯度范数与参数 param.g_epsilon，如果满足条件，终止循环。

3. 计算搜索方向

公式：
$\leftarrow -B^k g^k$
对应代码：

d = -g;  // 初始使用单位矩阵近似 Hessian

目标：使用当前梯度 $g^k )$ 和 Hessian 近似 $B^k )$ 计算搜索方向。
代码：初始迭代假设 $B^k = I )$ ，搜索方向为负梯度。

4. 线搜索

公式：
$\leftarrow \text{Lewis Overton line search}$
对应代码：

ls = line_search_lewisovertons(x, fx, g, step, d, xp, gp, step_min, step_max, cd, param);

目标：通过线搜索找到满足 Armijo 条件和弱 Wolfe 条件的步长 $(t)$ 。
代码：调用 line_search_lewisovertons 实现 Lewis-Overton 线搜索。

5. 更新变量和梯度

公式：
$x^{k+1} \leftarrow x^k + t \cdot d, \quad g^{k+1} \leftarrow \nabla f(x^{k+1})$
对应代码：

x = xp + step * d; // 更新变量
f = cd.proc_evaluate(cd.instance, x, g); // 更新目标函数值和梯度

目标：根据步长 $(t)$ 和方向 $(d)$ 更新变量，并重新计算目标函数值和梯度。

6. 更新 Hessian 信息

公式：
$B^{k+1} \leftarrow \text{Cautious-Limited-Memory-BFGS}(g^{k+1} - g^k, x^{k+1} - x^k)$
对应代码：

lm_s.col(end) = x - xp;    // s = x^{k+1} - x^k
lm_y.col(end) = g - gp;    // y = g^{k+1} - g^k

ys = lm_y.col(end).dot(lm_s.col(end)); // ys = y^T s
yy = lm_y.col(end).squaredNorm();      // yy = y^T y
lm_ys(end) = ys;

cau = lm_s.col(end).squaredNorm() * gp.norm() * param.cautious_factor;
if (ys > cau) {
    // 更新 BFGS 信息
    ...
}

目标：计算变量变化量 $(s)$ 和梯度变化量 $(y)$ ，并更新有限记忆存储。
代码：
- 通过 $ys = y^T s )$ 和 $yy = y^T y )$ 计算 Hessian 更新所需的比例因子。
- 判断 $\text{阈值} )$ ，执行谨慎更新。

7. 收敛和停止条件

公式：
$\text{若满足 } ||g^k|| \leq \delta \, \text{或其他停止准则}$
对应代码：

if (gnorm_inf / std::max(1.0, xnorm_inf) < param.g_epsilon) {
    ret = LBFGS_CONVERGENCE; // 梯度收敛
    break;
}

if (param.past > 0 && param.past <= k) {
    rate = std::fabs(pf(k % param.past) - fx) / std::max(1.0, std::fabs(fx));
    if (rate < param.delta) {
        ret = LBFGS_STOP; // 满足目标函数下降准则
        break;
    }
}

if (param.max_iterations != 0 && param.max_iterations <= k) {
    ret = LBFGSERR_MAXIMUMITERATION; // 达到最大迭代次数
    break;
}

目标：检测是否满足收敛条件（梯度、目标函数下降率、迭代次数）。
代码：
- 梯度收敛：梯度的无穷范数小于 param.g_epsilon。
- 目标函数下降率：变化量小于 param.delta。
- 最大迭代次数：超过 param.max_iterations。

3. 总结

L-BFGS 主要步骤

初始化：设置初始变量、梯度、目标函数值，假设初始 Hessian 为单位矩阵。
主循环：
- 计算搜索方向：通过负梯度和 Hessian 信息计算方向。
- 线搜索：使用 Lewis-Overton 方法找到满足条件的步长。
- 更新变量：通过步长更新变量和梯度。
- 更新 Hessian：根据变量和梯度变化更新近似信息。
检查停止准则：判断梯度、目标函数或迭代次数是否满足条件。
返回结果：输出最优解和目标函数值。

关键实现的亮点

线搜索：Lewis-Overton 方法保证适用于非光滑函数。
有限记忆：通过存储变量和梯度差向量节省内存。
谨慎更新：避免非正定的 Hessian 更新，提高稳定性。

L-BFGS 方法中 Hessian 更新

1. L-BFGS 方法中 Hessian 的作用

在优化问题中，Hessian 矩阵 $B^k )$ 是目标函数的二阶导数，用于描述函数的曲率信息。在 L-BFGS 方法中：

目标：通过近似 Hessian 矩阵 $B^k )$ ，找到更高效的搜索方向 $d = -B^k g )$ 。
特点：L-BFGS 不显式存储和计算完整的 Hessian 矩阵，而是通过有限记忆的递归公式更新搜索方向，从而降低内存和计算成本。

2. 公式回顾

更新公式

Hessian 更新公式：
$B^{k+1} \leftarrow \text{Cautious-Limited-Memory-BFGS}(g^{k+1} - g^k, x^{k+1} - x^k)$

$g^{k+1} - g^k )$ ：梯度变化。
$x^{k+1} - x^k )$ ：变量变化。
Cautious Update：通过添加阈值控制，确保更新的 Hessian 近似矩阵是正定的。

使用公式

在每次迭代中，使用更新后的 $B^{k+1} )$ 计算搜索方向：
$d = -B^{k+1} g^{k+1}$

3. 如何实现 Hessian 的更新

代码片段

/* 更新步长差向量和梯度差向量 */
lm_s.col(end) = x - xp;    // s = x^{k+1} - x^k
lm_y.col(end) = g - gp;    // y = g^{k+1} - g^k

/* 计算 ys 和 yy */
ys = lm_y.col(end).dot(lm_s.col(end)); // ys = y^T s
yy = lm_y.col(end).squaredNorm();      // yy = y^T y
lm_ys(end) = ys;

/* 如果 ys > 某一阈值，则进行更新 */
cau = lm_s.col(end).squaredNorm() * gp.norm() * param.cautious_factor;
if (ys > cau)
{
    ++bound; // 更新有限记忆的边界
    bound = m < bound ? m : bound; // 限制边界不超过存储容量 m
    end = (end + 1) % m; // 更新存储指针
}

实现解析

步长差向量 $s = x^{k+1} - x^k )$ ：
- 在每次迭代中，存储当前变量变化量 $(s)$ 。
- 对应代码：lm_s.col(end) = x - xp;。
梯度差向量 $y = g^{k+1} - g^k )$ ：
- 在每次迭代中，存储当前梯度变化量 $(y)$ 。
- 对应代码：lm_y.col(end) = g - gp;。
比例因子 $ys = y^T s )$ ：
- $(ys)$ 是 BFGS 更新的关键比例因子，必须满足 $(ys > 0)$ 才能确保 Hessian 近似矩阵正定。
- 对应代码：ys = lm_y.col(end).dot(lm_s.col(end));。
谨慎更新：
- 设置一个谨慎因子 $cautious_factor ) ( \text{阈值} = \|s\|^2 \|g^k\| \cdot \text{cautious\_factor} )$ ，避免 $(ys)$ 太小导致不稳定。
- 对应代码：
```
cau = lm_s.col(end).squaredNorm() * gp.norm() * param.cautious_factor;
if (ys > cau) {
    ...
}
```
有限记忆机制：
- 只保留最近 $(m)$ 次的 $(s)$ 和 $(y)$ 值，用于递归计算搜索方向。
- 对应代码：
```
++bound;
bound = m < bound ? m : bound;
end = (end + 1) % m;
```

4. 如何使用 $B^{k+1} )$

公式

在每次迭代中，根据有限记忆的 $(s)$ 和 $(y)$ 递归更新搜索方向 $(d)$ ：
$d = -B^{k+1} g$

代码片段

j = end;
for (i = 0; i < bound; ++i)
{
    j = (j + m - 1) % m; /* 逆序访问 */
    lm_alpha(j) = lm_s.col(j).dot(d) / lm_ys(j); // 计算 \alpha_j
    d += (-lm_alpha(j)) * lm_y.col(j); // 更新方向 d
}

d *= ys / yy; // 缩放方向 d

for (i = 0; i < bound; ++i)
{
    beta = lm_y.col(j).dot(d) / lm_ys(j); // 计算 \beta_j
    d += (lm_alpha(j) - beta) * lm_s.col(j); // 最终方向更新
    j = (j + 1) % m; /* 顺序访问 */
}

实现解析

逆序访问 $(s)$ 和 $(y)$ ：
- 遍历存储的历史值，从最近的一次开始，逐步递归调整搜索方向。
- 对应公式：
  $\alpha_j = \frac{s_j^T d}{y_j^T s_j}$
递归公式：
- 使用 $\alpha_j )$ 和 $\beta_j )$ 更新方向 $(d)$ ：
  $\sum_j (\alpha_j - \beta_j) s_j$
- 对应代码：
```
d += (-lm_alpha(j)) * lm_y.col(j);
d += (lm_alpha(j) - beta) * lm_s.col(j);
```
有限记忆约束：
- 只使用最近 $(m)$ 步的历史数据，避免存储完整的 Hessian 矩阵。