multivariate_newton and multivariate_halley creep to max_iter near a singular Jacobian
メンテナーはふだん 1 日以内に返信
まだ誰も着手していません。
評価
調査の方向性
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分
- マージ済み PR（30日）
- 66
環境構築
このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。
はじめの一歩
- issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
- 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
- リポジトリをフォークし、ブランチを切って変更します。
- issue 番号を参照したプルリクエストを送ります。
boostorg/math のほかの issue
-
難易度 5/5 1週間以上 初心者へのやさしさ 35/100
メンテナーはふだん 1 日以内に返信
-
Reverse-mode autodiff: order-2 tape slows order-1 gradients, rewound tapes keep allocating blocks, rvar copies silently share a node対応中かも @NAThompson が 2 日前に担当しました。 オープン
難易度 4/5 3〜5日 初心者へのやさしさ 20/100
メンテナーはふだん 1 日以内に返信
-
Feature Wishlistオープンenhancement
難易度 5/5 1週間以上 初心者へのやさしさ 20/100
boostorg/math#303 · コメント 20 件 ·
メンテナーはふだん 1 日以内に返信
似ている issue
-
難易度 1/5 1時間未満 初心者へのやさしさ 65/100
-
[XLA] hlo/pass/README.md says passes implement Run; since the NVI change they implement RunImplオープン
難易度 1/5 1時間未満 初心者へのやさしさ 91/100
メンテナーはふだん 1 日以内に返信
-
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
NVIDIA-RTX/NRI#312 ·
メンテナーはふだん 1 日以内に返信
-
permute_q_b_proj / permute_wv_proj accept scale factors from another device対応中かも @0z5a が今日担当しました。 オープン
難易度 2/5 1〜3時間 初心者へのやさしさ 78/100
deepseek-ai/FlashMLA#247 ·
メンテナーはふだん 1 日以内に返信
-
Linux npm packaging leaves workspace files that block subsequent CI checkout対応中かも @Spockkk0225 が今日担当しました。 オープンci
難易度 2/5 1〜3時間 初心者へのやさしさ 73/100
alibaba/neug#1195 · コメント 1 件 ·
メンテナーはふだん 1 日以内に返信