feat(phase4): HyperIdeal Newton solver, SparseQR fallback, examples, docs
Phase 4 complete — 87 CGAL tests pass, 2 skipped. Newton solver (phase4a): - hyper_ideal_hessian.hpp: symmetric FD Hessian (O(ε²), PSD by convexity) - newton_hyper_ideal(): Newton + backtracking for the HyperIdeal functional - detail::solve_with_fallback(): optional bool* fallback_used parameter - solve_linear_system(): public API exposing LDLT→SparseQR fallback SparseQR fallback tests (SparseQRFallback.*): - FullRankSystem_CorrectSolution: LDLT path, fallback_used=false - SingularMatrix_FallbackActivated: zero-pivot → QR activated, fallback_used=true - Euclidean_ClosedMeshNoPinConverges: gauge-mode null space handled via QR HyperIdeal Newton tests (NewtonSolver.HyperIdeal_*): - ConvergesTriangleAllVariable, ResultFieldsConsistent, ConvergesTetrahedron, SparseQRFallbackNoCrash - Natural-target base point (b=1.0, a=0.5) — x=0 is degenerate in log-space Pipeline tests (test_pipeline.cpp): - End-to-end: all three geometries, mesh I/O round-trip, solve+export Example programs (code/examples/): - example_euclidean.cpp: headless Euclidean pipeline - example_hyper_ideal.cpp: headless HyperIdeal pipeline - example_viewer.cpp: interactive libigl viewer with jet colour map README: - Mathematical scope table: C++ vs Java original (18 rows) - "For mathematicians" section: mental model, step-by-step new-functional guide, half-edge traversal snippets, recommended reading Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>
This commit is contained in:
Tarik Moussa
87
code/include/hyper_ideal_hessian.hpp
Normal file
87
code/include/hyper_ideal_hessian.hpp
Normal file
@@ -0,0 +1,87 @@
|
||||
#pragma once
|
||||
// hyper_ideal_hessian.hpp
|
||||
//
|
||||
// Phase 4a — Hessian of the hyper-ideal discrete conformal functional.
|
||||
//
|
||||
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||||
// │ Implementation strategy │
|
||||
// │ │
|
||||
// │ The hyper-ideal functional involves angle functions (ζ, σ, α, β) │
|
||||
// │ composed through several nested layers (lij → ζ13/14/15 → β/α). │
|
||||
// │ Deriving closed-form Hessian entries analytically through all these │
|
||||
// │ layers is feasible but lengthy; an analytical Hessian is left for a │
|
||||
// │ future phase. │
|
||||
// │ │
|
||||
// │ Here we compute the Hessian by symmetric finite differences of the │
|
||||
// │ gradient, which is exact to O(ε²) and sufficient for Newton's method │
|
||||
// │ at the meshes typical in Phase 4 (< 500 DOFs): │
|
||||
// │ │
|
||||
// │ H[i,j] = (G(x + ε·eⱼ)[i] − G(x − ε·eⱼ)[i]) / (2ε) │
|
||||
// │ │
|
||||
// │ The hyper-ideal energy is strictly convex (Springborn 2020), so H is │
|
||||
// │ positive semi-definite everywhere and Eigen::SimplicialLDLT applies │
|
||||
// │ directly. │
|
||||
// └──────────────────────────────────────────────────────────────────────────┘
|
||||
|
||||
#include "hyper_ideal_functional.hpp"
|
||||
#include <Eigen/Sparse>
|
||||
#include <vector>
|
||||
#include <cmath>
|
||||
|
||||
namespace conformallab {
|
||||
|
||||
// ── Numerical Hessian via symmetric finite differences ────────────────────────
|
||||
//
|
||||
// Returns the n×n sparse Hessian, where n = hyper_ideal_dimension(mesh, m).
|
||||
// eps: finite-difference step size (default 1e-5 gives ~1e-10 relative error).
|
||||
inline Eigen::SparseMatrix<double> hyper_ideal_hessian(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const HyperIdealMaps& m,
|
||||
double eps = 1e-5)
|
||||
{
|
||||
const int n = hyper_ideal_dimension(mesh, m);
|
||||
std::vector<Eigen::Triplet<double>> trips;
|
||||
trips.reserve(static_cast<std::size_t>(n * n)); // dense upper bound
|
||||
|
||||
std::vector<double> xp = x, xm = x;
|
||||
|
||||
for (int j = 0; j < n; ++j) {
|
||||
const std::size_t sj = static_cast<std::size_t>(j);
|
||||
xp[sj] = x[sj] + eps;
|
||||
xm[sj] = x[sj] - eps;
|
||||
|
||||
auto Gp = evaluate_hyper_ideal(mesh, xp, m, /*energy=*/false).gradient;
|
||||
auto Gm = evaluate_hyper_ideal(mesh, xm, m, /*energy=*/false).gradient;
|
||||
|
||||
xp[sj] = xm[sj] = x[sj]; // restore
|
||||
|
||||
for (int i = 0; i < n; ++i) {
|
||||
double val = (Gp[static_cast<std::size_t>(i)]
|
||||
- Gm[static_cast<std::size_t>(i)]) / (2.0 * eps);
|
||||
if (std::abs(val) > 1e-15)
|
||||
trips.emplace_back(i, j, val);
|
||||
}
|
||||
}
|
||||
|
||||
Eigen::SparseMatrix<double> H(n, n);
|
||||
H.setFromTriplets(trips.begin(), trips.end());
|
||||
return H;
|
||||
}
|
||||
|
||||
// ── Symmetrised Hessian ───────────────────────────────────────────────────────
|
||||
//
|
||||
// The FD Hessian is symmetric in exact arithmetic; floating-point rounding
|
||||
// can introduce tiny asymmetries. This helper returns (H + Hᵀ)/2.
|
||||
inline Eigen::SparseMatrix<double> hyper_ideal_hessian_sym(
|
||||
ConformalMesh& mesh,
|
||||
const std::vector<double>& x,
|
||||
const HyperIdealMaps& m,
|
||||
double eps = 1e-5)
|
||||
{
|
||||
auto H = hyper_ideal_hessian(mesh, x, m, eps);
|
||||
Eigen::SparseMatrix<double> Ht = H.transpose();
|
||||
return (H + Ht) * 0.5;
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
@@ -1,27 +1,39 @@
|
||||
#pragma once
|
||||
// newton_solver.hpp
|
||||
//
|
||||
// Phase 4a — Newton solver for the discrete conformal functionals.
|
||||
// Phase 4a — Newton solver for all three discrete conformal functionals.
|
||||
//
|
||||
// Solves G(x) = 0 where G is the gradient of the discrete conformal energy:
|
||||
// G_v = Θ_v − Σ_f α_v^f (angle-sum residual at each vertex DOF)
|
||||
// Solves G(x) = 0 where G is the gradient of the discrete conformal energy.
|
||||
//
|
||||
// Algorithm per iteration:
|
||||
// 1. Compute gradient G(x)
|
||||
// 2. Check convergence: max|G_i| < tol → done
|
||||
// 3. Compute sparse Hessian H(x)
|
||||
// 4. Factorize and solve the Newton system:
|
||||
// Euclidean: H · Δx = −G (H is PSD → SimplicialLDLT directly)
|
||||
// Spherical: (−H) · Δx = G (H is NSD → negate to get PSD matrix)
|
||||
// 5. Backtracking line search: halve α until ||G(x+α·Δx)|| < ||G(x)||
|
||||
// 6. x ← x + α·Δx, go to 1
|
||||
// ┌──────────────────────────────────────────────────────────────────────────┐
|
||||
// │ Gradient sign conventions │
|
||||
// │ Euclidean / Spherical: G_v = Θ_v − Σ α_v (target − actual) │
|
||||
// │ HyperIdeal: G_v = Σ β_v − Θ_v (actual − target) │
|
||||
// │ │
|
||||
// │ All solvers use the same Newton step Δx = −H⁻¹·G │
|
||||
// │ │
|
||||
// │ Hessian sign at equilibrium │
|
||||
// │ Euclidean: H PSD → SimplicialLDLT on H │
|
||||
// │ Spherical: H NSD → SimplicialLDLT on −H (solve (−H)Δx = G) │
|
||||
// │ HyperIdeal: H PSD → SimplicialLDLT on H (analytical H: future) │
|
||||
// └──────────────────────────────────────────────────────────────────────────┘
|
||||
//
|
||||
// SparseQR fallback:
|
||||
// When SimplicialLDLT reports a failure (e.g. singular H on a closed mesh
|
||||
// without a pinned vertex), the solver automatically retries with
|
||||
// Eigen::SparseQR, which finds the minimum-norm Newton step orthogonal to
|
||||
// the null space. This handles the gauge mode on closed surfaces without
|
||||
// requiring the caller to pin a vertex explicitly.
|
||||
//
|
||||
// Requires:
|
||||
// Eigen::SimplicialLDLT (part of Eigen's sparse Cholesky module)
|
||||
// Eigen::SimplicialLDLT, Eigen::SparseQR (Eigen sparse module)
|
||||
|
||||
#include "euclidean_hessian.hpp"
|
||||
#include "spherical_hessian.hpp"
|
||||
#include "hyper_ideal_hessian.hpp"
|
||||
#include <Eigen/SparseCholesky>
|
||||
#include <Eigen/SparseQR>
|
||||
#include <Eigen/OrderingMethods>
|
||||
#include <Eigen/Dense>
|
||||
#include <algorithm>
|
||||
#include <cmath>
|
||||
@@ -41,6 +53,58 @@ struct NewtonResult {
|
||||
|
||||
namespace detail {
|
||||
|
||||
// Solve A·Δx = rhs with SimplicialLDLT; on failure fall back to SparseQR.
|
||||
// Returns Δx. ok is set to false only if both solvers fail.
|
||||
// If fallback_used is non-null, it is set to true iff SparseQR was needed.
|
||||
inline Eigen::VectorXd solve_with_fallback(
|
||||
const Eigen::SparseMatrix<double>& A,
|
||||
const Eigen::VectorXd& rhs,
|
||||
bool& ok,
|
||||
bool* fallback_used = nullptr)
|
||||
{
|
||||
if (fallback_used) *fallback_used = false;
|
||||
|
||||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> ldlt(A);
|
||||
if (ldlt.info() == Eigen::Success) {
|
||||
Eigen::VectorXd dx = ldlt.solve(rhs);
|
||||
if (ldlt.info() == Eigen::Success) { ok = true; return dx; }
|
||||
}
|
||||
// Fallback: SparseQR — handles singular/rank-deficient H (gauge modes).
|
||||
if (fallback_used) *fallback_used = true;
|
||||
Eigen::SparseQR<Eigen::SparseMatrix<double>, Eigen::COLAMDOrdering<int>> qr(A);
|
||||
if (qr.info() == Eigen::Success) {
|
||||
Eigen::VectorXd dx = qr.solve(rhs);
|
||||
if (qr.info() == Eigen::Success) { ok = true; return dx; }
|
||||
}
|
||||
ok = false;
|
||||
return Eigen::VectorXd::Zero(rhs.size());
|
||||
}
|
||||
|
||||
} // namespace detail
|
||||
|
||||
// ── Public linear-system solver (SparseQR fallback) ──────────────────────────
|
||||
//
|
||||
// Solve A·x = rhs with Eigen::SimplicialLDLT; if that fails (singular or
|
||||
// rank-deficient A), retry with Eigen::SparseQR which finds the minimum-norm
|
||||
// solution orthogonal to the null space.
|
||||
//
|
||||
// This is the same primitive used internally by all three Newton solvers.
|
||||
// Exposing it publicly lets callers (tests, downstream code) reuse the logic
|
||||
// and — via the optional fallback_used pointer — verify which code path ran.
|
||||
//
|
||||
// fallback_used – if non-null, set to true iff SparseQR was invoked
|
||||
// Returns Eigen::VectorXd::Zero(rhs.size()) if both solvers fail.
|
||||
inline Eigen::VectorXd solve_linear_system(
|
||||
const Eigen::SparseMatrix<double>& A,
|
||||
const Eigen::VectorXd& rhs,
|
||||
bool* fallback_used = nullptr)
|
||||
{
|
||||
bool ok = false;
|
||||
return detail::solve_with_fallback(A, rhs, ok, fallback_used);
|
||||
}
|
||||
|
||||
namespace detail { // re-open for the remaining helpers
|
||||
|
||||
// Backtracking line search: find the largest α in {1, 0.5, 0.25, …} such that
|
||||
// ||G(x + α·Δx)||₂ < ||G(x)||₂. Returns the accepted step (α may stay 1).
|
||||
template <typename GradFn>
|
||||
@@ -109,14 +173,11 @@ inline NewtonResult newton_euclidean(
|
||||
return res;
|
||||
}
|
||||
|
||||
// ── Hessian + factorisation ───────────────────────────────────────────
|
||||
// ── Hessian + solve H·Δx = −G (SparseQR fallback for singular H) ──
|
||||
auto H = euclidean_hessian(mesh, x, m);
|
||||
solver.compute(H);
|
||||
if (solver.info() != Eigen::Success) break;
|
||||
|
||||
// ── Newton step: solve H·Δx = −G ────────────────────────────────────
|
||||
Eigen::VectorXd dx = solver.solve(-G);
|
||||
if (solver.info() != Eigen::Success) break;
|
||||
bool ok = false;
|
||||
Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
|
||||
if (!ok) break;
|
||||
|
||||
// ── Backtracking line search ──────────────────────────────────────────
|
||||
double norm0 = G.norm();
|
||||
@@ -157,8 +218,6 @@ inline NewtonResult newton_spherical(
|
||||
res.iterations = 0;
|
||||
res.grad_inf_norm = 0.0;
|
||||
|
||||
Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>> solver;
|
||||
|
||||
for (int iter = 0; iter < max_iter; ++iter) {
|
||||
// ── Gradient ──────────────────────────────────────────────────────────
|
||||
auto G_std = spherical_gradient(mesh, x, m);
|
||||
@@ -173,15 +232,12 @@ inline NewtonResult newton_spherical(
|
||||
return res;
|
||||
}
|
||||
|
||||
// ── Hessian: negate to get PSD matrix ────────────────────────────────
|
||||
auto H = spherical_hessian(mesh, x, m);
|
||||
auto negH = Eigen::SparseMatrix<double>(-H);
|
||||
solver.compute(negH);
|
||||
if (solver.info() != Eigen::Success) break;
|
||||
|
||||
// ── Newton step: solve (−H)·Δx = G ─────────────────────────────────
|
||||
Eigen::VectorXd dx = solver.solve(G);
|
||||
if (solver.info() != Eigen::Success) break;
|
||||
// ── Hessian: negate to get PSD; solve (−H)·Δx = G ──────────────────
|
||||
auto H = spherical_hessian(mesh, x, m);
|
||||
auto negH = Eigen::SparseMatrix<double>(-H);
|
||||
bool ok = false;
|
||||
Eigen::VectorXd dx = detail::solve_with_fallback(negH, G, ok);
|
||||
if (!ok) break;
|
||||
|
||||
// ── Backtracking line search ──────────────────────────────────────────
|
||||
double norm0 = G.norm();
|
||||
@@ -201,4 +257,70 @@ inline NewtonResult newton_spherical(
|
||||
return res;
|
||||
}
|
||||
|
||||
// ── HyperIdeal Newton solver ──────────────────────────────────────────────────
|
||||
//
|
||||
// Solves G(x) = 0 for the hyper-ideal discrete conformal functional.
|
||||
//
|
||||
// Gradient sign convention (opposite to Euclidean/Spherical):
|
||||
// G_v = Σ β_v − Θ_v, G_e = Σ α_e − θ_e (actual − target)
|
||||
//
|
||||
// The hyper-ideal energy is strictly convex (Springborn 2020), so H is PSD
|
||||
// and SimplicialLDLT (with SparseQR fallback) applies directly.
|
||||
//
|
||||
// The Hessian is computed by symmetric finite differences of G (see
|
||||
// hyper_ideal_hessian.hpp); replace with an analytical Hessian in Phase 5.
|
||||
inline NewtonResult newton_hyper_ideal(
|
||||
ConformalMesh& mesh,
|
||||
std::vector<double> x0,
|
||||
const HyperIdealMaps& m,
|
||||
double tol = 1e-8,
|
||||
int max_iter = 200,
|
||||
double hess_eps = 1e-5)
|
||||
{
|
||||
std::vector<double> x = x0;
|
||||
const int n = static_cast<int>(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).gradient;
|
||||
Eigen::Map<const Eigen::VectorXd> 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 (numerical FD) + solve H·Δx = −G ─────────────────────────
|
||||
auto H = hyper_ideal_hessian_sym(mesh, x, m, hess_eps);
|
||||
bool ok = false;
|
||||
Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
|
||||
if (!ok) break;
|
||||
|
||||
// ── Backtracking line search ──────────────────────────────────────────
|
||||
double norm0 = G.norm();
|
||||
x = detail::line_search(x, dx, norm0,
|
||||
[&](const std::vector<double>& xnew) {
|
||||
return evaluate_hyper_ideal(mesh, xnew, m, false).gradient;
|
||||
});
|
||||
|
||||
res.iterations = iter + 1;
|
||||
}
|
||||
|
||||
auto G_final = evaluate_hyper_ideal(mesh, x, m, false).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;
|
||||
}
|
||||
|
||||
} // namespace conformallab
|
||||
|
||||
Reference in New Issue
Block a user