diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index 4361e43..1be5041 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -48,12 +48,42 @@ namespace conformallab { // ── Result ──────────────────────────────────────────────────────────────────── -/// Result of `newton_solve(...)` — converged DOF vector + diagnostics. +/// Why a Newton solve terminated (I1 audit). Distinguishes the three +/// non-convergent exits that `converged == false` previously conflated. +enum class NewtonStatus { + Converged, ///< `grad_inf_norm < tol` — success. + MaxIterations, ///< hit `max_iter` without reaching `tol`. + LinearSolverFailed, ///< both LDLT and SparseQR failed on H·Δx = −G. + LineSearchStalled ///< the line search could not reduce the residual. +}; + +/// Human-readable name for a `NewtonStatus` (for logs / test messages). +inline const char* to_string(NewtonStatus s) noexcept +{ + switch (s) { + case NewtonStatus::Converged: return "Converged"; + case NewtonStatus::MaxIterations: return "MaxIterations"; + case NewtonStatus::LinearSolverFailed: return "LinearSolverFailed"; + case NewtonStatus::LineSearchStalled: return "LineSearchStalled"; + } + return "Unknown"; +} + +/// Result of a Newton solve — converged DOF vector + diagnostics. struct NewtonResult { std::vector x; ///< DOF vector at termination. - int iterations; ///< Newton steps taken. + int iterations; ///< Newton steps actually completed. double grad_inf_norm;///< max |Gᵢ| at termination. - bool converged; ///< `true` iff `grad_inf_norm < tol`. + bool converged; ///< `true` iff `status == Converged`. + + // I1: why it stopped (more informative than the `converged` bool alone). + NewtonStatus status = NewtonStatus::MaxIterations; + + // N7: linear-algebra conditioning diagnostics. + bool sparse_qr_fallback_used = false; ///< any iteration fell back to SparseQR. + double min_ldlt_pivot = 0.0; ///< smallest |Dᵢᵢ| of the last LDLT + ///< factorisation (≈ near-singularity + ///< proxy); 0 if LDLT never succeeded. }; // ── Internal helpers ────────────────────────────────────────────────────────── @@ -240,13 +270,16 @@ struct NewtonLinearSolver { Eigen::SimplicialLDLT> ldlt; bool analyzed = false; Eigen::Index last_nnz = -1; - bool fallback_used = false; ///< true iff the last solve used SparseQR + bool fallback_used = false; ///< true iff the last solve used SparseQR + double last_min_pivot = 0.0; ///< N7: smallest |Dᵢᵢ| of the last + ///< successful LDLT (0 if it fell back) Eigen::VectorXd solve(const Eigen::SparseMatrix& A, const Eigen::VectorXd& rhs, bool& ok) { - fallback_used = false; + fallback_used = false; + last_min_pivot = 0.0; if (!analyzed || A.nonZeros() != last_nnz) { ldlt.analyzePattern(A); analyzed = true; @@ -255,7 +288,14 @@ struct NewtonLinearSolver { ldlt.factorize(A); if (ldlt.info() == Eigen::Success) { Eigen::VectorXd x = ldlt.solve(rhs); - if (ldlt.info() == Eigen::Success) { ok = true; return x; } + if (ldlt.info() == Eigen::Success) { + ok = true; + // N7: smallest |Dᵢᵢ| in the LDLᵀ diagonal — a cheap + // near-singularity proxy (small ⇒ ill-conditioned step). + const auto D = ldlt.vectorD(); + last_min_pivot = (D.size() > 0) ? D.cwiseAbs().minCoeff() : 0.0; + return x; + } } // Fallback: SparseQR — handles singular/rank-deficient A (gauge modes). fallback_used = true; @@ -314,6 +354,7 @@ inline NewtonResult newton_core( const double inf_norm = G.cwiseAbs().maxCoeff(); if (inf_norm < tol) { res.converged = true; + res.status = NewtonStatus::Converged; // I1 res.grad_inf_norm = inf_norm; res.iterations = iter; res.x = x; @@ -328,7 +369,14 @@ inline NewtonResult newton_core( bool ok = false; Eigen::VectorXd dx = solver.solve(A, rhs, ok); - if (!ok) break; + // N7: accumulate conditioning diagnostics from this solve. + res.sparse_qr_fallback_used = res.sparse_qr_fallback_used || solver.fallback_used; + res.min_ldlt_pivot = solver.last_min_pivot; + if (!ok) { // I1/H1: linear solver failed this step + res.iterations = iter; // (iter steps actually completed) + res.status = NewtonStatus::LinearSolverFailed; + break; + } // Globalised line search (Armijo + steepest-descent fallback). // d_sd = −(H·G) is the merit-function steepest descent for f = ½‖G‖², @@ -338,12 +386,19 @@ inline NewtonResult newton_core( 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 + if (!improved) { // I1/H1: line search stalled — stop cleanly + res.iterations = iter; + res.status = NewtonStatus::LineSearchStalled; + break; + } G_std = std::move(G_next); // B2: reuse for the next convergence check res.iterations = iter + 1; } + // If we fell out of the loop without converging or breaking, it was max_iter + // (the default `res.status` is already MaxIterations; a break has set its own). + // Final gradient norm (reached only on max-iter / break — recomputed to be // correct after a stalled or failed step). auto G_final = grad(x); diff --git a/code/tests/cgal/test_newton_solver.cpp b/code/tests/cgal/test_newton_solver.cpp index 3649349..b60da71 100644 --- a/code/tests/cgal/test_newton_solver.cpp +++ b/code/tests/cgal/test_newton_solver.cpp @@ -500,6 +500,15 @@ TEST(SparseQRFallback, Euclidean_ClosedMeshNoPinConverges) << "Euclidean Newton on closed tetrahedron (no pin) must converge via SparseQR; " "grad_inf_norm = " << res.grad_inf_norm; EXPECT_LT(res.grad_inf_norm, 1e-8); + EXPECT_EQ(res.status, NewtonStatus::Converged); + // N7: the closed mesh with no pinned vertex has a gauge-singular Hessian. + // That near-singularity must be observable in the diagnostics — either the + // SparseQR fallback fired, or (when LDLT "succeeds" with a ~0 pivot, which is + // exactly the silent case N7 targets) the smallest LDLT pivot is tiny. + EXPECT_TRUE(res.sparse_qr_fallback_used || res.min_ldlt_pivot < 1e-6) + << "gauge singularity should surface via fallback flag or min_ldlt_pivot; " + << "fallback=" << res.sparse_qr_fallback_used + << " min_pivot=" << res.min_ldlt_pivot; } // ════════════════════════════════════════════════════════════════════════════ @@ -654,3 +663,77 @@ TEST(NewtonCore, ConcavePathConverges) EXPECT_NEAR(res.x[2], 2.0, 1e-9); EXPECT_EQ(grad_calls, 2); } + +// ════════════════════════════════════════════════════════════════════════════ +// I1 / H1 / N7 — newton_core termination status, iteration count, diagnostics +// ════════════════════════════════════════════════════════════════════════════ + +namespace { +// 1×1 sparse identity / scalar helpers for the synthetic newton_core tests. +inline Eigen::SparseMatrix diag_sparse(std::initializer_list d) +{ + const int n = static_cast(d.size()); + Eigen::SparseMatrix H(n, n); + int i = 0; + for (double v : d) { H.insert(i, i) = v; ++i; } + H.makeCompressed(); + return H; +} +} // namespace + +// I1: a converging solve reports status == Converged and the N7 pivot is the +// (well-conditioned) unit diagonal; no SparseQR fallback. +TEST(NewtonCore, Status_Converged_AndDiagnostics_N7) +{ + const std::vector xstar = {2.0, -1.0}; + auto grad = [&](const std::vector& x) { + return std::vector{ x[0] - xstar[0], x[1] - xstar[1] }; + }; + auto hess = [&](const std::vector&) { return diag_sparse({1.0, 1.0}); }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50); + + EXPECT_TRUE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::Converged); + EXPECT_EQ(res.iterations, 1); + EXPECT_FALSE(res.sparse_qr_fallback_used); // N7 + EXPECT_NEAR(res.min_ldlt_pivot, 1.0, 1e-12); // N7: D = I +} + +// I1: a slow (linearly-convergent cubic) problem that does not reach tol within +// max_iter reports MaxIterations, with iterations == max_iter (H1). +TEST(NewtonCore, Status_MaxIterations) +{ + // G(x) = x³, H = 3x² → Newton map x ← (2/3)x (linear convergence). + auto grad = [&](const std::vector& x) { + return std::vector{ x[0]*x[0]*x[0] }; + }; + auto hess = [&](const std::vector& x) { + return diag_sparse({ 3.0 * x[0] * x[0] }); + }; + + auto res = conformallab::detail::newton_core( + std::vector{1.0}, grad, hess, /*concave=*/false, /*tol=*/1e-8, /*max_iter=*/3); + + EXPECT_FALSE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::MaxIterations); + EXPECT_EQ(res.iterations, 3); +} + +// I1/H1: a constant non-zero gradient cannot be reduced by any step, so the +// line search stalls on the first iteration → LineSearchStalled, iterations == 0. +TEST(NewtonCore, Status_LineSearchStalled) +{ + auto grad = [&](const std::vector&) { + return std::vector{ 1.0, 1.0 }; // constant, independent of x + }; + auto hess = [&](const std::vector&) { return diag_sparse({1.0, 1.0}); }; + + auto res = conformallab::detail::newton_core( + std::vector{0.0, 0.0}, grad, hess, /*concave=*/false, 1e-9, 50); + + EXPECT_FALSE(res.converged); + EXPECT_EQ(res.status, conformallab::NewtonStatus::LineSearchStalled); + EXPECT_EQ(res.iterations, 0); // H1: no step completed +}