conformallab/ConformalLabpp
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 1 Wiki Activity

feat(solver): S2 — NewtonResult status enum (I1) + iteration fix (H1) + conditioning diagnostics (N7)

Browse Source 
Builds on the newton_core refactor; all three findings live in exactly that code.

I1 (test-coverage): add NewtonStatus { Converged, MaxIterations,
LinearSolverFailed, LineSearchStalled } and a `status` field to NewtonResult, so
the three non-convergent exits that `converged == false` previously conflated are
now distinguishable.  `converged` is kept (== status==Converged) for back-compat;
+ to_string(NewtonStatus) for logs/tests.  newton_core sets the status at each
exit point (centralised by the H2 refactor).

H1 (test-coverage): set res.iterations explicitly at the LinearSolverFailed and
LineSearchStalled breaks (= completed steps), instead of relying on the last
successful iteration's stale value.

N7 (numerical-stability): surface linear-algebra conditioning in the result —
`sparse_qr_fallback_used` (any iteration fell back to SparseQR) and
`min_ldlt_pivot` (smallest |Dᵢᵢ| of the last LDLT, a cheap near-singularity
proxy).  This catches the silent case the audit flagged: on a gauge-singular
Hessian SimplicialLDLT "succeeds" with a ~0 pivot and no fallback fires — now
min_ldlt_pivot is tiny and observable.

Tests (+3 synthetic newton_core status/diagnostic tests; +1 assertion on the
closed-mesh-no-pin fallback test).  All purely additive — no control-flow or
convergence behaviour changes; 301/301 CGAL tests pass incl. Java parity.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
2026-05-31 18:39:19 +02:00
parent 74d41e99a8
commit 1bf1defc70
2 changed files with 146 additions and 8 deletions
Download Patch File Download Diff File Expand all files Collapse all files

67
code/include/newton_solver.hpp
View File

@@ -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<double> 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 ──────────────────────────────────────────────────────────
@@ -241,12 +271,15 @@ struct NewtonLinearSolver {
bool analyzed = false;
Eigen::Index last_nnz = -1;
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<double>& A,
const Eigen::VectorXd& rhs,
bool& ok)
{
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<double> 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);

83
code/tests/cgal/test_newton_solver.cpp
View File

@@ -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<double> diag_sparse(std::initializer_list<double> d)
{
const int n = static_cast<int>(d.size());
Eigen::SparseMatrix<double> 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<double> xstar = {2.0, -1.0};
auto grad = [&](const std::vector<double>& x) {
return std::vector<double>{ x[0] - xstar[0], x[1] - xstar[1] };
};
auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
auto res = conformallab::detail::newton_core(
std::vector<double>{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<double>& x) {
return std::vector<double>{ x[0]*x[0]*x[0] };
};
auto hess = [&](const std::vector<double>& x) {
return diag_sparse({ 3.0 * x[0] * x[0] });
};
auto res = conformallab::detail::newton_core(
std::vector<double>{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<double>&) {
return std::vector<double>{ 1.0, 1.0 }; // constant, independent of x
};
auto hess = [&](const std::vector<double>&) { return diag_sparse({1.0, 1.0}); };
auto res = conformallab::detail::newton_core(
std::vector<double>{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
}