From 977d885d1a53129a264b1aa1d20d65b06841f289 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sun, 31 May 2026 16:37:26 +0200 Subject: [PATCH] refactor(solver): H2 unify 5 Newton loops into newton_core (+B2/B3/B4/B5) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The five DCE solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean, Inversive-Distance) carried near-identical Newton loops differing only in the gradient function, the Hessian function, and the spherical concave sign (test-coverage audit H2 — "5 near-identical Newton loops"). Extract one detail::newton_core(x, grad, hess, concave, tol, max_iter); each public solver is now a thin wrapper passing two lambdas. This lets the performance fixes land once instead of five times: - B2 (api-perf): line_search now optionally returns the gradient at the accepted point; newton_core reuses it as the next iteration's convergence gradient, removing ~1 redundant full-mesh gradient evaluation per iteration. - B4 (api-perf): NewtonLinearSolver keeps one persistent SimplicialLDLT and runs analyzePattern once (symbolic factorization cached across iterations), re-analyzing only when nnz changes (self-healing for the FD Hessians that prune near-zero triplets). SparseQR fallback preserved verbatim. - B3 (api-perf): the Euclidean has_edge_dof scan is hoisted out of the loop (it is layout-invariant) — now done once before newton_core. - B5 (api-perf): the dead SimplicialLDLT variable in newton_euclidean is gone. Behaviour is unchanged: concave spherical still factors −H and solves (−H)·Δx = G; the merit steepest-descent still uses the true H; HardJava clamp default preserved. Tests (+2): NewtonCore.ReusesLineSearchGradient_B2_Convex asserts exactly 2 gradient evals on a unit-Hessian quadratic (vs 3 pre-B2), and ConcavePathConverges covers the −H branch directly. 298/298 CGAL tests pass, including all Java golden-vector parity tests. Co-Authored-By: Claude Opus 4.8 --- code/include/newton_solver.hpp | 455 +++++++++++-------------- code/tests/cgal/test_newton_solver.cpp | 75 ++++ 2 files changed, 280 insertions(+), 250 deletions(-) diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 3f7cbb4..12a3b9c 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -142,6 +142,14 @@ namespace detail { // re-open for the remaining helpers // If neither phase satisfies Armijo, return the best (smallest-residual) point // visited. If nothing beat ‖G(x)‖, return x unchanged and set *improved=false, // so the caller can stop cleanly instead of taking the old divergent full step. +// +// B2 (api-performance audit): the gradient computed at the accepted trial point +// is normally discarded (only its norm is kept), forcing the next Newton +// iteration to recompute G(x) from scratch. When `accepted_grad` is non-null it +// receives the gradient vector at the returned point, so the caller can reuse it +// as the next iteration's convergence gradient. It is only written when the +// returned point differs from `x` (i.e. `*improved == true`); on a pure stall +// (`*improved == false`) the caller must recompute G(x) itself. template inline std::vector line_search( const std::vector& x, @@ -149,9 +157,10 @@ inline std::vector line_search( const Eigen::VectorXd& d_sd, double norm0, GradFn&& grad_fn, - bool* improved = nullptr, - int max_halvings = 20, - double c1 = 1e-4) + bool* improved = nullptr, + std::vector* accepted_grad = nullptr, + int max_halvings = 20, + double c1 = 1e-4) { const int n = static_cast(x.size()); const double norm0_sq = norm0 * norm0; @@ -160,14 +169,18 @@ inline std::vector line_search( std::vector best_x = x; double best_norm = norm0; - // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew`. + std::vector trial_grad; // gradient at the most recent trial point + std::vector best_grad; // gradient at best_x (B2) + + // Evaluate ‖G(x + α·dir)‖₂, leaving the trial point in `xnew` and the + // gradient there in `trial_grad`. auto eval = [&](const Eigen::VectorXd& dir, double alpha) -> double { for (int i = 0; i < n; ++i) xnew[static_cast(i)] = x[static_cast(i)] + alpha * dir[i]; - auto Gnew = grad_fn(xnew); + trial_grad = grad_fn(xnew); double s = 0.0; - for (double v : Gnew) s += v * v; + for (double v : trial_grad) s += v * v; return std::sqrt(s); }; @@ -175,9 +188,10 @@ inline std::vector line_search( double alpha = 1.0; for (int ls = 0; ls < max_halvings; ++ls) { double norm_new = eval(dx, alpha); - if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; } if (norm_new * norm_new <= (1.0 - 2.0 * c1 * alpha) * norm0_sq) { if (improved) *improved = true; + if (accepted_grad) *accepted_grad = trial_grad; return xnew; } alpha *= 0.5; @@ -190,9 +204,10 @@ inline std::vector line_search( for (int ls = 0; ls < max_halvings; ++ls) { double thresh_sq = norm0_sq - 2.0 * c1 * alpha * dsd_sq; double norm_new = eval(d_sd, alpha); - if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; } + if (norm_new < best_norm) { best_norm = norm_new; best_x = xnew; best_grad = trial_grad; } if (thresh_sq >= 0.0 && norm_new * norm_new <= thresh_sq) { if (improved) *improved = true; + if (accepted_grad) *accepted_grad = trial_grad; return xnew; } alpha *= 0.5; @@ -201,10 +216,144 @@ inline std::vector line_search( // ── Both phases failed Armijo — never take the divergent full step. ─────── // Return the best point seen; if none improved, stay put and signal stall. - if (improved) *improved = (best_norm < norm0); + const bool any_improvement = (best_norm < norm0); + if (improved) *improved = any_improvement; + if (accepted_grad && any_improvement) *accepted_grad = best_grad; return best_x; } +// ── Cached linear solver (B4) ───────────────────────────────────────────────── +// +// solve_with_fallback rebuilds a fresh SimplicialLDLT every call, re-running the +// symbolic analysis (fill-reducing reordering / COLAMD) each Newton iteration — +// even though the Hessian sparsity pattern is iteration-invariant for a fixed +// mesh + DOF layout. This object keeps one persistent LDLT and runs +// `analyzePattern` once, then only `factorize` per iteration. +// +// FD-built Hessians (hyper-ideal, inversive-distance) prune near-zero triplets, +// so their pattern can shift slightly between iterations; we detect that via a +// change in `nonZeros()` and re-analyze, which keeps the cache correct and +// self-healing. The SparseQR fallback for singular/rank-deficient systems is +// preserved verbatim (built fresh on the rare path). The fallback semantics are +// identical to `solve_with_fallback`, so behaviour is unchanged. +struct NewtonLinearSolver { + Eigen::SimplicialLDLT> ldlt; + bool analyzed = false; + Eigen::Index last_nnz = -1; + bool fallback_used = false; ///< true iff the last solve used SparseQR + + Eigen::VectorXd solve(const Eigen::SparseMatrix& A, + const Eigen::VectorXd& rhs, + bool& ok) + { + fallback_used = false; + if (!analyzed || A.nonZeros() != last_nnz) { + ldlt.analyzePattern(A); + analyzed = true; + last_nnz = A.nonZeros(); + } + ldlt.factorize(A); + if (ldlt.info() == Eigen::Success) { + Eigen::VectorXd x = ldlt.solve(rhs); + if (ldlt.info() == Eigen::Success) { ok = true; return x; } + } + // Fallback: SparseQR — handles singular/rank-deficient A (gauge modes). + fallback_used = true; + Eigen::SparseQR, Eigen::COLAMDOrdering> qr(A); + if (qr.info() == Eigen::Success) { + Eigen::VectorXd x = qr.solve(rhs); + if (qr.info() == Eigen::Success) { ok = true; return x; } + } + ok = false; + return Eigen::VectorXd::Zero(rhs.size()); + } +}; + +// ── Unified Newton core (H2) ────────────────────────────────────────────────── +// +// All five DCE Newton solvers (Euclidean, Spherical, HyperIdeal, CP-Euclidean, +// Inversive-Distance) ran a near-identical loop that differed only in (a) the +// gradient function, (b) the Hessian function, and (c) the spherical sign quirk +// (the concave spherical energy factorises −H). This template captures that one +// loop so a fix lands once instead of five times. +// +// `grad(x) -> std::vector` : the energy gradient G(x). +// `hess(x) -> Eigen::SparseMatrix`: the TRUE Hessian H(x) (un-negated). +// `concave` : if true, solve (−H)·Δx = G (spherical); +// otherwise H·Δx = −G. Both are the same +// Newton system; `concave` only selects +// the PSD matrix actually factorised. +// +// Folds in B2 (reuse the line-search gradient as the next convergence gradient), +// B4 (cached symbolic factorization), and B5 (no dead solver variable). +template +inline NewtonResult newton_core( + std::vector x, + GradFn&& grad, + HessFn&& hess, + bool concave, + double tol, + int max_iter) +{ + const int n = static_cast(x.size()); + + NewtonResult res; + res.converged = false; + res.iterations = 0; + res.grad_inf_norm = 0.0; + + NewtonLinearSolver solver; + + // B2: gradient carried across the loop; the first one is the only "extra" + // evaluation — every later iteration reuses the line-search gradient. + std::vector G_std = grad(x); + + for (int iter = 0; iter < max_iter; ++iter) { + Eigen::Map G(G_std.data(), n); + + const double inf_norm = G.cwiseAbs().maxCoeff(); + if (inf_norm < tol) { + res.converged = true; + res.grad_inf_norm = inf_norm; + res.iterations = iter; + res.x = x; + return res; + } + + Eigen::SparseMatrix H = hess(x); + + // Newton system: factor the PSD matrix (H for convex, −H for concave). + Eigen::SparseMatrix A = concave ? Eigen::SparseMatrix(-H) : H; + Eigen::VectorXd rhs = concave ? Eigen::VectorXd(G) : Eigen::VectorXd(-G); + + bool ok = false; + Eigen::VectorXd dx = solver.solve(A, rhs, ok); + if (!ok) break; + + // Globalised line search (Armijo + steepest-descent fallback). + // d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖², + // using the TRUE (un-negated) Hessian for both signs. + const double norm0 = G.norm(); + Eigen::VectorXd d_sd = -(H * G); + bool improved = true; + std::vector G_next; + x = detail::line_search(x, dx, d_sd, norm0, grad, &improved, &G_next); + if (!improved) break; // line search stalled — stop cleanly + + G_std = std::move(G_next); // B2: reuse for the next convergence check + res.iterations = iter + 1; + } + + // Final gradient norm (reached only on max-iter / break — recomputed to be + // correct after a stalled or failed step). + auto G_final = grad(x); + double inf_final = 0.0; + for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); + res.grad_inf_norm = inf_final; + res.x = x; + return res; +} + } // namespace detail // ── Euclidean Newton solver ──────────────────────────────────────────────────── @@ -238,63 +387,21 @@ inline NewtonResult newton_euclidean( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); + // Layout is loop-invariant — decide the Hessian variant once (B3): cyclic + // layout (edge DOFs present) → analytic cyclic Hessian, which covers the + // vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian omits. + bool has_edge_dof = false; + for (auto e : mesh.edges()) + if (m.e_idx[e] >= 0) { has_edge_dof = true; break; } - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - Eigen::SimplicialLDLT> solver; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = euclidean_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ── - // Cyclic layout (edge DOFs present) → analytic cyclic Hessian, covering - // the vertex-edge / edge-edge blocks the vertex-only cotangent Laplacian - // omits. Vertex-only layout → analytic cotangent Laplacian. - bool has_edge_dof = false; - for (auto e : mesh.edges()) - if (m.e_idx[e] >= 0) { has_edge_dof = true; break; } - auto H = has_edge_dof ? euclidean_hessian_analytic(mesh, x, m) - : euclidean_hessian(mesh, x, m); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); // merit-function steepest descent - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return euclidean_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; // line search stalled — stop cleanly - - res.iterations = iter + 1; - } - - // Report final gradient norm - auto G_final = euclidean_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return euclidean_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { + return has_edge_dof ? euclidean_hessian_analytic(mesh, xc, m) + : euclidean_hessian(mesh, xc, m); + }, + /*concave=*/false, tol, max_iter); } // ── Spherical Newton solver ─────────────────────────────────────────────────── @@ -327,55 +434,14 @@ inline NewtonResult newton_spherical( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = spherical_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian: negate to get PSD; solve (−H)·Δx = G ────────────────── - auto H = spherical_hessian(mesh, x, m); - auto negH = Eigen::SparseMatrix(-H); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - // d_sd uses the actual (un-negated) Hessian: ∇f = H·G for f = ½‖G‖². - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return spherical_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = spherical_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Concave energy: the Hessian H is NSD, so newton_core factors −H (PSD) and + // solves (−H)·Δx = G — the same Newton system, just the sign that keeps LDLT + // well-defined. The merit steepest-descent inside still uses the true H. + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return spherical_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { return spherical_hessian(mesh, xc, m); }, + /*concave=*/true, tol, max_iter); } // ── HyperIdeal Newton solver ────────────────────────────────────────────────── @@ -419,53 +485,17 @@ inline NewtonResult newton_hyper_ideal( double hess_eps = 1e-5, HyperIdealScaleClamp clamp = HyperIdealScaleClamp::HardJava) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - // ── Gradient ────────────────────────────────────────────────────────── - auto G_std = evaluate_hyper_ideal(mesh, x, m, /*energy=*/false, /*grad=*/true, clamp).gradient; - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // ── Hessian (block-FD, ~33–1166× faster than full-FD) + solve H·Δx = −G - auto H = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps, clamp); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - // ── Globalised line search (Armijo + steepest-descent fallback) ─────── - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return evaluate_hyper_ideal(mesh, xnew, m, false, true, clamp).gradient; - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = evaluate_hyper_ideal(mesh, x, m, false, true, clamp).gradient; - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // block-FD Hessian (~33–1166× faster than full-FD); `clamp` selects the + // vertex-scale floor mode (N3). Convex energy → factor H directly. + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { + return evaluate_hyper_ideal(mesh, xc, m, /*energy=*/false, /*grad=*/true, clamp).gradient; + }, + [&](const std::vector& xc) { + return hyper_ideal_hessian_block_fd_sym(mesh, xc, m, hess_eps, clamp); + }, + /*concave=*/false, tol, max_iter); } // ── CP-Euclidean Newton solver (Phase 9a.1) ─────────────────────────────────── @@ -499,50 +529,12 @@ inline NewtonResult newton_cp_euclidean( double tol = 1e-8, int max_iter = 200) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - auto G_std = cp_euclidean_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - auto H = cp_euclidean_hessian(mesh, x, m); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return cp_euclidean_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = cp_euclidean_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Convex energy with an exact analytic Hessian (BPS-2015 2×2-per-edge). + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return cp_euclidean_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { return cp_euclidean_hessian(mesh, xc, m); }, + /*concave=*/false, tol, max_iter); } // ── Inversive-Distance Newton solver (Phase 9a.2) ───────────────────────────── @@ -581,51 +573,14 @@ inline NewtonResult newton_inversive_distance( int max_iter = 200, double hess_eps = 1e-5) { - std::vector x = x0; - const int n = static_cast(x.size()); - - NewtonResult res; - res.converged = false; - res.iterations = 0; - res.grad_inf_norm = 0.0; - - for (int iter = 0; iter < max_iter; ++iter) { - auto G_std = inversive_distance_gradient(mesh, x, m); - Eigen::Map G(G_std.data(), n); - - double inf_norm = G.cwiseAbs().maxCoeff(); - if (inf_norm < tol) { - res.converged = true; - res.grad_inf_norm = inf_norm; - res.iterations = iter; - res.x = x; - return res; - } - - // Hessian (block-FD, ~n/6× faster than full-FD) + solve H·Δx = −G. - auto H = inversive_distance_hessian_block_fd_sym(mesh, x, m, hess_eps); - bool ok = false; - Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok); - if (!ok) break; - - double norm0 = G.norm(); - Eigen::VectorXd d_sd = -(H * G); - bool improved = true; - x = detail::line_search(x, dx, d_sd, norm0, - [&](const std::vector& xnew) { - return inversive_distance_gradient(mesh, xnew, m); - }, &improved); - if (!improved) break; - - res.iterations = iter + 1; - } - - auto G_final = inversive_distance_gradient(mesh, x, m); - double inf_final = 0.0; - for (double v : G_final) inf_final = std::max(inf_final, std::abs(v)); - res.grad_inf_norm = inf_final; - res.x = x; - return res; + // Convex energy; block-FD Hessian (~n/6× faster than full-FD). + return detail::newton_core( + std::move(x0), + [&](const std::vector& xc) { return inversive_distance_gradient(mesh, xc, m); }, + [&](const std::vector& xc) { + return inversive_distance_hessian_block_fd_sym(mesh, xc, m, hess_eps); + }, + /*concave=*/false, tol, max_iter); } } // namespace conformallab diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp index 78b42bf..223fa22 100644 --- a/code/tests/cgal/test_newton_solver.cpp +++ b/code/tests/cgal/test_newton_solver.cpp @@ -579,3 +579,78 @@ TEST(SparseQRFallback, OkFlag_NullPointerIsSafe) EXPECT_NEAR(x[1], 5.0, 1e-12); }); } + +// ════════════════════════════════════════════════════════════════════════════ +// H2 / B2 — newton_core reuses the line-search gradient +// +// All five solvers now share detail::newton_core. B2: the gradient computed at +// the accepted line-search point is carried into the next iteration's +// convergence check instead of being recomputed at the loop top. On a +// unit-Hessian quadratic, exact Newton reaches the minimum in one step, so the +// total gradient-eval count is exactly 2 (1 initial + 1 accepted trial); the +// pre-B2 loop would have recomputed G at the top of iteration 1 → 3. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonCore, ReusesLineSearchGradient_B2_Convex) +{ + const std::vector xstar = {1.0, -2.0, 3.5}; + int grad_calls = 0; + + auto grad = [&](const std::vector& x) { + ++grad_calls; + std::vector g(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) g[i] = x[i] - xstar[i]; + return g; // G(x) = x − x*, H = +I + }; + auto hess = [&](const std::vector&) { + Eigen::SparseMatrix H(3, 3); + for (int i = 0; i < 3; ++i) H.insert(i, i) = 1.0; + H.makeCompressed(); + return H; + }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0, 0.0}, grad, hess, + /*concave=*/false, /*tol=*/1e-9, /*max_iter=*/50); + + ASSERT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_NEAR(res.x[0], 1.0, 1e-9); + EXPECT_NEAR(res.x[1], -2.0, 1e-9); + EXPECT_NEAR(res.x[2], 3.5, 1e-9); + EXPECT_EQ(grad_calls, 2) + << "B2: the loop-top gradient must be reused from the line search"; +} + +// Concave path (spherical-style): H is NSD, newton_core factors −H and solves +// (−H)·Δx = G. G(x) = x* − x with H = −I makes x* a maximiser; the core must +// still converge in one step. +TEST(NewtonCore, ConcavePathConverges) +{ + const std::vector xstar = {0.5, -1.5, 2.0}; + int grad_calls = 0; + + auto grad = [&](const std::vector& x) { + ++grad_calls; + std::vector g(x.size()); + for (std::size_t i = 0; i < x.size(); ++i) g[i] = xstar[i] - x[i]; + return g; // G(x) = x* − x, H = −I + }; + auto hess = [&](const std::vector&) { + Eigen::SparseMatrix H(3, 3); + for (int i = 0; i < 3; ++i) H.insert(i, i) = -1.0; + H.makeCompressed(); + return H; + }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0, 0.0}, grad, hess, + /*concave=*/true, /*tol=*/1e-9, /*max_iter=*/50); + + ASSERT_TRUE(res.converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_NEAR(res.x[0], 0.5, 1e-9); + EXPECT_NEAR(res.x[1], -1.5, 1e-9); + EXPECT_NEAR(res.x[2], 2.0, 1e-9); + EXPECT_EQ(grad_calls, 2); +}