Solution (source code)

= Solution

Use the extended target <probability density function>
$$
\rho(x,p)=\mu(x)(2\pi)^{-d/2}e^{-\|p\|^2/2}.
$$
The desired <stationary distribution> for positions is $\mu$; the full position-momentum <joint probability distribution> is $\rho$, not $\mu$ alone. Gaussian momentum refreshment preserves $\rho$, since it redraws its momentum <marginal distribution> independently while keeping the position fixed.

Let $R(x,p)=(x,-p)$ and $T=T_{\varepsilon,L}$ be the surrogate <leapfrog integration> map. The proposal in the paper is $S=R\circ T$: its notation $(X',-P')=T(x,p)$ includes the final momentum flip. Reversibility gives $RTR=T^{-1}$, so $S^2=I$. Both $R$ and $T$ preserve volume. Thus $S$ is an <involutive Metropolis proposal>, and the required acceptance <probability> is
$$
\boxed{\alpha(x,p,x',p')=1\wedge
\frac{\mu(x')e^{-\|p'\|^2/2}}{\mu(x)e^{-\|p\|^2/2}}.}
$$
Equivalently, in terms of the surrogate <Hamiltonian> $H_\nu$,
$$
\alpha=1\wedge\left\{
 e^{H_\nu(x,p)-H_\nu(x',p')}
 \frac{\mu(x')\nu(x)}{\mu(x)\nu(x')}
\right\}.
$$
Only evaluations of the true target <probability density function> are needed at endpoints; the trajectory uses the surrogate <gradient>. Target and surrogate <normalizing constants> cancel.

For $z=(x,p)$, the accepted flux satisfies
$$
\rho(z)\alpha(z)=\min\{\rho(z),\rho(Sz)\}=\rho(Sz)\alpha(Sz).
$$
Changing variables $z\mapsto Sz$ has unit absolute <Jacobian determinant>, so this identity proves <detailed balance> for accepted moves. The rejection mass stays at the same point and is also reversible. The accept/reject step therefore preserves $\rho$. Since momentum refreshment also preserves $\rho$, their composition and either phase of the alternating process preserve $\rho$. Projecting onto positions proves \b[the position <stationary distribution> is exactly the original target]. Stationarity alone does not establish uniqueness or convergence from every start; those require additional <irreducible Markov chain> and <aperiodic Markov chain> hypotheses.

There is a genuine defect in the printed smoothing claim. For the positive, nondifferentiable <Laplace distribution> $\mu(y)=e^{-|y|}/2$ on the real line, direct minimization gives, for every $\lambda>0$,
$$
\min_y\{\log\mu(y)+\lambda(x-y)^2\}
=-\log2-|x|-\frac1{4\lambda}.
$$
For $x>0$ the minimizer is $y=x+1/(2\lambda)$; for $x<0$ it is $y=x-1/(2\lambda)$; at zero both signs minimize. Hence normalization leaves $\nu=\mu$, still nondifferentiable at zero. The displayed construction does not in general produce the asserted smooth surrogate. For example, the also nondifferentiable target $\mu(y)\propto e^{-y^4-|y|}$ makes that infimum $-\infty$ for every $x$, since the negative quartic term dominates the quadratic penalty. The invariance proof above is valid conditional on actually having a usable smooth surrogate, as the subsequent algorithm assumes.

A corrected sufficient construction is <Moreau smoothing of a negative log-density>. For a <proper convex function> with <sequential lower semicontinuity> $U=-\log\mu$, set
$$
U_\lambda(x)=\min_y\{U(y)+\lambda\|x-y\|^2\},\qquad
\log\nu(x)=C-U_\lambda(x).
$$
The <Moreau envelope> is differentiable, with $\nabla U_\lambda(x)=2\lambda(x-\operatorname{prox}_{U/(2\lambda)}(x))$; require $e^{-U_\lambda}$ to be integrable to normalize it. The signs differ from the printed formula. For the <Laplace distribution>, this gives
$$
U_\lambda(x)-\log2=
\begin{cases}
\lambda x^2,&|x|\leq1/(2\lambda),\\
|x|-1/(4\lambda),&|x|>1/(2\lambda),
\end{cases}
$$
a genuinely differentiable potential with integrable exponential tails. Using that surrogate with the boxed acceptance <probability> still targets the original <Laplace distribution>.