fix(mixed-mrf): port full-space Jacobian + Pfaffian split to Kyy block#104
Merged
Conversation
Completes the Fix 2 port from GGMGradientEngine to the mixed-MRF Kyy
block. Replaces the previous Roverato-style Cholesky Jacobian
log|det J| = q log 2 + sum_j (deg_higher_yy(j) + 2) psi_j
with the graph-agnostic split plus a per-column constraint correction:
ldj = q log 2 + sum_j (q + 1 - j) psi_j
pfaffian = 0.5 * sum_qq log det(A_qq diag(M_qq^{-1}) A_qq^T)
logp += ldj - pfaffian
For the full Kyy graph all A_qq are empty and the Pfaffian collapses to
zero, so the new and old formulas agree there. For sparse Kyy the
previous formula did not represent any RATTLE measure correction; the
new one is the correct manifold-marginal under any diagonal mass.
Applied to both gradient paths:
- logp_and_gradient (theta-space): identity mass; Pfaffian = 0 on full
Kyy, which is the only regime the dispatcher routes here.
- logp_and_gradient_full (full-space): mass-weighted Pfaffian using
this->inv_mass_ (empty => identity). Backward pass adds the matching
adjoint to R_bar via two triangular solves through chol(G_qq).
ensure_constraint_structure() is now called at the top of each gradient
function so chol_constraint_structure_.columns[col].excluded_indices is
available when building per-column A_qq. The projection helpers
(MixedMRFModel::project_position/_momentum) already built mass-weighted
A_qq matrices for RATTLE; the new gradient code consumes the same
constraint structure.
No regressions in `devtools::test()` (6989 pass, 0 fail).
Adds tests/testthat/test-mixed-gradient-pfaffian.R covering the new Cholesky-to-K + Pfaffian split in MixedMRFModel::logp_and_gradient_full. Four cases at max relative error < 1e-5: - full Kyy graph (refactor sanity; Pfaffian = 0) - sparse Kyy, 1 excluded edge, identity mass - sparse Kyy, 1 excluded edge, random diagonal mass - sparse Kyy, 2 excluded edges, random diagonal mass (multi-column) Plumbs an optional inv_mass_diag arg through mixed_test_logp_and_gradient_full so the integrator's mass-weighted Pfaffian can be exercised from R; the arg falls back to identity in the gradient engine when NULL. RcppExports regenerated via Rcpp::compileAttributes().
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #104 +/- ##
==========================================
- Coverage 86.51% 86.49% -0.03%
==========================================
Files 77 77
Lines 13269 13356 +87
==========================================
+ Hits 11480 11552 +72
- Misses 1789 1804 +15 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Completes the Fix 2 port from
GGMGradientEngineto the mixed-MRF Kyy block. Replaces the previous Roverato-style Cholesky Jacobianlog|det J| = q log 2 + sum_j (deg_higher_yy(j) + 2) psi_jwith the graph-agnostic split plus a per-column constraint correction:ldj = q log 2 + sum_j (q + 1 - j) psi_j
pfaffian = 0.5 * sum_qq log det(A_qq diag(M_qq^{-1}) A_qq^T)
logp += ldj - pfaffian
A_qqare empty and the Pfaffian collapses to zero, so the new and old formulas agree there.Applied to both gradient paths in
mixed_mrf_gradient.cpp. Theta-space uses identity mass (Pfaffian = 0 on full Kyy, which is the only regime the dispatcher routes here). Full-space usesthis->inv_mass_(empty ⇒ identity) and accumulates the matching adjoint intoR_barvia two triangular solves throughchol(G_qq).Mirrors the GGM Fix 2 pattern in
ggm_gradient.cppline-for-line.Test plan
devtools::test()passes (no regressions in the existing 6989 tests)mixed_test_logp_and_gradient_fullagrees with analytic gradient (max relative error < 1e-5 across 4 configurations: full Kyy, sparse Kyy with 1 or 2 excluded edges, identity and random diagonal mass) — new tests intest-mixed-gradient-pfaffian.Rmodel_type = "mixed"with sparse Kyy