From b29eb4e0c4837880ccf62e7c6ead4f8fa5041e03 Mon Sep 17 00:00:00 2001 From: Tarik Moussa Date: Sun, 31 May 2026 16:03:12 +0200 Subject: [PATCH] =?UTF-8?q?perf(inv-dist):=20B1=20port=20=E2=80=94=20block?= =?UTF-8?q?-FD=20Hessian=20for=20newton=5Finversive=5Fdistance?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The Inversive-Distance solver built its Hessian inline by full finite differences: n perturbations × a full O(F) gradient eval = O(n·F) per Newton iteration (quadratic in mesh size), with no fast path at all (api-performance audit B1, second half). Port the per-face block-FD scheme already used by HyperIdeal (Phase 9b): the gradient decomposes by face (G_v = Θ_v − Σ_{f∋v} α_v, and each face's angles depend only on its 3 vertex DOFs), so the Hessian decomposes into per-face 3×3 blocks. Cost drops to O(F) face evaluations, a ≈ n/6 speed-up. - inversive_distance_functional.hpp: add the pure 3→3 kernel inversive_distance_face_grad_contribs (returns the per-face contribution −α to G; mirrors the gradient's face-skip on ℓ²≤0 exactly). - inversive_distance_hessian.hpp (new): full-FD baseline + block-FD + sym variants, mirroring hyper_ideal_hessian.hpp. - newton_solver.hpp: drop the inline full-FD lambda; call inversive_distance_hessian_block_fd_sym. - test: InversiveDistance_BlockFDHessianMatchesFullFD cross-validates the two Hessians entry-wise on a perturbed (off-equilibrium) config. 291/291 CGAL tests pass; all Inversive-Distance convergence tests unchanged. Co-Authored-By: Claude Opus 4.8 --- .../include/inversive_distance_functional.hpp | 37 ++++ code/include/inversive_distance_hessian.hpp | 187 ++++++++++++++++++ code/include/newton_solver.hpp | 44 +---- code/tests/cgal/test_newton_phase9a.cpp | 44 +++++ 4 files changed, 276 insertions(+), 36 deletions(-) create mode 100644 code/include/inversive_distance_hessian.hpp diff --git a/code/include/inversive_distance_functional.hpp b/code/include/inversive_distance_functional.hpp index 8d67a32..c973cf5 100644 --- a/code/include/inversive_distance_functional.hpp +++ b/code/include/inversive_distance_functional.hpp @@ -315,6 +315,43 @@ inline std::vector inversive_distance_gradient( return G; } +/// Per-face contribution to the Inversive-Distance gradient, as a pure +/// 3→3 kernel of the three local vertex DOFs `(u1,u2,u3)` and the three +/// edge inversive distances `(I12,I23,I31)`. No mesh, no property maps — +/// used by the per-face block-FD Hessian in `inversive_distance_hessian.hpp`. +/// +/// Returns the three values this face *adds* to the global gradient at +/// `(v1,v2,v3)`. Since `G_v = Θ_v − Σ_faces α_v`, the per-face contribution +/// is the **negative** corner angles `(−α₁,−α₂,−α₃)`. A non-real circle +/// configuration (any `ℓ² ≤ 0`) contributes nothing — exactly mirroring the +/// face-skip (`continue`) in `inversive_distance_gradient`, so the block-FD +/// Hessian and the full-FD Hessian see the same per-face support. +struct IDFaceGradContribs { + double g1; ///< contribution to G at v₁ (= −α₁) + double g2; ///< contribution to G at v₂ (= −α₂) + double g3; ///< contribution to G at v₃ (= −α₃) +}; + +inline IDFaceGradContribs inversive_distance_face_grad_contribs( + double u1, double u2, double u3, + double I12, double I23, double I31) +{ + double l12sq = id_detail::edge_length_squared(u1, u2, I12); + double l23sq = id_detail::edge_length_squared(u2, u3, I23); + double l31sq = id_detail::edge_length_squared(u3, u1, I31); + + // Same face-skip as inversive_distance_gradient: a non-real circle + // configuration has no limiting angle, so the face contributes 0. + if (l12sq <= 0.0 || l23sq <= 0.0 || l31sq <= 0.0) + return {0.0, 0.0, 0.0}; + + // euclidean_angles deliberately keeps its limiting (valid=false) angles + // here — the convex C¹ extension — matching the gradient's choice not to + // skip on !fa.valid. Corner at v_k = fa.alpha_k (see gradient Pass-2 trace). + auto fa = euclidean_angles(std::log(l12sq), std::log(l23sq), std::log(l31sq)); + return {-fa.alpha1, -fa.alpha2, -fa.alpha3}; +} + /// Inversive-Distance energy `E(u) = ∫₀¹ ⟨G(t·u), u⟩ dt`, evaluated /// with 10-point Gauss-Legendre (constants shared with `euclidean_energy`). inline double inversive_distance_energy( diff --git a/code/include/inversive_distance_hessian.hpp b/code/include/inversive_distance_hessian.hpp new file mode 100644 index 0000000..80a6d5e --- /dev/null +++ b/code/include/inversive_distance_hessian.hpp @@ -0,0 +1,187 @@ +#pragma once +// Copyright (c) 2024-2026 Tarik Moussa. +// SPDX-License-Identifier: MIT + +// inversive_distance_hessian.hpp +// +// Phase 9a.2 — Hessian of the inversive-distance circle-packing functional +// (Luo 2004 / Bowers-Stephenson 2004). +// +// ┌──────────────────────────────────────────────────────────────────────────┐ +// │ Implementation strategy │ +// │ │ +// │ TWO finite-difference Hessian implementations are provided here, │ +// │ mirroring the HyperIdeal pair in `hyper_ideal_hessian.hpp`: │ +// │ │ +// │ 1. `inversive_distance_hessian` — full finite-difference baseline. │ +// │ Cost ≈ n × (cost of a full gradient evaluation) = O(n · F). │ +// │ Kept as the cross-validation reference for the block-FD variant. │ +// │ │ +// │ 2. `inversive_distance_hessian_block_fd` — per-face block-FD. │ +// │ Each face contributes to the gradient through exactly 3 vertex │ +// │ DOFs (u₁,u₂,u₃); we FD the 3×3 local Jacobian of that face's │ +// │ gradient contribution and scatter it. Cost ≈ F × 6 face-angle │ +// │ evaluations = O(F). Speed-up factor ≈ n/6 over full-FD. │ +// │ │ +// │ Why the block-FD is correct (locality lemma): │ +// │ G_v = Θ_v − Σ_{f ∋ v} α_v(f), and α_v(f) depends ONLY on the 3 │ +// │ vertex DOFs of face f. Hence ∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)} │ +// │ ∂(−α_x)/∂y at f, so accumulating per-face 3×3 blocks reproduces the │ +// │ full Hessian (identical to O(ε²)). │ +// │ │ +// │ An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in │ +// │ `doc/roadmap/research-track.md` as Phase 9a.2-analytic; it would take │ +// │ the cost from O(F)·(FD constant) to a single O(F) analytic pass. │ +// └──────────────────────────────────────────────────────────────────────────┘ + +#include "inversive_distance_functional.hpp" +#include +#include +#include + +namespace conformallab { + +/// Full finite-difference Inversive-Distance Hessian (baseline). +/// Cost: `n` full-gradient evaluations ≈ `O(n·F)`. Use for small meshes +/// or as a correctness reference for the block-FD variant. +inline Eigen::SparseMatrix inversive_distance_hessian( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + const int n = inversive_distance_dimension(mesh, m); + std::vector> trips; + trips.reserve(static_cast(n) * 16); + + std::vector xp = x, xm = x; + for (int j = 0; j < n; ++j) { + const std::size_t sj = static_cast(j); + xp[sj] = x[sj] + eps; + xm[sj] = x[sj] - eps; + + auto Gp = inversive_distance_gradient(mesh, xp, m); + auto Gm = inversive_distance_gradient(mesh, xm, m); + + xp[sj] = xm[sj] = x[sj]; // restore + + for (int i = 0; i < n; ++i) { + double val = (Gp[static_cast(i)] + - Gm[static_cast(i)]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(i, j, val); + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +/// Symmetrised full-FD Inversive-Distance Hessian: `(H + Hᵀ)/2`. +inline Eigen::SparseMatrix inversive_distance_hessian_sym( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + auto H = inversive_distance_hessian(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + +// ── Block-FD Hessian ────────────────────────────────────────────────────────── +// +// The 3 local DOFs of a face f are (u_{v1}, u_{v2}, u_{v3}). For each free +// local DOF we recompute the face's gradient contribution (−α₁,−α₂,−α₃) at +// x ± ε along that axis and read off the 3×3 Jacobian. The result scatters +// into the global Hessian via the DOF-index lookup. +// +// Cost: F × 6 face-angle evaluations (3 DOFs × 2 directions) vs n×F for +// full-FD — a speed-up of ≈ n/6, i.e. ~hundreds× on large closed meshes. +/// Per-face block-FD Inversive-Distance Hessian. Uses the locality lemma +/// `∂G_x/∂y = Σ_{f: x,y ∈ {v1,v2,v3}(f)} ∂(−α_x)/∂y` to perturb only the 3 +/// face-local DOFs at a time, giving an `F·6` face-evaluation budget vs `n·F` +/// for full-FD. Mathematically equivalent to `inversive_distance_hessian` +/// up to O(ε²) FD rounding. +inline Eigen::SparseMatrix inversive_distance_hessian_block_fd( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + const int n = inversive_distance_dimension(mesh, m); + std::vector> trips; + trips.reserve(9 * mesh.number_of_faces()); + + for (auto f : mesh.faces()) { + Halfedge_index h0 = mesh.halfedge(f); + Halfedge_index h1 = mesh.next(h0); + Halfedge_index h2 = mesh.next(h1); + + Vertex_index v1 = mesh.source(h0); + Vertex_index v2 = mesh.source(h1); + Vertex_index v3 = mesh.source(h2); + Edge_index e12 = mesh.edge(h0); + Edge_index e23 = mesh.edge(h1); + Edge_index e31 = mesh.edge(h2); + + // Local DOF indices: (u1, u2, u3). Pinned slots = -1. + const int idx[3] = { m.v_idx[v1], m.v_idx[v2], m.v_idx[v3] }; + + const double I12 = m.I_e[e12]; + const double I23 = m.I_e[e23]; + const double I31 = m.I_e[e31]; + + // Local DOF values (0 for pinned). + const double vals[3] = { + id_detail::dof_val(idx[0], x), + id_detail::dof_val(idx[1], x), + id_detail::dof_val(idx[2], x) + }; + + for (int j = 0; j < 3; ++j) { + if (idx[j] < 0) continue; // never perturb a pinned DOF + + double vp[3], vm[3]; + for (int k = 0; k < 3; ++k) { vp[k] = vm[k] = vals[k]; } + vp[j] += eps; + vm[j] -= eps; + + auto Cp = inversive_distance_face_grad_contribs( + vp[0], vp[1], vp[2], I12, I23, I31); + auto Cm = inversive_distance_face_grad_contribs( + vm[0], vm[1], vm[2], I12, I23, I31); + + const double Gp[3] = { Cp.g1, Cp.g2, Cp.g3 }; + const double Gm[3] = { Cm.g1, Cm.g2, Cm.g3 }; + + for (int i = 0; i < 3; ++i) { + if (idx[i] < 0) continue; // pinned: contributes nothing + const double val = (Gp[i] - Gm[i]) / (2.0 * eps); + if (std::abs(val) > 1e-15) + trips.emplace_back(idx[i], idx[j], val); + } + } + } + + Eigen::SparseMatrix H(n, n); + H.setFromTriplets(trips.begin(), trips.end()); + return H; +} + +/// Symmetrised block-FD Inversive-Distance Hessian: `(H + Hᵀ)/2` of +/// `inversive_distance_hessian_block_fd(...)` for solvers requiring strict +/// symmetry. +inline Eigen::SparseMatrix inversive_distance_hessian_block_fd_sym( + const ConformalMesh& mesh, + const std::vector& x, + const InversiveDistanceMaps& m, + double eps = 1e-5) +{ + auto H = inversive_distance_hessian_block_fd(mesh, x, m, eps); + Eigen::SparseMatrix Ht = H.transpose(); + return (H + Ht) * 0.5; +} + +} // namespace conformallab diff --git a/code/include/newton_solver.hpp b/code/include/newton_solver.hpp index da2efe3..fd61248 100644 --- a/code/include/newton_solver.hpp +++ b/code/include/newton_solver.hpp @@ -36,6 +36,7 @@ #include "hyper_ideal_hessian.hpp" #include "cp_euclidean_functional.hpp" #include "inversive_distance_functional.hpp" +#include "inversive_distance_hessian.hpp" #include #include #include @@ -546,10 +547,11 @@ inline NewtonResult newton_cp_euclidean( /// domain where every triangle satisfies the inequalities. Luo's 1-form is /// closed there, so the path-integral energy is well-defined. /// -/// MVP implementation: the Hessian is computed by **finite differences** of -/// the analytic gradient (same pattern as the Phase 4a HyperIdeal solver). -/// An analytic Hessian via Glickenstein 2011 eq. (4.6) is tracked in -/// `doc/roadmap/research-track.md` as Phase 9a.2-analytic. +/// The Hessian is computed by **per-face block finite differences** +/// (`inversive_distance_hessian_block_fd_sym`, same pattern as the HyperIdeal +/// solver after Phase 9b) — O(F) face evaluations instead of the O(n·F) of the +/// full-FD baseline. An analytic Hessian via Glickenstein 2011 eq. (4.6) is +/// tracked in `doc/roadmap/research-track.md` as Phase 9a.2-analytic. /// /// \param mesh Input triangle mesh. /// \param x0 Initial DOF vector (length = number of free vertices). @@ -580,37 +582,6 @@ inline NewtonResult newton_inversive_distance( res.iterations = 0; res.grad_inf_norm = 0.0; - // Local FD Hessian builder — n × (cost of gradient eval). - auto build_hessian = [&](const std::vector& xc) -> Eigen::SparseMatrix { - std::vector> trips; - trips.reserve(static_cast(n) * 16); // sparse heuristic - - std::vector xp = xc, xm = xc; - for (int j = 0; j < n; ++j) { - const std::size_t sj = static_cast(j); - xp[sj] = xc[sj] + hess_eps; - xm[sj] = xc[sj] - hess_eps; - - auto Gp = inversive_distance_gradient(mesh, xp, m); - auto Gm = inversive_distance_gradient(mesh, xm, m); - - xp[sj] = xm[sj] = xc[sj]; // restore - - for (int i = 0; i < n; ++i) { - double val = (Gp[static_cast(i)] - - Gm[static_cast(i)]) - / (2.0 * hess_eps); - if (std::abs(val) > 1e-15) - trips.emplace_back(i, j, val); - } - } - Eigen::SparseMatrix H(n, n); - H.setFromTriplets(trips.begin(), trips.end()); - // Symmetrise — FD rounding may introduce tiny asymmetries. - Eigen::SparseMatrix Ht = H.transpose(); - return (H + Ht) * 0.5; - }; - for (int iter = 0; iter < max_iter; ++iter) { auto G_std = inversive_distance_gradient(mesh, x, m); Eigen::Map G(G_std.data(), n); @@ -624,7 +595,8 @@ inline NewtonResult newton_inversive_distance( return res; } - auto H = build_hessian(x); + // 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; diff --git a/code/tests/cgal/test_newton_phase9a.cpp b/code/tests/cgal/test_newton_phase9a.cpp index 03e5a1e..689e9d1 100644 --- a/code/tests/cgal/test_newton_phase9a.cpp +++ b/code/tests/cgal/test_newton_phase9a.cpp @@ -16,6 +16,7 @@ #include "newton_solver.hpp" #include "cp_euclidean_functional.hpp" #include "inversive_distance_functional.hpp" +#include "inversive_distance_hessian.hpp" #include "mesh_builder.hpp" #include "conformal_mesh.hpp" @@ -232,6 +233,49 @@ TEST(NewtonPhase9a, InversiveDistance_PerturbedTetrahedron_Converges) EXPECT_LT(res.grad_inf_norm, 1e-8); } +// ════════════════════════════════════════════════════════════════════════════ +// 6b. Inversive-Distance block-FD Hessian == full-FD Hessian (B1 port) +// +// The solver was switched from the O(n·F) full-FD Hessian to the O(F) +// per-face block-FD Hessian. The two must agree to FD rounding on any +// configuration — including a perturbed (non-equilibrium) one, where the +// off-diagonal coupling is non-trivial. This is the locality-lemma +// cross-validation that licenses the faster path. +// ════════════════════════════════════════════════════════════════════════════ + +TEST(NewtonPhase9a, InversiveDistance_BlockFDHessianMatchesFullFD) +{ + auto mesh = make_tetrahedron(); + auto m = setup_inversive_distance_maps(mesh); + compute_inversive_distance_init_from_mesh(mesh, m); + + // Pin one vertex; index the rest. + auto vit = mesh.vertices().begin(); + m.v_idx[*vit++] = -1; + int n = 0; + for (; vit != mesh.vertices().end(); ++vit) m.v_idx[*vit] = n++; + + // Evaluate the two Hessians away from equilibrium so off-diagonals are live. + std::vector x(static_cast(n), 0.0); + for (int i = 0; i < n; ++i) x[static_cast(i)] = 0.07 * (i + 1); + + auto H_full = inversive_distance_hessian_sym(mesh, x, m); + auto H_block = inversive_distance_hessian_block_fd_sym(mesh, x, m); + + ASSERT_EQ(H_full.rows(), H_block.rows()); + ASSERT_EQ(H_full.cols(), H_block.cols()); + + Eigen::MatrixXd Df = Eigen::MatrixXd(H_full); + Eigen::MatrixXd Db = Eigen::MatrixXd(H_block); + double max_abs_diff = (Df - Db).cwiseAbs().maxCoeff(); + EXPECT_LT(max_abs_diff, 1e-7) + << "block-FD and full-FD Hessians must agree to FD rounding;\n" + << "max |Δ| = " << max_abs_diff; + + // Sanity: the matrices are non-trivial (not both accidentally zero). + EXPECT_GT(Df.cwiseAbs().maxCoeff(), 1e-3); +} + // ════════════════════════════════════════════════════════════════════════════ // 7. CP-Euclidean Newton — uses analytic Hessian (NOT FD) //