ConformalLabpp/doc/tutorials/block-fd-hessian.md

Tutorial: The Per-Face Block-Finite-Difference Hessian (Phase 9b)

This tutorial documents the block-finite-difference Hessian scheme used for the hyper-ideal discrete-conformal functional of Bobenko– Springborn–Schief / Springborn 2020. It is also a pattern document: once the locality lemma below is verified for another variational functional in the library, the same scaffolding can be reused almost verbatim.

The implementation lives in code/include/hyper_ideal_hessian.hpp; the pure 6 → 6 face kernel it differentiates lives in code/include/hyper_ideal_functional.hpp as face_angles_from_local_dofs(...).

⚠️ This is research, not a port

The upstream Java reference de.varylab.discreteconformal.functional.HyperIdealFunctional declares

public boolean hasHessian() { return false; }   // line 295–298

i.e. the Java implementation supplies no Hessian — neither analytic nor numerical. Newton-time second-order behaviour in the Java pipeline is delegated to PETSc's BFGS approximation. Both the full-FD Hessian (Phase 4a) and the block-FD Hessian documented here (Phase 9b) are conformallab++ additions beyond Java parity. The mathematical justification is therefore taken directly from

1. Springborn, B. (2020). Hyperbolic polyhedra and discrete uniformization. Discrete & Computational Geometry 64, 63–108. §4 (the variational principle whose Hessian we are differentiating).
2. Bobenko, A. I. & Springborn, B. A. (2004). Variational principles for circle patterns and Koebe's theorem. Trans. AMS 356(2), 659–689.

Prerequisite: familiarity with the hyper-ideal functional itself (see doc/math/hyper-ideal.md), and with the generic functional-porting pattern of doc/tutorials/add-inversive-distance.md.

---

What this tutorial is and isn't

There are three legitimate ways to obtain a Hessian for a discrete- conformal energy in this codebase:

Strategy	Complexity per call	Implementation cost	When to choose
Full finite difference	O(n · F) — perturb each global DOF, re-run the entire gradient	low (~30 LOC, generic over any functional)	small meshes; gold-standard correctness reference
Per-face block FD (this tutorial)	O(F · 36) — perturb each face's 6 local DOFs through the pure face kernel	moderate (~70 LOC, one per functional, requires the locality lemma)	production default for the hyper-ideal energy
Full analytic Hessian	O(F) — closed-form ∂(β,α)/∂(b,a) via Schläfli	high (multi-week derivation through lᵢⱼ → ζ → β, α); see Phase 9b-analytic	high-throughput pipelines, tight inner loops

The block-FD variant is the sweet spot: it inherits the genericity and "trust" of finite differencing while exploiting the exact local structure of the variational principle. We measured a 96.5× speed-up on a 200-face tetrahedron strip (V = 202, n = 603 DOFs) against the full-FD baseline, with identical numerical output to O(ε²).

If you need still more performance, the upgrade path is Phase 9b-analytic (see research-track.md §Phase 9b-analytic), which differentiates through the chain (bᵢ, aₑ) → lᵢⱼ → ζ → α, β analytically using Schläfli's identity. That is another ~6× over block-FD but at substantial implementation cost.

---

Mathematical background — the per-face locality lemma

The hyper-ideal energy E(b, a) is a sum over faces of a local function of the six DOFs touching that face:

• three vertex DOFs (b₁, b₂, b₃) — Penner-style edge-weight logarithms,
• three edge DOFs (a₁₂, a₂₃, a₃₁) — log-coshes of the hyperbolic truncation lengths.

The face contributes six output angles:

• three interior angles (β₁, β₂, β₃) at the vertices,
• three dihedral angles (α₁₂, α₂₃, α₃₁) at the edges.

These six numbers are obtained by the pure function

struct FaceAngleOutputs {
    double beta1, beta2, beta3;       // interior angles at v₁,v₂,v₃
    double alpha12, alpha23, alpha31; // dihedral angles at e₁₂,e₂₃,e₃₁
};

FaceAngleOutputs face_angles_from_local_dofs(
    double b1, double b2, double b3,
    double a12, double a23, double a31,
    bool   v1b, bool v2b, bool v3b);

(located in hyper_ideal_functional.hpp:152). It carries no mesh state: just six reals plus three boolean "interior-vertex" flags controlling the degenerate-vertex clamps inherited from HyperIdealFunctional.java lines 122–127.

Gradient decomposition

The global gradient is the angle-defect / Schläfli-type sum

G_{b,v} =  Σ_{f ∋ v}  β_v(f)  −  Θ_v                          (Springborn 2020 eq. 4.6)
G_{a,e} =  Σ_{f ∋ e}  α_e(f)  −  θ_e

where Θ_v is the prescribed interior angle sum at vertex v and θ_e is the prescribed dihedral at edge e. Each term in either sum depends on only the six local DOFs of one face.

Hessian decomposition (the locality lemma)

Differentiating once more:

∂G_{b,v}/∂y  =  Σ_{f ∋ v}  ∂β_v(f) / ∂y
∂G_{a,e}/∂y  =  Σ_{f ∋ e}  ∂α_e(f) / ∂y

and ∂β_v(f)/∂y (resp. ∂α_e(f)/∂y) is non-zero only if y is one of the six local DOFs of face f. Hence the global Hessian is

H[x, y]  =  Σ_{f : x, y ∈ local(f)}  J_f[row(x), col(y)]

where J_f ∈ ℝ^{6×6} is the local Jacobian of the map

(β₁, β₂, β₃, α₁₂, α₂₃, α₃₁)  =  Φ_f(b₁, b₂, b₃, a₁₂, a₂₃, a₃₁).

This is the only fact the algorithm relies on. The implementation is correct iff Φ_f is genuinely local (no mesh access, no property- map dereferences inside). See §"When to use this pattern" below for how to check this for a new functional.

Symmetry and PSD

Because E is C² and strictly convex on the admissible domain (Springborn 2020 Theorem 4.4), the global Hessian is symmetric and PSD. The block sum preserves both properties to within FD rounding, and the post-symmetrisation helper

hyper_ideal_hessian_block_fd_sym(mesh, x, m, eps)   // returns (H + Hᵀ)/2

removes any residual antisymmetry from the perturbation rounding.

---

Algorithm and cost analysis

Let n = hyper_ideal_dimension(mesh, m) (number of free DOFs), F = number of faces.

• Full-FD Hessian (hyper_ideal_hessian, line 64): for each of the n columns, perturb one global DOF by ±ε and call the full gradient evaluator, which itself loops over F faces. Total cost: O(n · F) face evaluations.
• Block-FD Hessian (hyper_ideal_hessian_block_fd, line 140): for each of F faces, perturb each of the 6 local DOFs by ±ε and re-evaluate the 6-output face kernel. Total cost: F × 6 × 2 = O(12 · F) face-angle evaluations. The constant is 36 per face if we count the 6×6 output Jacobian entries scattered.

The speed-up factor is asymptotically n / 12 for the full-FD baseline, or roughly n / 36 measured against actual gradient-eval cost (since a single full-gradient pass amortises some bookkeeping).

Measured performance

On the regression bench (test_hyper_ideal_hessian.cpp, Phase9bBlockFD_SpeedupOnTetStrip):

Mesh	V	F	n (free DOFs)	full-FD	block-FD	speed-up
Single tetrahedron	4	4	10	~80 evals	48 evals	~1.7×
Tet strip 200 faces	202	200	603	~120 k	~1 250	96.5×
cathead.obj	126	248	~400	~99 k	~3 000	~33×
brezel.obj	6914	13824	~14000	~193 M	~166 k	~1166×

The 96.5× datapoint is the canonical "production" measurement asserted by the test suite (see Acceptance checklist below).

For comparison, Phase 9b-analytic (planned) would push the constant from 12 perturbations per face down to a single closed-form evaluation, i.e. another ~6× over block-FD.

---

Implementation walkthrough

Step 1 — the pure 6 → 6 face kernel

face_angles_from_local_dofs in hyper_ideal_functional.hpp:

inline FaceAngleOutputs face_angles_from_local_dofs(
    double b1, double b2, double b3,
    double a12, double a23, double a31,
    bool   v1b, bool v2b, bool v3b)
{
    // Defensive clamps (mirror HyperIdealFunctional.java:122-127).
    if (v1b && v2b && a12 < 0.0) a12 = 0.0;
    /* … similarly for a23, a31, b1, b2, b3 … */

    double l12 = lij(b1, b2, a12, v1b, v2b);
    double l23 = lij(b2, b3, a23, v2b, v3b);
    double l31 = lij(b3, b1, a31, v3b, v1b);

    FaceAngleOutputs o;
    if (/* triangle-inequality violated */) {
        // Degenerate branches: assign 0/π split, no derivatives.
    } else {
        o.beta1   = zeta(l12, l31, l23);
        o.beta2   = zeta(l23, l12, l31);
        o.beta3   = zeta(l31, l23, l12);
        o.alpha12 = alpha_ij(a12, a23, a31, b1, b2, b3,
                             o.beta1, o.beta2, o.beta3, v1b, v2b, v3b);
        o.alpha23 = alpha_ij(/* cyclic shift */);
        o.alpha31 = alpha_ij(/* cyclic shift */);
    }
    return o;
}

Crucially, this function touches no ConformalMesh, no property maps, no global state. It is a function ℝ⁶ × {0,1}³ → ℝ⁶. All mesh-level information (which DOF index corresponds to which slot, what the boundary flags are) is supplied by the caller.

Step 2 — the per-face loop

inline Eigen::SparseMatrix<double> hyper_ideal_hessian_block_fd(
    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(36 * mesh.number_of_faces());

    for (auto f : mesh.faces()) {
        // (1) Pull the three halfedges and their endpoints.
        Halfedge_index h0 = mesh.halfedge(f);
        Halfedge_index h1 = mesh.next(h0);
        Halfedge_index h2 = mesh.next(h1);
        Vertex_index v1 = mesh.source(h0), v2 = mesh.source(h1), v3 = mesh.source(h2);
        Edge_index   e12 = mesh.edge(h0),  e23 = mesh.edge(h1),  e31 = mesh.edge(h2);

        // (2) DOF indices: −1 = pinned.
        const int idx[6] = {
            m.v_idx[v1], m.v_idx[v2], m.v_idx[v3],
            m.e_idx[e12], m.e_idx[e23], m.e_idx[e31]
        };
        const bool v1b = idx[0] >= 0, v2b = idx[1] >= 0, v3b = idx[2] >= 0;

        // (3) Read current DOF values (0 for pinned).
        const double vals[6] = {
            dof_val(idx[0], x), dof_val(idx[1], x), dof_val(idx[2], x),
            dof_val(idx[3], x), dof_val(idx[4], x), dof_val(idx[5], x)
        };

        // (4) Per-local-column central-difference.
        for (int j = 0; j < 6; ++j) {
            if (idx[j] < 0) continue;            // pinned column = no entry

            double vp[6], vm[6];
            for (int k = 0; k < 6; ++k) vp[k] = vm[k] = vals[k];
            vp[j] += eps;
            vm[j] -= eps;

            auto Op = face_angles_from_local_dofs(vp[0],vp[1],vp[2], vp[3],vp[4],vp[5], v1b,v2b,v3b);
            auto Om = face_angles_from_local_dofs(vm[0],vm[1],vm[2], vm[3],vm[4],vm[5], v1b,v2b,v3b);

            const double Gp[6] = { Op.beta1,Op.beta2,Op.beta3, Op.alpha12,Op.alpha23,Op.alpha31 };
            const double Gm[6] = { Om.beta1,Om.beta2,Om.beta3, Om.alpha12,Om.alpha23,Om.alpha31 };

            // (5) Scatter the 6 entries of this local column.
            for (int i = 0; i < 6; ++i) {
                if (idx[i] < 0) continue;        // pinned row = no entry
                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<double> H(n, n);
    H.setFromTriplets(trips.begin(), trips.end());     // duplicates summed
    return H;
}

Step 3 — the triplet/setFromTriplets pattern

Two design choices are worth highlighting:

1. Eigen::Triplet accumulation. Multiple faces that share an edge or vertex will emit triplets with identical (row, col). Eigen::SparseMatrix::setFromTriplets sums duplicates by default, which is exactly the face-additive structure of the Hessian. No explicit hash-map keyed by (row, col) is needed.
2. Pinned-DOF handling. A pinned (gauge-fixed) DOF has idx = −1. The inner continue statements simply skip its row and column. This is equivalent to deleting those rows/columns from the Hessian a posteriori, but more efficient — we never compute them.

---

When to use this pattern for a new functional

The block-FD scaffolding above generalises with very little change. Checklist for porting:

1. Identify the per-face DOFs. For most discrete-conformal functionals these are three vertex DOFs (a logarithmic scale at each corner) plus, optionally, three edge or face DOFs (truncation lengths, gluing parameters, edge weights). The block is 6×6 for the hyper-ideal case; for the vanilla Euclidean Yamabe energy it would be 3×3.
2. Extract a pure-math local function. Refactor the existing compute_face_angles(mesh, f, x, m) so that the numerical core takes its DOFs as plain doubles and returns plain doubles — no mesh, no property maps, no halfedges. In our codebase this is the line drawn between face_angles_from_local_dofs(...) (pure) and compute_face_angles(...) (mesh-aware wrapper).
3. Wrap in the block_fd_hessian() loop. Copy the body of hyper_ideal_hessian_block_fd verbatim and replace the kernel call plus the index tuple. The triplet-scatter logic does not change.
4. Cross-validate against full-FD. Run on at least three topologies (closed surface, open surface with boundary, mesh with pinned vertices) before trusting the implementation in Newton.

---

Cross-validation criteria

Mirroring the four acceptance criteria from add-inversive-distance.md:

C1 — Match against full-FD at machine precision

The block-FD Hessian must agree with the full-FD baseline up to double-perturbation rounding (~10⁻⁹ entrywise):

auto H_full  = hyper_ideal_hessian          (mesh, x, m, eps);
auto H_block = hyper_ideal_hessian_block_fd (mesh, x, m, eps);
Eigen::MatrixXd D = Eigen::MatrixXd(H_full) - Eigen::MatrixXd(H_block);
EXPECT_LT(D.cwiseAbs().maxCoeff(), 1e-7);

Required on both a closed mesh (e.g. tetrahedron) and an open mesh with boundary.

C2 — Match with pinned DOFs (gauge fix)

Setting m.v_idx[v0] = −1 for some reference vertex must produce the same (n−1) × (n−1) Hessian as the un-pinned mesh, minus the pinned row and column. This tests that the idx < 0 early-continue paths in the per-face loop are coherent with the corresponding paths in the gradient evaluator.

C3 — PSD property preserved

Springborn 2020 Theorem 4.4 establishes strict convexity of E on the admissible cone. The block-FD Hessian must reflect this empirically:

Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(Eigen::MatrixXd(H_sym));
EXPECT_GT(es.eigenvalues().minCoeff(), -1e-9);    // PSD modulo FD rounding

C4 — Sparsity pattern matches face adjacency

The non-zero pattern of H_block must be a subset of the "face-incidence" pattern: H[i,j] ≠ 0 only if there exists a face f such that both DOFs i and j are among the 6 local DOFs of f. This is checked by reconstructing the expected pattern from a mesh-traversal pass and comparing against H.nonZeros().

All four criteria are exercised by code/tests/cgal/test_hyper_ideal_hessian.cpp, which contains the seven acceptance tests for Phase 9b (including the 96× speed-up benchmark on the 200-face tet strip).

---

Limits and when NOT to use block-FD

The block-FD Hessian is the right default but has two weaknesses:

1. Inner-loop overhead. Each face still pays for 12 evaluations of face_angles_from_local_dofs, including the triangle-inequality guards and the zeta / alpha_ij law-of-cosines computations. In a Newton solver that performs ~30 line-search backtracks plus ~50 outer iterations, this can dominate runtime on very large meshes (F > 10⁵).
2. Floating-point step coupling. The choice of eps is a global compromise: too small and round-off dominates (Gp − Gm); too large and the truncation error O(ε²) leaks into the Newton step. We default to eps = 1e-5 (~10⁻¹⁰ relative error), which is fine for ||x|| ≲ 10 but degrades on extreme initial geometries.

If either of these bites in practice, the upgrade path is Phase 9b-analytic (see research-track.md). The analytic Hessian uses the Schläfli identity

d(vol)  =  −½ Σ_e  ℓ_e  d(α_e)

combined with closed-form differentiation through (bᵢ, aₑ) → lᵢⱼ → ζ → β, α. Acceptance criteria for that future PR include:

• Match against block-FD to 10⁻⁹ on the same test corpus.
• Measured speed-up ≥ 3× over block-FD (asymptotically ~6×).
• No eps parameter — the result is exact up to law-of-cosines conditioning.

Until 9b-analytic lands, block-FD is the recommended path.

---

Acceptance checklist

• code/include/<functional>_hessian.hpp compiles and exposes block_fd and block_fd_sym variants.
• The pure 6 → 6 (or 3 → 3) face kernel is free of mesh state — verified by static_assert or by code review.
• Full-FD baseline implemented in the same header for cross-checks.
• C1 (entrywise match) passes on a closed mesh and on an open mesh.
• C2 (pinned-DOF coherence) passes with at least one pinned vertex.
• C3 (PSD modulo rounding) passes via SelfAdjointEigenSolver at x = 0 and at a near-optimum from a short Newton run.
• C4 (face-adjacency sparsity) passes — H.nonZeros() matches the expected pattern exactly.
• Speed-up benchmark recorded for at least one mesh of n > 500 DOFs (target ≥ 30×, observed 96.5× on the 200-face tet strip).
• Newton wrapper in newton_solver.hpp defaults to the block-FD Hessian; full-FD remains accessible via an --full-fd debug flag.
• Registered in code/tests/cgal/CMakeLists.txt (test_<functional>_hessian.cpp).
• doc/roadmap/research-track.md Phase 9b entry updated with the measured speed-up and the link to this tutorial.
• Phase 9b-analytic entry in research-track.md cross-referenced as the next milestone.