更新：问题原来是我从分子和分母的向量中分解出一个小值然后计算点积/范数平方的步骤（参见下面的倒数第二段）。以直接的方式计算这些项导致 CG 步数与 Eigen Solver 的步数几乎相同。

我正在 C++ 中实现 Steihaug 方法，用于解决大规模无约束凸优化问题。当前实例，涉及大约 33000 个变量；理想情况下，问题大小将是 50 万到 100 万个变量。Hessian 矩阵是对称的且非常稀疏（约 30% 的非零项），并且具有关于主对角线的块结构。Hessian 是 PSD（从不严格 PD）并且有一个巨大的（~$10^8$) 条件编号。我正在使用对角线预处理器。科学问题是关于将基于 log-sum-exp 的平滑应用于组合优化问题的 LP 松弛。

随着该方法接近最优值，每次外部迭代需要越来越多的内部 CG 迭代：在最后一次外部迭代中达到约 7700 次 CG 迭代，达到全局最优。

当我将最后一次外部迭代的数据插入到 Eigen::BiCGSTAB 求解器中时，使用对角线预处理器，它会在 100 次迭代中收敛到所需的最优值。

我将不胜感激有关使实现在数字上稳健的一些指示。

编辑：我想一个重要的问题是梯度和粗麻布中的所有非零数字的数量级都很小。因为，我们正在接近全局最优值。

对于计算$\alpha$（最陡下降步）和$\beta$（方向更新系数），我正在分解最小的非零幅度（并在分子和分母之间取消它）并计算点积。这没有帮助。（更新：结果证明这是导致实施错误的不必要步骤。）

另外，我正在计算残差（如$Ax - b$) 每五十次迭代以补偿漂移。