where \(M(\theta )\)  is a user-chosen mass matrix, \(\lambda \)  is a vector of Legrange multipliers, and \(\pi (\theta )\)  is the unnormalized target density. Their constrained HMC algorithm draws proposal momentum values from \(\mathcal {N}\left(p_0 | 0, M(\theta _0)\right)\) , based on the current state of the Hamiltonian system \((p_0, \theta _0),\)  and then projects this momentum onto the constrained space such that (REF ) is satisfied. Movement through the system is then simulated using the RATTLE algorithm, which involves alternating between steps of the Leapfrog algorithm and updating steps via Newton's method to ensure the constraints continue to be satisfied [1]}. After simulating the Hamiltonian's movement, a Metropolis-Hastings acceptance step is taken, and the process is repeated. For further details of Brubaker's Constrained HMC, see the original paper [2]}.
