Gradient-based MCMC requires computing the gradient of the unnormalized posterior density. In modern probabilistic programming frameworks, gradients are computed using automatic differentiation (AD) [1]}, [2]}. This involves building a differentiable computation graph for the target density, where all operations on inputs are recorded. However, computing the sensitivities efficiently and reliably in AD frameworks is not straightforward when iterative numerical solvers are involved [3]}, [2]}, [5]}.
