multivariate_newton and multivariate_halley creep to max_iter near a singular Jacobian

未关闭
  • 哪里不对？
#1,628 0 条评论 0 个 reaction 已指派 0 人 在 GitHub 查看

维护者通常 1 天内回复

还没有人认领这个 Issue。

评估

难度
5/5
预计耗时
一周以上
新手友好度
12/100
Issue 类型
缺陷
描述清晰度
基本清楚
活跃度
活跃
技术栈
cpp
领域
backend

调研方向

The solvers live in include/boost/math/tools/multivariate_roots.hpp, with multivariate_newton and multivariate_halley sharing the Newton step, backtracking line search and convergence tests. Start from the reproducer, which runs Freudenstein–Roth from (15, -2) and reports 1000 iterations and about 6000 function calls. Done would mean the solver stops early near the singular Jacobian with a clear non-root result, and the existing tests still pass. The proposed fix is a design change to the step and line search, so it needs maintainer guidance first.

由索引模型根据 Issue 内容生成。

描述

On Freudenstein–Roth, multivariate_newton and multivariate_halley can zigzag along a line where J is nearly singular for all of max_iter, returning a point that is neither a root nor a local minimum of ‖F‖:

F1 = -13 + x1 + ((5 - x2) x2 - 2) x2,   F2 = -29 + x1 + ((x2 + 1) x2 - 14) x2

J depends only on x2, and det J = 6x2² − 8x2 − 12 vanishes on the whole line x2 ≈ −0.8968.

method, start iterations calls of f result ‖F‖
Halley, (15, −2) 1000 (limit) 6098 (9.15095, −0.896394) 7.698
Newton, (8.424166, −0.933377) 1000 (limit) 6249 (9.13988, −0.897481) 7.696
Newton, (9.15, −0.8964) 1000 (limit) 7001 (9.16012, −0.896196) 7.695
Newton, (15, −2) throws evaluation_error 96

At the end, φ = ½‖F‖² falls by 10⁻⁵ to 10⁻⁶ (relative) per iteration. The points are not stationary: the scaled gradient max|gᵢ|·max(|xᵢ|, 1)/max(φ, 1), with g = JᵀF, is about 2. But the cosine between the Newton step and −g is only 6×10⁻⁵ to 9×10⁻⁵. The backtracking line search always finds a tiny decrease, so neither the singular-J test nor the line-search failure test fires. The non-root minimum is at about (11.4128, −0.8968), with ‖F‖ = 6.9989.

This doesn't contradict the docs, which allow stopping at max_iter. But it costs thousands of evaluations to return a point with no meaning.

Possible fix:

  • Compute g = JᵀF before the LU overwrites J.
  • When −gᵀa/(‖g‖‖a‖) is below a threshold, backtrack along a scaled −g instead, using the actual slope gᵀp in the Armijo test and the quadratic interpolation. The current −2φ is right only for the Newton step.
  • Report a small g with a nonzero F as a non-root stationary point.

A dogleg step (as in MINPACK's hybrj) is the larger alternative. It would also help on the Rosenbrock gradient system, where the docs note that Newton takes 253 iterations.

Reproducer 
#include <boost/math/tools/multivariate_roots.hpp>
#include <array>
#include <cmath>
#include <cstdio>
#include <tuple>
using V = std::array<double, 2>;
static long calls = 0;
auto fr = [](V const& x) {
   ++calls;
   double y = x[1];
   V F{-13 + x[0] + ((5 - y) * y - 2) * y, -29 + x[0] + ((y + 1) * y - 14) * y};
   std::array<V, 2> J{V{1, 10 * y - 3 * y * y - 2}, V{1, 3 * y * y + 2 * y - 14}};
   auto d2 = [y](V const& v) { return V{(10 - 6 * y) * v[1] * v[1], (6 * y + 2) * v[1] * v[1]}; };
   return std::make_tuple(F, J, d2);
};
int main()
{
   std::uintmax_t it = 1000;
   V x = boost::math::tools::multivariate_halley(fr, V{15, -2}, 53, it);
   auto F = std::get<0>(fr(x));
   std::printf("%ju iterations, %ld calls, x = (%.6g, %.6g), |F| = %.4g\n", it, calls - 1, x[0], x[1], std::hypot(F[0], F[1]));
   // 1000 iterations, 6098 calls, x = (9.15095, -0.896394), |F| = 7.698
}

Found with the WebAssembly demo's basins-of-attraction plots.

主要语言
C++
星标
378
派生
271
平均合并
12 小时 37 分钟
30 天内合并 PR
66

环境准备

这个项目没有提供开发容器、Dockerfile 或贡献指南，环境需要你自己搭建：先看它的 README，通用步骤见我们的新手贡献指南。

从这里开始

  1. 先读完整个 Issue，再读项目的贡献指南。
  2. 在 Issue 下留言说明你要接手 —— 这能避免两个人做同样的事。
  3. Fork 仓库，在一个分支上完成修改。
  4. 提交 Pull Request，并在描述里引用这个 Issue 编号。

boostorg/math 的其他 Issue

查看 boostorg/math 的全部 Issue

相似的 Issue

更多 C++ Issue

把新 issue 发到你的邮箱

精选适合新手参与的 GitHub issue 摘要。