text:
\begin{document}

\title{Hessian-Free High-Resolution Nesterov Acceleration for Sampling}

\author{\name Ruilin Li \email ruilin.li@gatech.edu \\
       \addr School of Mathematics\\
       School of Computational Science and Engineering\\
       Georgia Institute of Technology
       \AND
       \name Hongyuan Zha \email zhahy@cuhk.edu.cn \\
       \addr School of Data Science \\
       The Chinese University of Hong Kong, Shenzhen \\
       Shenzhen Institute of Artificial Intelligence and Robotics for Society
       \AND
       \name Molei Tao \email mtao@gatech.edu \\
       \addr School of Mathematics\\
       Georgia Institute of Technology
}

\editor{

}

\maketitle

\begin{abstract}
We propose an accelerated-gradient-based MCMC method. It relies on a modification of the Nesterov's accelerated gradient method for strongly convex functions (NAG-SC): We reformulate NAG-SC as a Hessian-Free High-Resolution ODE, lift its high-resolution coefficient to be a hyperparameter,

and then inject appropriate noise and discretize the resulting diffusion process. Accelerated sampling enabled by the new hyperparameter is theoretically quantified. At continuous-time level, for log-concave/log-strongly-concave target measures,  exponential convergence in $\chi^2$ divergence/2-Wasserstein distance is proved, with a rate analogous to the state-of-the-art results of underdamped Langevin dynamics, plus an \textbf{additional acceleration}. At discrete algorithm level, a dedicated discretization algorithm is proposed to simulate the Hessian-Free High-Resolution SDE in a cost-efficient manner. For log-strong-concave target measures, the proposed algorithm achieves $\tilde{\mathcal{O}}(\frac{\sqrt{d}}{\epsilon})$ iteration complexity in 2-Wasserstein distance, same as underdamped Langevin dynamics, but with a reduced constant.
Empirical experiments are conducted to numerically verify our theoretical results and demonstrate accelerated sampling in practice.

\end{abstract}

\begin{keywords}
Markov Chain Monte Carlo, Nesterov's accelerated gradient, accelerated sampling, Langevin dynamics, optimization in the space of probabilities
\end{keywords}

\section{Introduction}
Optimization methods have been a major algorithmic machinery that drives both the theory and practice of machine learning in recent years. Since the seminal work of Nesterov \citep{nesterov1983method}, acceleration has played a key role in gradient-based optimization methods. One notable example is the Nesterov's Accelerated Gradient (NAG) method, which is an instance of a more general family of ``momentum methods". NAG in fact consists of multiple methods, including NAG-C for convex functions, and NAG-SC for strongly convex functions, both of which have provably faster convergences than the vanilla gradient descent (GD) method in their corresponding setups \citep{nesterov1983method, nesterov2013introductory}. Although they are classical methods, significant new perspectives of acceleration have recently been studied, e.g., \cite{su2014differential, wibisono2016variational, wilson2016lyapunov, hu2017dissipativity, attouch2018fast,shi2018understanding}. This work aims at turning NAG-SC into a sampler, and `NAG' from hereon will refer to NAG-SC unless confusion arises.

In fact, approaches for sampling statistical distributions, such as gradient-based Markov Chain Monte Carlo (MCMC) methods, also remain of great importance in machine learning, for example due to their links to statistical inference and abilities to represent uncertainties lacking in optimization-based methods. Although they are not entirely the same thing, there are profound and interesting interplays between optimization and sampling. For example, the perspective of viewing sampling as optimization in probability space dates back to at least the 90s \citep{jordan1998variational}, and is gaining increasing attention in the machine learning commmunity (e.g., \cite{liu2016stein, dalalyan2017further, wibisono2018sampling, zhang2018policy, frogner2018approximate, chizat2018global, chen2018unified, ma2019there, erdogdu2020convergence}). In fact, an unadjusted Euler-Maruyama discretization of overdamped Langevin dynamics (abbreviated as OLD in this article) is commonly considered as the analog of GD in sampling (although many other discretizations are also possible), and oftentimes referred to as Unadjusted Langevin Algorithm (ULA) \citep{roberts1996exponential} and/or Langevin Monte Carlo (LMC). The convergence properties of the continuous dynamics of OLD, as well as asymptotic and non-asymptotic analyses of its discretizations have been continuously studied (e.g., \cite{roberts1996exponential, villani2008optimal, pavliotis2014stochastic,  dalalyan2017theoretical, durmus2016sampling, dalalyan2017further, durmus2019analysis, durmus2019high, vempala2019rapid, cheng2018convergence, dwivedi2019log, ma2019sampling, chewi2020optimal, erdogdu2020convergence}).

Meanwhile, the notion of acceleration is less quantified in sampling compared to that in optimization, although attention has been rapidly building up. Along this direction, one line of research is based on diffusion processes such as underdamped Langevin dynamics (ULD), which is closely related to OLD. For example, the convergence and nonasymptotics of discretized ULD have been studied in \cite{cheng2017underdamped, dalalyan2018sampling, ma2019there}, and were demonstrated provably faster than discretized OLD in suitable setups. These are not only great progresses but also forming perspectives complementary to the extensive studies of the convergence of continuous ULD in the mathematical community (e.g, \cite{ mattingly2002ergodicity, cao2019explicit, dolbeault2009hypocoercivity, dolbeault2015hypocoercivity, villani2009hypocoercivity, eckmann2003spectral, baudoin2013bakry, eberle2019couplings}). Another important line of research is related to accelerating particle-based approaches for optimization in probability spaces \citep{liu2019understanding, pmlr-v97-taghvaei19a, wang2019accelerated}, although we note there is no clear boundary between these two lines (e.g., \cite{leimkuhler2018ensemble}). More generally, it has been known that adding an irreversible part to the reversible dynamics of OLD\footnote{For irreversibility-induced-acceleration \emph{not} based on OLD, see e.g., \cite{bierkens2019zig} and references therein.} accelerates its convergence (e.g., \cite{hwang2005accelerating,lelievre2013optimal,ohzeki2015langevin,rey2015irreversible,duncan2016variance}), and this work can be viewed to be under this umbrella, although the discretization of an accelerated continuous process is also important, and will be analyzed.

More precisely, we proposed an accelerated gradient-based MCMC algorithm termed HFHR, that is based on diffusion process and inspired by a simple yet natural motivation: how to appropriately inject noise to NAG algorithm \textbf{in discrete time}, so that it is turned into an algorithm for momentum-accelerated sampling? Note we don't want to add noise to the learning-rate$\to 0$ limit of NAG, which has been well studied \citep{ma2019there}, because the discretization of this low-resolution ODE by a finite step size may not converge as fast as NAG with the same learning rate. However, we will still use continuous dynamics as intermediate steps, and our roadmap is the following: first, view NAG as the discretization of a high-resolution ordinary differential equation (ODE), then convert it to a stochastic differential equation (SDE) by injecting noise appropriately, and finally discretize the SDE appropriately to obtain a fast and efficient HFHR algorithm.

More precisely, the first step combines ingredients from the existing literature so that  a non-asymptotic formulation can be prepared for the later steps. The goal is to better account for NAG's behavior when a finite (not infinitesimal) learning rate is used. As pointed out in \cite{shi2018understanding}, a low-resolution limiting ODE \citep{su2014differential}, albeit being a milestone leading to a new venue of research (e.g, \cite{wibisono2016variational}), does not fully capture the acceleration enabled by NAG --- for example, it can't distinguish between NAG and other momentum methods such as heavy ball \citep{polyak1964some}. The main reason is, the low-resolution ODE describes the $h\to 0$ limit of NAG, but in practice NAG uses a finite (nonzero) $h$. High-resolution ODE was thus proposed to include additional $\mathcal{O}(h)$ terms to account for the finite $h$ effect \citep{shi2018understanding}, which led to a better characterization of NAG. The original form of the high-resolution ODE involves the Hessian of the objective function, which is computationally expensive to evaluate and store for high-dimensional problems, but this difficulty can be bypassed via a change of variable (e.g., \citealp{alvarez2002second, attouch2020first}), which allows us to derive a High-Resolution and Hessian-Free limiting ODE for NAG.

Then we make the following algorithmic innovations. One is to replace a specific coefficient in the HFHR ODE by a hyperparameter $\alpha \ge 0$, which, as we will demonstrate both theoretically and empirically, can lead to accelerated convergence of the eventual sampling algorithm. The other is to add noise to the hyperparametrized HFHR ODE in a specific way, which turns it into an SDE suitable for the sampling purpose. This SDE will be termed as HFHR dynamics. In order to obtain an actual algorithm, the HFHR SDE needs to be discretized, and we'll propose a relatively simple discretization and analyze it to  demonstrate that the accelerated convergence is \textbf{not} an artificial consequence of time-rescaling, which would not give acceleration after discretization with an appropriate step size. Meanwhile, we'd like to point out that our discretization is just one of the many possible schemes. It was known that high-order discretizations can improve statistical accuracy and even the speed of convergence (see e.g., \cite{chen2015convergence,li2019stochastic};

the analogue in (stochastic) optimization has also been studied, e.g., \cite{zhang2018direct}), although such improvements often come with a cost of more computations per iteration. The discretization considered here will just be a simple first-order scheme that uses one (batch-)gradient evaluation per step, but it better utilizes the structure of HFHR dynamics than Euler-Maruyama.

Our presentation will be structured as follows. After detailing the construction of HFHR, we will theoretically analyze the convergence of HFHR, at both the continuous level (HFHR dynamics) and the discrete level (HFHR algorithm). We will first show that HFHR dynamics admits the target distribution as its invariant distribution (Theorem \ref{thm:invariant_distribution}) and converges exponentially fast to it when the target distribution satisfies a Poincar\'e's inequality (Theorem \ref{thm:simple_result}). Then we move on to a more specific setup of log-concave / log-strongly-concave target distributions, which are commonly considered in the literature \citep{kim2016global, bubeck2018sampling, dalalyan2017theoretical, dalalyan2018sampling, dwivedi2019log, shen2019randomized},
and demonstrate explicitly an additional acceleration of HFHR when compared to ULD (Theorems \ref{thm:main} and \ref{thm:exp_convergence_coupling}). For our discretized HFHR algorithm, a non-asymptotic error bound will be obtained (Theorem \ref{thm:discretization}), which confirms that the acceleration of HFHR over ULD in continuous time carries through to the discrete territory, at least for log-strongly-concave target distributions. Finally, the theoretical analysis is complemented by numerical experiments, demonstrating the performance of HFHR for a series of representative target distributions as well as Bayesian neural network learning tasks. The validity and tightness of our theoretical results will also be numerically verified.

The main contribution of this article is the idea of adding noise to the NAG-SC optimization algorithm to turn it into a sampler, which provides a new perspective that is neither overdamped or underdamped Langevin. Theoretical analyses and numerical experiments are provided just to validate this new perspective. For instance, the discretized algorithm we provide (note it's not unique as different discretizations can be used) achieves an $\polylog(\frac{\sqrt{d}}{\epsilon})$ iteration complexity in 2-Wasserstein distance, which has the same order as the best-known results for KLMC; meanwhile, we show theoretically that the scaling constant in this bound can be reduced by HFHR (Theorem \ref{thm:wasserstein} and Corollary \ref{corollary:iteration complexity}), suggesting acceleration, and accelerated convergence was empirically and consistently observed as well (Section \ref{sec:experiment}).

\section{Background: Langevin Dynamics for Sampling}

Consider the sampling from Gibbs measure $\mu$ whose density is
\begin{equation}\label{eq:target_distribution}
    d\mu = \frac{1}{\int e^{- f(\bs{y})}d\bs{y}} e^{-f(\bs{x})} d\bs{x} \quad \mbox{ w.r.t. Lebesgue measure in } \mathbb{R}^d
\end{equation}
where $f: \mathbb{R}^d \mapsto \mathbb{R}$ will be called the potential function.

Many gradient-based MCMC methods have been proposed for this sampling problem. A large proportion of them involve two steps: first, choose a Markov process whose stationary distribution is the target distribution $\mu$, and then discretize the process and simulate for long enough until the generated samples follow an empirical distribution sufficiently close to the target distribution. We will focus on methods based on Langevin dynamics.

\paragraph{Continuous Process.}
At the continuous level, two diffusion processes popular for sampling (and modeling important physical processes too) are both named after Langevin. One is called overdamped Langevin dynamics (OLD), and the other is called kinetic Langevin dynamics (abbreviated as ULD to comply with a convention of calling it underdamped Langevin). They are respectively given by

\begin{align}
    &\mbox{(OLD)} \qquad d\qbold_t = -\nabla f(\qbold_t)dt + \sqrt{2}d\bs{W}_t \label{eq:LD} \\
    &\mbox{(ULD)} \quad
    \begin{cases}
        d\qbold_t = \pbold_t dt\\
        d\pbold_t = -\gamma \pbold_t dt - \nabla f(\qbold_t) dt + \sqrt{2\gamma} d\bs{B}_t
    \end{cases}
\end{align}

where $\qbold_t, \pbold_t \in \mathbb{R}^d$, $\bs{W}_t, \bs{B}_t$ are i.i.d. Wiener processes in $\mathbb{R}^d$, and $\gamma > 0$ is a friction coefficient. Under mild conditions (e.g., \cite{pavliotis2014stochastic}), OLD converges to $\mu$ and ULD converges to
\begin{equation}\label{eq:joint_invariant_distribution}
    d\pi(\qbold, \pbold) = d\mu(\qbold) \nu(\pbold) d\pbold, \, \text{ where } \nu(\pbold) = (2\pi)^{-\frac{d}{2}} e^{ -\nicefrac{\|\pbold\|^2}{2}}
\end{equation}
so its $\qbold$ marginal is the target distribution.

OLD and ULD are closely related. In fact, OLD is the $\gamma\to\infty$ overdamping limit of ULD after time dilation (e.g., \cite{pavliotis2014stochastic}; hence the name). However, OLD is a reversible Markov process but ULD is irreversible. The non-equilibrium statistical mechanics of both are interestingly different, however related too (see e.g., \cite{souza2019metastable}), but this paper will only focus on the convergence to statistical equilibrium.

Many celebrated approaches exist for establishing the exponential convergence (a.k.a. geometric ergodicity) of OLD, including the seminal work of \cite{roberts1996exponential}
, the ones using spectral gap (e.g., \citealp[Lemma 1]{dalalyan2017theoretical}), synchronous coupling \cite[p33-35]{villani2008optimal}\cite[Proposition 1]{durmus2019high}, functional inequalities such as Poincar\'e's inequality (PI) \cite[Theorem  4.4]{pavliotis2014stochastic} and logarithmic Sobolev inequality (LSI) \cite[Theorem  1]{vempala2019rapid}. There are also fruitful exponential convergence results for ULD, including the ones leveraging Lyapunov function \cite[Theorem  3.2]{mattingly2002ergodicity}, hypocoercivity \citep{villani2009hypocoercivity,dolbeault2009hypocoercivity, dolbeault2015hypocoercivity, roussel2018spectral}, coupling \cite[Theorem  5]{cheng2017underdamped}\cite[Theorem  1]{dalalyan2018sampling}\cite[Theorem  2.3]{eberle2019couplings}, LSI \cite[Section 3.1]{ma2019there}, modified Poincar\'e's inequality \cite[Theorem  1]{cao2019explicit}, and spectral analysis \citep{kozlov1989effective, eckmann2003spectral}. Generally speaking, due to technical difficulty related to lack of uniform ellipticity, the exponential convergence of ULD takes more effort to establish than OLD, especially when the potential $f$ is not strongly convex.

\paragraph{Discretized Process.}
However, studying the continuous processes is often not enough and they need to be discretized in order to be implemented as algorithms.
The study of asymptotic convergence of discretized OLD dates back to at least the 1990s \citep{meyn1994computable, roberts1996exponential}. The non-asymptotic convergence analysis of LMC discretization of OLD
can be found in \cite{dalalyan2017theoretical} and it shows the discretization achieves $\epsilon$ error, in total variation distance, in $\polylog(\nicefrac{d}{\epsilon^2})$ steps. Following this, iteration complexity of discretized OLD was also quantified in different metrics, $\polylog(\nicefrac{d}{\epsilon^2})$ in 2-Wasserstein distance \citep{durmus2016sampling} and $\polylog(\nicefrac{d}{\epsilon})$ in KL divergence \citep{cheng2018convergence}. For discretized ULD, one has improved  $\polylog(\nicefrac{\sqrt{d}}{\epsilon})$ iteration complexity in 2-Wasserstein distance \citep{cheng2017underdamped, dalalyan2018sampling} and   $\polylog(\nicefrac{\sqrt{d}}{\sqrt{\epsilon}})$ in KL divergence \citep{ma2019there}.

Better dimension dependence is generally conceived as a major advantage of ULD over OLD.

\section{Terminology and Notations}
The following conditions will be frequently referred to in various parts of this paper.
\begin{assumption}\label{asp:smooth} (\textbf{Standard Smoothness Condition})
    Assume the potential function $f: \mathbb{R}^d \mapsto \mathbb{R}$ is $\mathcal{C}^2$  and $L$-smooth, i.e., there exists a constant $L > 0$ such that $\forall  \bs{x}, \bs{y} \in \mathbb{R}^d$, we have
    \[
        \|\nabla f(\bs{y}) - \nabla f(\bs{x})\| \le L\|\bs{y} - \bs{x}\|.
    \]
\end{assumption}
(Note the above assumptions are equivalent to $\nabla^2 f \preceq LI$.)

\begin{assumption}\label{asp:convex} (\textbf{Convexity})
    The potential function $f: \mathbb{R}^d \mapsto \mathbb{R}$ is convex, if $\forall  \bs{x}, \bs{y} \in \mathbb{R}^d$, we have
    \[
        f(\bs{y}) \ge f(\bs{x}) + \innerprod{\nabla f(\bs{x})}{ \bs{y} - \bs{x}}.
    \]
\end{assumption}

\begin{assumption}\label{asp:strong-convex} (\textbf{Strong-convexity})
    The potential function $f: \mathbb{R}^d \mapsto \mathbb{R}$ is $m$-strongly-convex, if there exists $m > 0$ such that $\forall  \bs{x}, \bs{y} \in \mathbb{R}^d$, we have
    \[
        f(\bs{y}) \ge f(\bs{x}) + \innerprod{\nabla f(\bs{x})}{ \bs{y} - \bs{x}} + \frac{m}{2}\norm{\bs{y} - \bs{x}}^2.
    \]
\end{assumption}

Two quantities we will use to quantify convergence are $\chi^2$ divergence and 2-Wasserstein distance
\begin{align}\label{eq:chi_squared}
    &(\chi^2 \text{ divergence}) \quad && \chi^2(\mu_1 \| \mu_2) = \int \left(\frac{d\mu_1}{d\mu_2} - 1\right)^2 d\mu_2
\\
    &(\text{2-Wasserstein distance }) \quad && W_2(\mu_1, \mu_2) = \left(\inf_{\pi \in \Pi(\mu_1, \mu_2)} \mathbb{E}_{(\bs{X}, \bs{Y}) \sim \pi} \norm{\bs{X} - \bs{Y}}^2 \right)^\frac{1}{2}
\end{align}
where $\Pi(\mu_1, \mu_2)$ is the set of all couplings of $\mu_1$ and $\mu_2$.

\section{The Construction of HFHR}
HFHR is motivated by Nesterov Accelerated Gradient descent algorithm for Strongly Convex function (NAG-SC) in optimization \citep{nesterov1998introductory}.
It is obtained by formulating NAG-SC as a Hessian free high-resolution ODE (based on \cite{shi2018understanding} and \cite{alvarez2002second, attouch2020first}), lifting the high-resolution correction's coefficient as a free parameter, and adding appropriate Gaussian noises.

More precisely, let's start with NAG-SC algorithm:
\begin{align}
    \bs{x}_{k+1} &= \bs{y}_k - s \nabla f(\bs{y}_k) \label{eq_NAG_SC_1} \\
    \bs{y}_{k+1} &= \bs{x}_{k+1} + c(\bs{x}_{k+1} - \bs{x}_k) \label{eq_NAG_SC_2}
\end{align}
where $s$ is the learning rate (also known as step size), and $c = \frac{1 - \sqrt{ms}}{1 + \sqrt{ms}}$ is a constant chosen according to step size $s$ and the strong convexity coefficient $m$ of $f$, although the method also works for non-strongly-convex $f$.

A high-resolution ODE description of Equation  \eqref{eq_NAG_SC_1} and \eqref{eq_NAG_SC_2} is obtained in \cite[Section 2]{shi2018understanding}
\begin{equation}\label{eq:high_resolution_hessian}
    \ddot{\bs{y}} + \sqrt{s}\left(\frac{2(1-c)}{s(1+c)} + \nabla^2 f(\bs{y})\right)\dot{\bs{y}} + \frac{2}{1+c}\nabla f(\bs{y}) = \bs{0},
\end{equation}

which can better account for the effect of non-infinitesimal $s$ than the $s\to 0$ limit (note $c$ depends on $s$). However, in this original form, Hessian of $f$ is involved, which is expensive to compute and store especially for high-dimensional problems.

To obtain a Hessian-free high-resolution ODE description of Equation  \eqref{eq_NAG_SC_1} and \eqref{eq_NAG_SC_2},  we first turn the iteration into a `mechanical' version by introducing position variable $\qbold_k = \bs{y}_k$ and momentum variable $\bs{p}_k=\nicefrac{(\bs{y}_k-\bs{x}_k)}{h}$. Replacing $\bs{x}_{k+1}$ in \eqref{eq_NAG_SC_1} and the first $\bs{x}_{k+1}$ in \eqref{eq_NAG_SC_2} by $\qbold_{k+1}$ and $\bs{p}_{k+1}$, the second $\bs{x}_{k+1}$ in \eqref{eq_NAG_SC_2} by $\bs{q}_k - s\nabla f(\bs{q}_k)$, and the $\bs{x}_k$ in \eqref{eq_NAG_SC_2} by $\bs{q}_k$ and $\bs{p}_k$, we obtain
\[
    \begin{cases}
        \bs{q}_{k+1} = \bs{q}_k + h \bs{p}_{k+1} - s \nabla f(\bs{q}_k) \\
        \bs{p}_{k+1} = c \bs{p}_k - c \frac{s}{h} \nabla f(\bs{q}_k)
    \end{cases}
\]
Now, choose $\gamma$, $\alpha$ and $h$ as
$
    h = \sqrt{cs},\gamma = \frac{1-c}{h}, \alpha = \frac{s}{h}

$.
It is easy to see that $\gamma>0$, $\alpha >0$, then NAG-SC exactly rewrites as
\begin{equation}\label{eq_NAG_SC_rewritten}
    \begin{cases}
        \qbold_{k+1} = \qbold_k + h \pbold_{k+1} - h \alpha \nabla f(\qbold_k) \\
        \pbold_{k+1} = \pbold_k -h \gamma \pbold_k - h \nabla f(\qbold_k)
    \end{cases}.
\end{equation}
Note the technique for bypassing the Hessian without introducing any approximation is already well studied in the literature (e.g., \cite{alvarez2002second, attouch2020first}).

So far, both $h$ and $\alpha$ are actually determined by the hyperparameter $s$ of NAG-SC. However, if we now consider $\alpha$ as an independent variable (i.e., `lift' it) and let $h\to 0$, we see \eqref{eq_NAG_SC_rewritten} is a 1st-order discretization (with step size $h$) of the dynamics
\begin{equation}
\begin{cases}
    \dot{\qbold} = \pbold - \alpha \nabla f(\qbold) \\
    \dot{\pbold} = -\gamma \pbold - \nabla f(\qbold)
\end{cases}.
\label{eq_HFHRconti}
\end{equation}
Note $\alpha$, if inherited from NAG-SC, should be $\alpha=\sqrt{s/c}=\mathcal{O}(h)$, which, in a low-resolution ODE will be discarded, and this eventually leads to ULD rather than HFHR. However, we now allow it to be a free parameter and will see that $\alpha\neq \mathcal{O}(h)$ can be advantageous.

Before quantifying these advantages later on, we finish the construction by appropriately injecting Gaussian noises to Equation \eqref{eq_HFHRconti}. This is just like how OLD can be obtained by adding noise to the gradient flow. The right amount and structure of noise turn the ODE into a Markov process that can serve the purpose of sampling, and the detailed form of our noise is given by:

\begin{equation}\label{eq:HFHR}
    \begin{cases}
        d\qbold_t = \pbold_t - \alpha \nabla f(\qbold_t) + \sqrt{2\alpha} d\bs{W}_t \\
        d\pbold_t = -\gamma \pbold_t - \nabla f(\qbold_t) + \sqrt{2\gamma} d\bs{B}_t
    \end{cases}.
\end{equation}
Here $\alpha \ge 0 , \gamma > 0$ are constant parameters, and $\bs{W}_t, \bs{B}_t$ are independent standard Brownian motions in $\mathbb{R}^d$. This irreversible diffusion process will be named as \textbf{Hessian-Free High-Resolution}(\textbf{HFHR}) dynamics. We will write it as HFHR$(\alpha,\gamma)$ to emphasize the dependence on $\alpha$ and $\gamma$ when needed.

Note that HFHR can be intuitively understood as a mixture of an ULD and a rescaled OLD
\begin{equation}\label{eq:decomposition}
    \underbrace{
    \begin{cases}
        \frac{d}{dt}\qbold_t = \pbold_t \\
        \frac{d}{dt}\pbold_t = -\gamma \pbold_t - \nabla f(\qbold_t) + \sqrt{2\gamma} d\bs{B}_t
    \end{cases}}_{\text{ULD}}
    ,\quad
    \underbrace{
    \begin{cases}
        \frac{d}{dt}\qbold_t = -\alpha \nabla f(\qbold_t) + \sqrt{2\alpha} d\bs{W}_t \\
        \frac{d}{dt}\pbold_t = \bs{0}
    \end{cases}}_{\text{rescaled OLD}}
\end{equation}
As both OLD and ULD have $\pi$ as invariant distribution (for OLD, it is more precisely just the $\qbold$ marginal of $\pi$), it is not surprising that the invariant distribution of HFHR is also $\pi$ as shown in the following theorem.

\begin{theorem}\label{thm:invariant_distribution}
    $\pi$ is the invariant distribution of HFHR described in Equation \eqref{eq:HFHR}.
\end{theorem}

\begin{remark}[HFHR is not ULD+OLD]
    However, that the `decomposition' in Equation \eqref{eq:decomposition} is only formal, and the convergence of \eqref{eq:HFHR} can be very different from that of ULD or OLD. In fact, even for a linear system $\dot{x}=Ax+Bx$, its dynamics can already be very different from any finite composition of the flow maps of $\dot{x}=Ax$ and $\dot{x}=Bx$ unless $[A,B]=0$. The high nonlinearity only makes the mixture \eqref{eq:HFHR} even more different from ULD and OLD.
\end{remark}

\begin{remark}[HFHR's acceleration is not due to time rescaling]
    The OLD part in Equation \eqref{eq:decomposition} is a time-rescaled version of the original OLD in Equation \eqref{eq:LD}. If we only had the OLD part, choosing a large $\alpha$ would seemingly accelerate the convergence, but this is meaningless from an algorithmic point of view, because if one discretizes the original OLD in Equation \eqref{eq:LD} using a step size $h$, the rescaled version should use a step size $h/\alpha$. However, this is no longer the case when ULD and OLD are combined, as the nonlinear interaction between them will change both the convergence rate and stability limit in nonlinear ways. The acceleration enabled by HFHR is genuine; see Section \ref{sec:discretization} (and Appendix \ref{sec_onAlpha} too).
\end{remark}

\section{Theoretical Analysis of the Continuous HFHR Dynamics}
\label{sec:theoryConti}
In this section, we establish exponential convergence guarantees of HFHR in several different setups. The most general setup corresponds to Theorem \ref{thm:simple_result}, where the convexity of $f$ is not required and we only require the target measure to satisfy an isoperimetric inequality. However, the bound there might be loose due to its generality. When the target measure is, respectively, log-concave and log-strongly-concave, Theorems \ref{thm:main} and \ref{thm:exp_convergence_coupling} demonstrate additional acceleration of HFHR when compared to ULD.

\bigskip
\noindent
More precisely, we first show HFHR converges exponentially fast as long as $\mu$ satisfies a Poincar\'e's inequality.

\begin{theorem}\label{thm:simple_result}
    Suppose $\alpha > 0$ and the target measure $\mu$ satisfies a Poincar\'e's inequality
    \begin{equation}
            \int (g - \int g d\mu)^2 \le \frac{1}{\PIconstant{\mu}} \int \|\nabla g\|^2 d\mu
    \end{equation}
    for any $g \in C^2(\mathbb{R}^d) \cap L^2(\mathbb{R}^d, \mu)$ with some positive constant $\PIconstant{\mu}> 0$. Then we have
    \[
      \chi^2(\rho_t, \pi) \le \chi^2(\rho_0, \pi) \exp\left(-2\min\{\PIconstant{\mu}, 1\} \min\{\alpha, \gamma\}  t\right)
    \]
    where $\rho_t$ is the joint law of $(\qbold_t, \pbold_t)$ in HFHR.
\end{theorem}

\begin{remark}
One always has $\chi^2(\mu_1 \| \mu_2) \ge \text{KL}(\mu_1 \| \mu_2) \ge \frac{1}{2}\|\mu_1 - \mu_2\|^2_{\text{TV}}$, due to the relation to KL divergence (see e.g. \cite{sason2016f}) and Pinsker's inequality, hence exponential convergence in $\chi^2$ divergence also implies that in  KL divergence and total variation distance.
\end{remark}

The Poincar\'e's inequality assumption holds for a large family of measures, including log-concave ones as shown in the next two propositions.

\begin{proposition}\cite[Theorem A.19]{villani2009hypocoercivity}
  If $\displaystyle \lim_{\|\bs{x}\|\to\infty} \left(\frac{\|\nabla f(\bs{x})\|^2}{2}-\Delta f(\bs{x})\right)=\infty$, then $\mu$ satisfies the Poincar\'e's inequality.
\end{proposition}

\begin{proposition}\label{lemma:PI}
\citep{cheeger1969lower} Every log-concave measure $\mu \propto e^{-f(\bs{x})}$ satisfies the Poincar\'e's inequality.

\end{proposition}

However, $f$ doesn't have to be convex for Poincar\'e's inequality to hold, and thus Theorem \ref{thm:simple_result} is a non-convex result.

\bigskip
\noindent
We now focus on convex and $L$-smooth $f$ (i.e., under Conditions A\ref{asp:smooth} and A\ref{asp:convex}) for which we will obtain much tighter bounds. Proposition \ref{lemma:PI} ensures log-concave target distribution $\mu \propto e^{-f(\qbold)}$ satisfies PI with some positive constant $\PIconstant{\mu} \triangleq \lambda > 0$. Then:

\begin{theorem}\label{thm:main}
Assuming Conditions A\ref{asp:smooth}, A\ref{asp:convex}, and in addition that $\gamma^2 \ge \max\{2\lambda, L\}$ and $\alpha \le \frac{\gamma}{\lambda} - \frac{2}{\gamma}$, we have
\[
    \chi^2(\rho_t \| \pi) \le C e^{- (\frac{\sqrt{\lambda}}{2\gamma} + \frac{\sqrt{\lambda}}{16} \alpha )t}
\]
where $C=\left\{ \chi^2(\rho_0 \| \pi) + \Einfty \left[ \innerprod{\nablax \frac{\rho_0}{\pi}}{S\nablax \frac{\rho_0}{\pi}} \right] \right\}$ is a constant determined by the initial condition, $\rho_t$ is the law of $(\qbold_t, \pbold_t)$, $\nablax = (\nablaq, \nablap)$, and $S = \frac{1}{\gamma}\m{(\frac{2}{\gamma} + \alpha)I & I \\ I & \gamma I}$.
\end{theorem}

\begin{remark}
Assumptions on a lower bound of $\gamma$, such as $\gamma^2 > \max\{2\lambda, L\}$, were also made in existing literature \citep{dalalyan2018sampling, eberle2019couplings}. The assumption $\alpha \le \nicefrac{\gamma}{\lambda} - \nicefrac{2}{\gamma}$ is to ensure additional linear acceleration, i.e., a coefficient $\nicefrac{\sqrt{\lambda}}{16}$
to be independent of $\alpha$. For large $\alpha$, there will still be additional acceleration, not necessarily linear in $\alpha$ any more though.

\end{remark}

\noindent
When $f$ is not only convex but also $m$-strongly convex, we have the following result:
\begin{theorem}\label{thm:exp_convergence_coupling}
Assume Conditions A\ref{asp:smooth}, A\ref{asp:strong-convex} and further assume $\gamma^2 > L + m$ and $\alpha \le \frac{\gamma^2 - L - m}{m\gamma}$. Denote the law of $\qbold_t$ by $\mu_t$. Then there exists a constant $\kappa^\prime > 0$ depending only on $\alpha$ and $\gamma$, such that
\[
    W_2(\mu_t, \mu) \le \kappa^\prime e^{-(\frac{m}{\gamma} + m\alpha)t} W_2(\mu_0, \mu).
\]
The definition of $\kappa^\prime$ can be found in Appendix \ref{app:notation}.
\end{theorem}

\begin{table}
    \centering
    \begin{tabular}{ccccc}
        \toprule
        \multirow{2}{*}
        {Dynamics} & \multicolumn{2}{c}{Setup}  & {$\gamma$} & {Metric}\\
        \cmidrule{2-3} \\
        {} & log-concave & log-strongly-concave  & {}\\
        \midrule
        ULD \\
        \cite[Theorem  1]{dalalyan2018sampling} & N/A & $\frac{\sqrt{m}}{\sqrt{\kappa} + \sqrt{\kappa - 1}}$ & $2\sqrt{L}$ & $W_2$\\
        ULD \cite[Theorem  1]{cao2019explicit} & $\mathcal{O}(\sqrt{\lambda})$ & $\mathcal{O}(\sqrt{m})$ & $\sqrt{\lambda}$ & $\chi^2$ \\
        ULD \cite[Theorem  1]{ma2019there} & $\frac{\rho}{10}$ & $\frac{\rho}{10}$ & 2 & KL \\
        HFHR (Theorem \ref{thm:main}) & $\frac{\sqrt{\lambda}}{2\sqrt{L +2\lambda}} + \frac{\sqrt{\lambda}}{16}\alpha$ & $\frac{1}{2\sqrt{\kappa +2}} + \frac{\sqrt{m}}{16}\alpha$ & $\sqrt{L+2\lambda}$ & $\chi^2$ \\
        HFHR (Theorem \ref{thm:exp_convergence_coupling}) & N/A & $\frac{\sqrt{m}}{2\sqrt{\kappa}} + m\alpha$ & $2\sqrt{L}$ & $W_2$ \\
        \bottomrule
    \end{tabular}
    \caption{Comparison of convergence rate of HFHR and ULD with known dependence on parameters of dynamics. In log-strongly-concave setup, we write $m = \lambda$ due to Bakry-{\'E}mery condition \cite{bakry1985diffusions} and denote condition number $\kappa = \nicefrac{L}{m}$. $\rho > 0$ is the LSI constant assmued in \cite{ma2019there}. The column of $\gamma$ contains the values of $\gamma$ corresponding to the best rate.} \label{tab:comparison}
\end{table}

Respectively, Theorem \ref{thm:main} and \ref{thm:exp_convergence_coupling} state that HFHR converges to the target distribution exponentially fast in log-concave and log-strongly-concave setups.

Before demonstrating the advantage of HFHR over ULD, we will first need to inspect the bound for HFHR when $\alpha=0$, i.e. the ULD case, to ensure it is a tight bound for ULD. To this end, we compare our bound with several existing convergence results for ULD with known dependence on the parameters of dynamics in Table \ref{tab:comparison}.

In the log-concave setup, the $\mathcal{O}(\sqrt{\lambda})$ rate from \cite{cao2019explicit} is optimal when $\lambda \to 0$ and can be realized by isotropic quadratic potential \cite[Remark 1.1]{cao2019explicit}. This new result is enabled by assuming growth condition on the Hessian of $f$ and a compact embedding condition, in additional to PI assumption. Our result for ULD $\frac{\sqrt{\lambda}}{2\sqrt{L+2\lambda}} = \mathcal{O}(\sqrt{\lambda})$ is comparable to the optimal one in the same $\lambda\to 0$ regime as long as $L = \Theta(1)$.

Our proof of Theorem \ref{thm:main} is based on the machinery introduced in a recent work on ULD \citep{ma2019there}; however,
in \cite{ma2019there}, it is assumed that $\mu$ satisfies logarithmic Soblev inequality (LSI) with constant $\rho$, which is a stronger condition than PI \citep{deuschel2001large, wang2005generalization}, and it does not necessarily hold for generic log-concave measures \citep{bobkov1999isoperimetric}. The rate $\mathcal{O}(\rho)$, however, is not directly comparable with other results as \cite{ma2019there} works with a rescaled ULD \footnote{

How a rescaling affects convergence rate can be found in \cite[Lemma 8]{cheng2017underdamped} and \cite[Theorem  1]{dalalyan2018sampling}.

}.

In the log-strongly-concave setup, \cite[Theorem \ 1]{dalalyan2018sampling} obtained exponential convergence result in 2-Wasserstein distance with rate $\frac{\sqrt{m}}{\sqrt{\kappa} + \sqrt{\kappa - 1}}$ using a simple and elegant coupling approach, and showed this rate is optimal as it is achieved by the bivariate function $f(x,y) = \frac{m}{2}x^2 + \frac{L}{2}y^2$. In Theorem \ref{thm:exp_convergence_coupling}, we use the same coupling approach to obtain an (asymptotically) equivalent rate $\frac{\sqrt{m}}{2\sqrt{\kappa}}$.

After showing that our convergence rate results for HFHR(0, $\gamma$) are comparable to the optimal rates for ULD in many cases, the acceleration of HFHR immediately becomes evident. For example, if we set $\gamma = 2\sqrt{L}$ and push $\alpha$ to the upper bound specified in Theorem \ref{thm:exp_convergence_coupling}, we obtain an $\mathcal{O}(\sqrt{L})$ rate in the log-strongly-concave setup.

Compared with the rate in \cite{dalalyan2018sampling}, this is a speed-up of order $\kappa$.

\section{Discretization}\label{sec:discretization}

In order to have an implementable algorithm, we now discretize the time of HFHR dynamics. For simplicity, we'll work with constant step size $h$. Inspired by Strang splitting for differential equations \citep{strang1968construction, mclachlan2002splitting}, we consider a 2nd-order symmetric composition for updating over each time interval $[kh, (k+1)h]$:
\begin{equation}\label{eq:exact_strang}
    \bs{x}_{k+1} := \phi^{\frac{h}{2}} \circ \psi^h \circ \phi^{\frac{h}{2}} (\bs{x}_{k})
\end{equation}
where $\bs{x}_{k} = \begin{bmatrix} \bs{q}_{kh} \\ \bs{p}_{kh} \end{bmatrix}$, $\phi$ and $\psi$ correspond to solution flows of split SDEs, respectively given by
\begin{equation*}
    \phi:
    \begin{cases}
        d\bs{q} = \bs{p} dt \\
        d\bs{p} = -\gamma \bs{p} dt + \sqrt{2\gamma} d\bs{B}
    \end{cases}
    \qquad
    \psi:
    \begin{cases}
        d\bs{q} = -\alpha \nabla f(\bs{q}) dt + \sqrt{2\alpha} d\bs{W} \\
        d\bs{p} = -\nabla f(\bs{q}) dt
    \end{cases}
    ,
\end{equation*}
and $\phi^t(\bs{x}_0)$ and  $\psi^t(\bs{x}_0)$ mean $x$'s value after evolving the $\phi$ and $\psi$ stochastic dynamics for $t$ time with initial condition  $\bs{x}_0$.

Note that $\phi$ flow can be solved explicitly since the second equation is an Ornstein-Unlenbeck process and integrating the second equation followed by integrating the first one gives us an explicit solution
\begin{equation}\label{eq:phi}
    \begin{cases}
        \bs{q}_t  = \bs{q}_0 + \frac{1 - e^{-\gamma t}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^t \frac{1 - e^{-\gamma(t - s)}}{\gamma} d\bs{B}(s),\\
        \bs{p}_t = e^{-\gamma t}\bs{p}_0 + \sqrt{2\gamma} \int_0^t e^{-\gamma (t - s)} d\bs{B}(s).
    \end{cases}
\end{equation}
For an implementation of the stochastic integral part in Equation \ref{eq:phi}, denoting $\bs{X} = \sqrt{2\gamma} \int_0^t \frac{1 - e^{-\gamma(t - s)}}{\gamma} d\bs{B}(s)$ and $\bs{Y} = \sqrt{2\gamma} \int_0^t e^{-\gamma (t - s)} d\bs{B}(s)$, we have
\[
    \mbox{Cov}(\bs{X}, \bs{Y}) = \begin{bmatrix}  \frac{\gamma h + 4e^{-\gamma \frac{h}{2} } - e^{-\gamma h} -3 }{\gamma^2} I_d & \frac{(1 - e^{-\gamma \frac{h}{2}})^2}{\gamma} I_d \\  \frac{(1 - e^{-\gamma \frac{h}{2}})^2}{\gamma} I_d & (1 - e^{-\gamma h}) I_d  \end{bmatrix}.
\]
As mean and covariance fully determine a Gaussian distribution, $\begin{bmatrix} \bs{X} \\ \bs{Y}\end{bmatrix} = M \bs{\xi}$ where $M$ is the Cholesky decomposition of $\mbox{Cov}(\bs{X}, \bs{Y})$, $\bs{\xi}$ is a $2d$ standard Gaussian random vector, i.i.d. at each step, and $\phi^t$ can thus be exactly simulated.

However, $\psi$ flow is generally not explicitly solvable unless $f$ is a quadratic function in $\bs{q}$. We simply choose to approximate $\psi^h(\bs{x}_0)$ with one-step Euler-Maruyama integration
\[
    \psi^h(\bs{x}_0) \approx \widetilde{\psi}^h(\bs{x}_0) \text{ given by } \begin{cases}
    \bs{q}_h = \bs{q}_0 -\alpha \nabla f(\bs{q}_0) h + \sqrt{2\alpha h} \bs{\eta} \\
    \bs{p}_h = \bs{p}_0 -\nabla f(\bs{q}_0) h
    \end{cases}
\]
where $\bs{\eta}$ is a standard $d$-dimensional Gaussian random vector, again i.i.d. each time $\tilde{\psi}$ is called.

Altogether, one step of an implementable Strang's splitting of HFHR is hence
\begin{equation}\label{eq:implementable_strang}
    \phi^\frac{h}{2} \circ \widetilde{\psi}^h \circ \phi^\frac{h}{2}
\end{equation}
We call this numerical method defined by Eq.\eqref{eq:implementable_strang} the HFHR algorithm, which is summarized in Algorithm \ref{alg:HFHR}.

As $\psi$ in the Strang splitting is replaced by a 1st-order approximation $\tilde{\psi}$, it is not surprising to expect a reduced convergence order even if one considered the splitting of a deterministic system. For numerical SDEs, the notions of convergence are more subtle (e.g., \cite{milstein2013stochastic, kloeden1999, telatovich2017strong}) and there is even an order barrier, namely constant $h$ methods cannot exceed strong order 1 \citep{ruemelin1982numerical}. The following theorem shows that the HFHR algorithm has a strong error of order 1.

\begin{algorithm}[h]
\caption{HFHR Algorithm}\label{alg:HFHR}
\begin{algorithmic}[1]
\State \textbf{Input}: potential function $f$ and its gradient $\nabla f$, damping coefficients $\alpha$ and $\gamma$, step size $h$, initial condition $(\qbold_0, \pbold_0)$
\Procedure{Discretized HFHR}{$f, \nabla f, \alpha, \gamma, h, \qbold_0, \pbold_0$}
    \State $k=0$ and initialize $\begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$
    \While{not converge}
        \State Generate independent standard Gaussian random vectors $\bs{\eta}_{k+1}\in \mathbb{R}^d, \bs{\xi}^1_{k+1}, \bs{\xi}^2_{k+1} \in \mathbb{R}^{2d}$
        \State Run $\phi^\frac{h}{2}$ :  $ \begin{bmatrix} \bs{q}_1 \\ \bs{p}_1 \end{bmatrix} = \begin{bmatrix} \bs{q}_{kh} + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_{kh}  \\ e^{-\gamma \frac{h}{2}}\bs{p}_{kh} \end{bmatrix} + M \bs{\xi}^1_{k+1} $
        \State Run $\widetilde{\psi}^h$ : $ \begin{bmatrix} \bs{q}_2 \\ \bs{p}_2 \end{bmatrix} = \begin{bmatrix} \bs{q}_1 - \alpha \nabla f(\bs{q}_1) h + \sqrt{2\alpha h} \bs{\eta}_{k+1}  \\ \bs{p}_1 - \nabla f(\bs{q}_1) h \end{bmatrix} $
        \State Run $\phi^\frac{h}{2}$ :  $ \begin{bmatrix} \bs{q}_3 \\ \bs{p}_3 \end{bmatrix} = \begin{bmatrix} \bs{q}_2 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_2  \\ e^{-\gamma \frac{h}{2}}\bs{p}_2 \end{bmatrix} + M \bs{\xi}^2_{k+1} $
        \State $ \begin{bmatrix} \bs{q}_{(k+1)h} \\ \bs{p}_{(k+1)h} \end{bmatrix} \gets \begin{bmatrix} \bs{q}_3 \\ \bs{p}_3 \end{bmatrix} $
        \State $k \gets k + 1$
    \EndWhile  \label{roy's loop}
\EndProcedure
\end{algorithmic}
\end{algorithm}

\begin{theorem}\label{thm:discretization}
[\textbf{arbitrary long time discretization error of Algorithm \ref{alg:HFHR}} in $L_2$]

Under Conditions A\ref{asp:smooth}, A\ref{asp:strong-convex}, and further assume the function $\nabla \Delta f$ grows at most linearly, i.e., $\norm{\nabla \Delta f(\bs{q})} \le G\sqrt{1 + \norm{\bs{q}}^2 }, \forall \bs{q} \in \mathbb{R}^d$.

Assume without loss of generality that  $ \bs{0} \in \argmin_{\bs{x} \in \mathbb{R}^d} f(\bs{x})$. Also suppose $\gamma$ in HFHR dynamics satisfy $\gamma^2 > L$.

Then there exists $C > 0$, such that for $0 < h \le h_0  \triangleq \min\{\frac{1}{4 \kappa^\prime L^\prime}, h_1, h_2, h_3 \}$ where

\begin{align*}
    h_1 =& \frac{\sqrt{\lambda^\prime}}{4\sqrt{2} \kappa^\prime L \max\{\alpha + 1.25, \gamma + 1\} (1.92 + 2.30\alpha L)}, \\
    h_2 =& \frac{\lambda^\prime}{16\sqrt{2} \kappa^\prime (L+G) \max\{\alpha+1.25, \gamma + 1\}(1.74 + 0.71\alpha)}, \\
    h_3 =& \frac{\lambda^\prime}{8\kappa^\prime L \max\{\alpha + 1.25, \gamma + 1\} (1.92 + 2.30\alpha L) },
\end{align*}
we have
\[
    \left(\mathbb{E}\norm{\bs{x}_k - \bar{\bs{x}}_k}^2\right)^\frac{1}{2}  \le C h
\]
where $\bar{\bs{x}}_k$ is the $k$-th iterate of Algorithm \ref{alg:HFHR} with step size $h$ starting from $\bs{x}_0$, $\bs{x}_k$ is the solution of HFHR dynamics at time $kh$, starting from $\bs{x}_0$. This result holds uniformly for all $k \ge 0$ and $k$ can go to $\infty$. In particular, $C = \mathcal{O}(\sqrt{d})$ and if $\gamma - \frac{L+m}{\gamma}\ge m\alpha$, then there exist $b_1, b_2 > 0$, both independent of $\alpha$ and are of order $\mathcal{O}(\sqrt{d})$,
such that
\begin{equation}
    C \le \frac{b_1 \alpha^3 + b_2}{\frac{m}{\gamma} + m\alpha}.
    \label{eq:CdependenceOnAlpha}
\end{equation}
$\kappa^\prime, L^\prime, \lambda^\prime$ are constants depending only on $L, m, \gamma, \alpha$ and details about them can be found in Appendix \ref{app:notation}.
\end{theorem}

\begin{remark}
The linear growth (at infinity) condition on $\nabla \Delta f$ is actually not as restrictive as it appears, and in some sense even weaker than some classical condition for the existence of solutions to SDE. For example, a standard condition for ensuring the existence and uniqueness of a global solution to SDE is at most a linear growth (at infinity) of the drift (e.g., \cite{pavliotis2014stochastic}). If we consider monomial potentials, i.e.,  $f(x) = f_p(x) = x^p, p \in \mathbb{Z}_+$, then the linear growth condition on $\nabla \Delta f$ is met when $p \le 4$, whereas the classical condition for the existence of solutions holds only when $p \le 2$.
\end{remark}
\begin{remark}
The proof is based on the powerful mean-squared analysis framework (e.g., \cite{milstein2013stochastic}) in numerical SDE. We extend it to an infinite time horizon by utilizing an additional contractivity. Moreover, we carefully track all constants in order to quantify the dependence on important parameters such as step size $h$ and dimension $d$.  As a side note, the proof technique in \cite[Theorem 2]{dalalyan2018sampling} does not directly apply to Algorithm \ref{alg:HFHR} since the diffusion on position variable in HFHR lowers the regularity of sample paths, which is critical in their proof.
\end{remark}

Two results follow from Theorem \ref{thm:discretization}: Theorem \ref{thm:wasserstein} characterizes the difference between the law of position variable $\bs{q}$ produced by Algorithm \ref{alg:HFHR} and the target distribution $\mu$ in 2-Wasserstein distance. Corollary \ref{corollary:iteration complexity} gives the iteration complexity of Algorithm \ref{alg:HFHR}, i.e., a bound on the number of steps needed for sampling the target distribution with $\epsilon$ error.

\begin{theorem}\label{thm:wasserstein}
Consider the same assumption as in Theorem \ref{thm:discretization} and further assume $\gamma - \frac{L+m}{\gamma} \ge m\alpha$. If $(\bs{q}_0, \bs{p}_0) \sim \pi_0$, then there exists $h_0, C > 0$ (same as that in Theorem \ref{thm:discretization}; recall $C =\mathcal{O}(\sqrt{d})$) such that when $0 < h < h_0$, we have
\begin{equation}\label{eq:discretization_wasserstein}
     W_2(\mu_k, \mu) \le \sqrt{2} \kappa^\prime e^{-(\frac{m}{\gamma} + m\alpha) kh }W_2(\pi_0, \pi) + \sqrt{2}Ch
\end{equation}
where $\mu_k$ is the law of the $q$ marginal of the $k$-th iterate in Algorithm \ref{alg:HFHR}, and $\mu$ is the $q$ marginal of the invariant distribution $\pi$.
\end{theorem}

Inspecting the role of $\alpha$ in Equation \eqref{eq:discretization_wasserstein}, we see that $\alpha$ clearly increases the rate of exponential decay, but at the same time it can also increase the discretization error (see \eqref{eq:CdependenceOnAlpha}; assuming $h$ is fixed). However, as the following Corollary \ref{corollary:iteration complexity} and its remark will show, the net effect of having a positive $\alpha > 0$, at least for some $\alpha^\star$, is reduced iteration complexity.

\begin{corollary}\label{corollary:iteration complexity}
Consider the same assumption as in Theorem \ref{thm:discretization} and further assume $\gamma - \frac{L+m}{\gamma} \ge m\alpha$. If $(\bs{q}_0, \bs{p}_0) \sim \pi_0$, then there exists $h_0, C > 0$ (same as that in Theorem \ref{thm:discretization}; recall $C =\mathcal{O}(\sqrt{d})$) such that for any target error tolerance $\epsilon > 0$, if we choose $h = h^\star \triangleq \min\{ h_0, \frac{\epsilon}{2\sqrt{2}C}\}$, then after
\begin{equation}\label{eq:iteration_complexity}
    k=k^\star \triangleq \frac{1}{\frac{m}{\gamma} + m\alpha} \max\{\frac{1}{h_0}, \frac{2\sqrt{2}C}{\epsilon}\} \log \frac{2\sqrt{2}\kappa^\prime W_2(\pi_0, \pi)}{\epsilon}
\end{equation}
steps, we have $W_2(\mu_k, \mu) \le \epsilon$. When high accuracy is needed, i.e., $\epsilon < 2\sqrt{2}C h_0$, this iteration complexity becomes
\[
k^\star =  2\sqrt{2}\frac{C}{\frac{m}{\gamma} + m\alpha} \frac{1}{\epsilon} \log \frac{2\sqrt{2} \kappa^\prime W_2(\pi_0, \pi) }{\epsilon} = \tilde{\mathcal{O}}(\frac{\sqrt{d}}{\epsilon}).
\]
\end{corollary}

\begin{remark}\label{remark:iteration_complexity}
Recall from Theorem \ref{thm:discretization} that $C \le \frac{b_1 \alpha^3 +b_2}{\frac{m}{\gamma} + m\alpha}$ when $\gamma - \frac{L+m}{\gamma}\ge m\alpha$, so if we consider the minimizer $\alpha^\star$ of the upper bound of $\frac{C}{\frac{m}{\gamma} + m\alpha}$, namely
\[
    \alpha^\star = \argmin_{\alpha \ge 0} \frac{b_1 \alpha^3 + b_2}{ (\frac{m}{\gamma} + m\alpha)^2},
\]
we have $\alpha^\star > 0$. This suggests the possibility of having an optimal choice of $\alpha$, although the optimal value, if existent, may still not be analytically known as the bounds may not be tight. How the bound of the constant in the iteration complexity depends on $\alpha$ is illustrated in Figure \ref{fig:illustration}

\end{remark}

\begin{figure}[h!]
    \centering
    \includegraphics[scale=0.4]{img/best_alpha_illustration.pdf}
    \caption{Illustration of the effect of $\alpha$ on iteration complexity}
    \label{fig:illustration}
\end{figure}

\begin{remark}\label{remark:iteration_complexity_2}
The result of Theorem \ref{thm:wasserstein} agrees with Theorem 2 for a discretization of ULD (termed as 1st-order KLMC) in \cite{dalalyan2018sampling}, in that both results have an exponentially decaying first term due to converging continuous dynamics, and a second term corresponding to discretization error that scales linearly with $h \sqrt{d}$. Consequently, the iteration complexity of $\polylog(\frac{\sqrt{d}}{\epsilon})$ in Corollary \ref{corollary:iteration complexity} matches that in the discussion after Theorem 2 in \cite{dalalyan2018sampling}, with an improved constant for $\alpha > 0$.
It is worth mentioning that a higher-order discretization of ULD, termed as 2nd-order KLMC, was also proposed in \cite{dalalyan2018sampling}, and its reduced discretization error beautifully leads to a reduced iteration complexity of $\polylog(\frac{\sqrt{d}}{\sqrt{\epsilon}})$. However, that was at the expense of using the Hessian of $f$ in the iterations, which can be computationally expensive for high-dimensional problems. To focus on the main idea of turning NAG-SC into a sampler and the fact that the gained acceleration is not a fake one created by  time-rescaling, we do not consider high-order discretizations in this article.
\end{remark}

Remark \ref{remark:iteration_complexity} and \ref{remark:iteration_complexity_2} demonstrate that HFHR algorithm can lead to a similar bound on iteration complexity as ULD algorithm but with an improved constant, which establishes HFHR as a more efficient algorithm than ULD algorithm. This improvement also shows that the acceleration of HFHR can be carried through from continuous time to discrete time. Appendix \ref{sec_onAlpha} presents an additional case study that further demonstrates this point analytically. The same conclusion has been consistently observed in numerical experiments too.

\section{Numerical Experiments}\label{sec:experiment}
This section empirically investigates the performance of HFHR algorithm and compares it with ULD algorithm (1st-order KLMC in \cite{dalalyan2017theoretical}).

In Section \ref{subsec:simple}, we test both algorithms on a collection of target distributions corresponding to simple but representative potential $f$'s, which cover strongly convex (low/high dimensional Gaussian, small/large strongly-convex coefficient $m$), convex, and non-convex cases (bi-modal, perturbed Gaussian, and 2D Rosenbrock's function \citep{rosenbrock1960automatic}).

Sections \ref{subsec:verify_dynamics} and \ref{subsec:verify_algorithm} are devoted to systematic numerical verifications of our theoretical analyses, for which we work with quadratic potentials (corresponding to Gaussian targets) and a nonlinear potential when results become trivially true for Gaussians. The latter will also be used to demonstrate the improved iteration complexity of HFHR algorithm.

In Section \ref{subsec:bnn}, we demonstrate how HFHR performs in a downstream learning task based on a high-dimensional, non-convex, multi-modal Bayesian neural network model.
In all experiments, we use the same $\gamma$ and step size $h$ for ULD and HFHR. All experiments are conducted on a machine with a 2.20GHz Intel(R) Xeon(R) E5-2630 v4 CPU and an Nvidia GeForce GTX 1080 GPU.

\subsection{Simple Target Distributions}\label{subsec:simple}
In this subsection, we test 8 target distributions with simple, yet representative potential functions, summarized in Table \ref{tab:test_function}. They are classified into two groups, Gaussians and non-Gaussians. For Gaussian, smoothness coefficient $L$ is available, hence we take $\gamma=2\sqrt{L}$ in accordance with Table \ref{tab:comparison}. To be consistent with Theorem \ref{thm:discretization}, we measure closeness to the target distribution using $W_2$ which has closed-form expression between two Gaussians. For non-Gaussians, we empirically set $\gamma=2$ and measure sample quality by $\chi^2$ divergence with densities empiricially approximated by histograms. For the special case of $f_8(x,y)$, note approximating its density using a uniform-mesh-based histogram is either inaccurate or requiring the mesh to be very fine due to high nonconvexity, and we thus report the error in the $x$ component $|\mathbb{E}x - \mu|$ instead, where $\mu$ is the true mean of $x$-component.

Step size $h$ is tuned so that it is near the \textit{stability limit} of ULD algorithm. In each experiment, we sample 10,000 independent realizations

of each algorithm for empirical estimations.

Figure \ref{fig:simple} shows the experimental results. Each result contains two phases. In the first phase, error has a trend of smooth decrease, corresponding to the exponentially decay bound in Equation  \eqref{eq:discretization_wasserstein}; in the second phase, the curve is noisy since accuracy is saturated by discretization error. Across all experiments, the additional acceleration of HFHR over ULD is clearly seen in the first phase and can be very significant, for example in Figure \ref{fig:gaussian_100d} and \ref{fig:perturbed_1d}.

\begin{table}
\caption{Test potential functions. We use the shorthand notation $G^{d}_{m, \kappa}(\bs{x}) = \frac{m}{2}(\kappa x_d^2 + \sum_{i=1}^{d-1} x_i^2)$. Letters `S', `C' and `N' represent strongly convex, convex, and non-convex, respectively.} \label{tab:test_function}
    \centering
    \begin{tabular}{ccccccccc}
        \toprule
        $f$ & $f_1$ & $f_2$ & $f_3$ & $f_4$ & $f_5$ & $f_6$ & $f_7$ & $f_8$\\
        Expresion & $\frac{1}{2}x^2$ & $G^2_{0.1,10}$ & $G^2_{10,10}$ & $G^{100}_{1,100}$
        & $\frac{1}{4}x^4$ & $\frac{5x^2 + \sin(10x)}{10}$ & $5(x^4 - 2x^2)$ & $\frac{(x-1)^2 + 10(y-x^2)^2}{2}$\\
        Convexity & S & S & S & S & C & N & N & N \\
        \bottomrule
    \end{tabular}
\end{table}

\begin{figure}
\centering
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gaussian_1d_splitting.pdf}
		\caption{1D Gaussian} \label{fig:gaussian_1d}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gaussian_2d_small_m_splitting.pdf}
		\caption{2D Gaussian($m=0.1$)} \label{fig:gaussian_2d_small_m}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gaussian_2d_large_m_splitting.pdf}
		\caption{2D Gaussian($m=10$)} \label{fig:gaussian_2d_large_m}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gaussian_100d_splitting.pdf}
		\caption{100D Gaussian} \label{fig:gaussian_100d}
	\end{subfigure}

    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/convex_1d_splitting.pdf}
		\caption{Convex} \label{fig:convex_1d}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/perturbed_1d_splitting.pdf}
		\caption{Perturbed} \label{fig:perturbed_1d}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/bimodal_1d_splitting.pdf}
		\caption{Bi-modal} \label{fig:bimodal_1d}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/rosenbrock_2d_splitting.pdf}
		\caption{2D Rosenbrock} \label{fig:rosenbrock_2d}
	\end{subfigure}

    \caption{(a) $f_1\, (h=2)$. (b) $f_2\, (h=0.2)$. (c) $f_3\, (h=2.5)$. (d)$f_4\, (h=0.2)$. (e) $f_5\, (h=0.5)$. (f) $f_6\, (h=0.001)$. (g) $f_7\, (h=0.1)$. (h) $f_8\, (h=0.005)$.  $y$-axes are in log scale.

    }\label{fig:simple}

\end{figure}

\subsection{A Case Study on Gaussian: Empirical Performances versus Theoretical Guarantees for HFHR Dynamics (Continuous Time)}\label{subsec:verify_dynamics}
In this subsection, we numerically verify the tightness of theoretical results for the continuous HFHR dynamics in Section \ref{sec:theoryConti}.
To this end, we work with Gaussian distributions because (i) Gaussians have log-strongly-concave densities and hence satisfy all the assumptions in our theory; (ii) Gaussians are analytically tractable, all parameters needed for the theory, including Lipschitz constant $L$, strongly-convex coefficient $m$, etc., are known to us, also the metrics used in our theory (e.g. 2-Wasserstein distance, $\chi^2$ divergence) have closed form expression for Gaussians; (iii) Gaussians are widely used in various machine learning studies and applications.

The first verification starts with Theorem \ref{thm:simple_result}. We experiment with 1-dimensional standard Gaussian distribution and use Algorithm \ref{alg:HFHR} with tiny step size ($h=10^{-4}$, small enough so that results about the continuous dynamics can be probed) to simulate HFHR dynamics and compare with the theoretical bound in Theorem \ref{thm:simple_result}. We set $\gamma=2$ which corresponds to the critical damping of ULD and generate 10,000 samples to estimate $\chi^2$ divergence to the target distribution.

\begin{figure}[h]
\centering
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_4_alpha=1e-1.pdf}
		\caption{$\alpha=0.1$} \label{fig:verify_theorem_4_alpha=0.1}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_4_alpha=5e-1.pdf}
		\caption{$\alpha=0.5$} \label{fig:verify_theorem_4_alpha=0.5}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_4_alpha=1.pdf}
		\caption{$\alpha=1$} \label{fig:verify_theorem_4_alpha=1}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_4_alpha=large.pdf}
		\caption{$\alpha \geq \gamma $} \label{fig:verify_theorem_4_alpha_large}
	\end{subfigure}

    \caption{ Illustration of the consistency between the theoretical bound in Theorem \ref{thm:simple_result} and experiment results.
    } \label{fig:verify_theorem_4}

\end{figure}

We run Algorithm \ref{alg:HFHR} with $\alpha \in \{0.1, 0.5, 1, 2, 5, 10\}$ and plot its results along with the upper bound in Theorem \ref{thm:simple_result} in Figure \ref{fig:verify_theorem_4}. In all four subfigures, HFHR outperforms ULD in terms of convergence speed, even when $\alpha$ is small (e.g. Figure \ref{fig:verify_theorem_4_alpha=0.1}), and the gap between convergence speeds becomes larger as we increase $\alpha$ (e.g., Figure \ref{fig:verify_theorem_4_alpha=0.1} $\sim$ \ref{fig:verify_theorem_4_alpha=1}). Note also the bound in Theorem \ref{thm:simple_result} may not be tight due to its generality. For example, the rate of the exponential convergence bound in Theorem \ref{thm:simple_result} reads $-2\min\{\alpha, \gamma\}$ and is $-2\gamma$ whenever $\alpha \ge \gamma$, but as Figure \ref{fig:verify_theorem_4_alpha_large} shows, when  $\alpha \ge 2$, larger $\alpha$ still introduces further acceleration in practice (in the $\alpha=2$ case, the bound appears to be rather tight though). On a side note, we observe that the convergence of HFHR and ULD in $\chi^2$ divergence consists of two phases, the convergence in the first phase is super-linear and gradually transitions to roughly linear in the second phase. Theorem \ref{thm:simple_result} characterizes the linear convergence rate in the second phase.

Next we empirically verify Theorem \ref{thm:main}. The experiment setup is identical to that for Figure \ref{fig:verify_theorem_4}. The results are plotted in Figure \ref{fig:verify_theorem_8} along with the theoretical bounds in Theorem \ref{thm:main} for both HFHR and ULD. In all four plots, HFHR converges faster than ULD. In particular, in Figure \ref{fig:verify_theorem_8_alpha=5} and \ref{fig:verify_theorem_8_alpha=10}, HFHR reaches $10^{-3}$-accuracy in $\chi^2$ divergence in less than 0.75 unit time while ULD just reaches $10^{-1}$-accuracy in the experiment time frame. In addition, we make two observations: (i) In Figure \ref{fig:verify_theorem_8_alpha=0.1} and \ref{fig:verify_theorem_8_alpha=1}, the two upper bounds for HFHR and ULD respectively in Theorem \ref{thm:main} are hardly distinguishable because $\alpha$ is small, however, the actual acceleration created by $\alpha$ is clearly manifested, even when $\alpha$ is tiny, i.e., $\alpha=0.1$ in Figure \ref{fig:verify_theorem_8_alpha=0.1}. Therefore, in practice, HFHR dynamics can introduce even more acceleration than Theorem \ref{thm:main} ensures. (ii) Note that the Lipschitz constant $L=1$ and Poincar\'e constant $\lambda=1$, so the assumption $\gamma^2 \ge \max\{2\lambda, L\}$ is satisfied and
the assumption $\alpha \le \frac{\gamma}{\lambda} - \frac{2}{\gamma}$ is satisfied when $\alpha \le 1$. In Figure \ref{fig:verify_theorem_8_alpha=5} and \ref{fig:verify_theorem_8_alpha=10}, $\alpha$ is too large to satisfy the assumption, but it still creates significant acceleration. This observation demonstrates the applicability of HFHR is wider than conditions needed in Theorem \ref{thm:main}. In fact, Theorem \ref{thm:main} only requires convexity, but Gaussian targets correspond to strongly convex potentials, so its bound being loose should not be surprising.

\begin{figure}[h]
\centering
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_8_alpha=1e-1.pdf}
		\caption{$\alpha=0.1$} \label{fig:verify_theorem_8_alpha=0.1}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_8_alpha=1.pdf}
		\caption{$\alpha=1$} \label{fig:verify_theorem_8_alpha=1}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_8_alpha=5.pdf}
		\caption{$\alpha=5$} \label{fig:verify_theorem_8_alpha=5}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_8_alpha=10.pdf}
		\caption{$\alpha =10$} \label{fig:verify_theorem_8_alpha=10}
	\end{subfigure}

    \caption{ Illustration of the consistency between the theoretical bound in Theorem \ref{thm:main} and experiment results.  }\label{fig:verify_theorem_8}

\end{figure}

Finally, we verify Theorem \ref{thm:exp_convergence_coupling}. The experimental setup is the same as that for Figure \ref{fig:verify_theorem_4}.  We plot the $W_2$ errors of HFHR and ULD, together with their respective upper bounds in Theorem \ref{thm:exp_convergence_coupling} in Figure \ref{fig:verify_theorem_10}. We observe acceleration of HFHR over ULD in all four subfigures of Figure \ref{fig:verify_theorem_10}. In addition, our theoretical upper bounds are nearly parallel to the performance of actual HFHR algorithm, which empirically demonstrates that our theory in Theorem \ref{thm:exp_convergence_coupling} is tight up to a constant. Note that the two assumptions in Theorem \ref{thm:exp_convergence_coupling}, namely $\gamma^2 > L + m$ and $\alpha \le \frac{\gamma^2 - L - m}{m\gamma}$.
In Figures \ref{fig:verify_theorem_10_alpha=2} and \ref{fig:verify_theorem_10_alpha=10}, the values of $\alpha$ break the assumption; nevertheless, Algorithm \ref{alg:HFHR} empirically works well and creates significant acceleration over ULD.

\begin{figure}[h]
\centering
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_10_alpha=5e-1.pdf}
		\caption{$\alpha=0.5$} \label{fig:verify_theorem_10_alpha=0.5}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_10_alpha=1.pdf}
		\caption{$\alpha=1$} \label{fig:verify_theorem_10_alpha=1}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_10_alpha=2.pdf}
		\caption{$\alpha=2$} \label{fig:verify_theorem_10_alpha=2}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_theorem_10_alpha=10.pdf}
		\caption{$\alpha =10$} \label{fig:verify_theorem_10_alpha=10}
	\end{subfigure}

    \caption{Illustration of the consistency between the theoretical bound in Theorem \ref{thm:exp_convergence_coupling} and experiment results.  }\label{fig:verify_theorem_10}

\end{figure}

\subsection{A Case Study on a Nonlinear Problem: Empirical Performances versus Theoretical Guarantees for HFHR Algorithm (Discrete Time)}\label{subsec:verify_algorithm}

We now numerically verify the tightness of theoretical results in Section \ref{sec:discretization}.  For this purpose, we no longer use Gaussian targets for the following reason:

A key feature of HFHR algorithm is that both its discretization error and iteration complexity depends on the ambient dimension like $\mathcal{O}(\sqrt{d})$ (similiar to ULD algorithm). This will be trivially true for standard Gaussian, because in this case HFHR will decouple across different dimensions, and hence its discretization error having a $\mathcal{O}(\sqrt{d})$ dependence would just be a consequence of using 2-Wasserstein distance for quantifying statistical accuracy.

Therefore, to inspect a more interesting example, we consider a potential that is no longer additive across different dimensions, namely
c
It's not hard to see that all dimensions will be coupled in the HFHR dynamics (and ULD too). Moreover, the new potential $f$ is still a strongly convex function and satisfies the assumption in Theorem \ref{thm:wasserstein}. When the target measure is non-Gaussian, we no longer have a closed form expression for 2-Wasserstein distance and it is computationally expensive to approximate the 2-Wasserstein distance by samples. Therefore, we use the error of mean instead as a surrogate because
\[
    \displaystyle
    \norm{\mathbb{E}_{\mu_k} \bs{q} - \mathbb{E}_{\mu} \bs{q} } \le W_2(\mu_k, \mu)
\]
and hence the bound in Equation \eqref{eq:discretization_wasserstein} also applies to the error in mean, and so does the iteration complexity bound in Equation \eqref{eq:iteration_complexity}.

We first verify Theorem \ref{thm:wasserstein}, which says the sampling error is upper bounded by two terms, the first corresponding to the exponential decay of the continuous dynamics, whose rate is characterized in Theorem \ref{thm:exp_convergence_coupling} and already numerically verified previously, and the second term corresponding to discretization error. We theoretically proved in Theorem \ref{thm:wasserstein} that the discretization error is linear in step size $h$ and square root of dimension $\sqrt{d}$.

To numerically verify the linear dependence on $h$, we work with $d=2$ and ran ULD algorithm for sufficiently long time with tiny step size ($h=0.0005$) to obtain $10^8$ independent realizations and use them (as benchmark) to empirically estimate $\mathbb{E}_\mu \bs{q}$. We then set $\gamma = 2, \alpha=1$ and
use Algorithm \ref{alg:HFHR}, with $h \in \{2^k| -7 \le k \le 0\}$. For each $h$, we run $\frac{T}{h}$ (with $T=50$) iterations in Algorithm \ref{alg:HFHR} to ensure the Markov chains are well-mixed and the contribution to final error from exponential decay is order-of-magnitude smaller than discretization error. The results are plotted in Figure \ref{fig::verify_corollary_14_h}.

The observed linear dependence on the step size $h$ is consistent with Theorem \ref{thm:wasserstein}.

\begin{figure}[h]
    \centering
    \begin{subfigure}{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_dependence_on_h.pdf}
    	\caption{Linear dependence of discretization error of Algorithm \ref{alg:HFHR} on $h$} \label{fig::verify_corollary_14_h}
	\end{subfigure}
    \begin{subfigure}{0.48\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/verify_dependence_on_d.pdf}
		\caption{Linear dependence of discretization error of Algorithm \ref{alg:HFHR} on $\sqrt{d}$ } \label{fig::verify_corollary_14_d}
	\end{subfigure}
    \caption{ Illustration of the consistency between the theoretical bound in Theorem \ref{thm:wasserstein} and experiment results.  }

\end{figure}

To numerically verify the $\mathcal{O}(\sqrt{d})$ dependence on dimension $d$, we extensively experiment with $d \in \left\{1, 2, 5, 10, 20, 50, 100, 200, 500, 1000\right\}$. For each $d$, we run 1,000 independent realizations of ULD algorithm until well converged
with tiny step size ($h=0.005$) and use their empirical average as the `true' mean. Then we fix $\gamma=2, \alpha=1, h=0.1, T=10$ and for each $d$, we run 1,000 independent realizations of HFHR algorithm for $\frac{T}{h}=100$ iterations. Experiment results are plotted  in Figure \ref{fig::verify_corollary_14_d}. The linear dependence of error on $\sqrt{d}$ shown in Figure \ref{fig::verify_corollary_14_d} demonstrates that the bound in Theorem \ref{thm:wasserstein} is tight in terms of $d$-dependence.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.7\textwidth]{img/verify_iteration_complexity.pdf}
	\caption{Improvement of Algorithm \ref{alg:HFHR} over ULD algorithm in iteration complexity. (vertical bar stands for one standard deviation.) } \label{fig::verify_iteration_complexity}
\end{figure}

The final experiment in this subsection is to compare Algorithm \ref{alg:HFHR} with ULD algorithm in terms of iteration complexity. The goal is to demonstrate that the genuine acceleration of HFHR is not an artifact due to time rescaling, which would disappear after discretization as the stability limit changes accordingly. To do so, we push both ULD and HFHR to their respective largest $h$ values that still allow monotonic convergence at a large scale, and compare their mixing times. For Gaussian targets, this is already analytically studied in Appendix \ref{sec_onAlpha} (see e.g., Fig.\ref{fig_discountFactors}). Thus, we again work with the potential in Equation \eqref{eq:nonlinear-potential}. For general nonlinear problems like this one,
Remark \ref{remark:iteration_complexity} and the illustration of Figure \ref{fig:illustration} suggest that with appropriately chosen $\alpha$, HFHR algorithm can effectively reduce the constant factor of the iteration complexity, which implies accelerated sampling. To just provide one empirical verification of this improvement over ULD algorithm, we choose $d=10$ and use the error of mean $\norm{\mathbb{E}_{\mu_k} \bs{q} - \mathbb{E}_{\mu} \bs{q} } $ to measure sampling accuracy. The benchmark, i.e., $\mathbb{E}_{\mu} \bs{q}$, is again obtained from 1,000 independent realizations of ULD algorithm with tiny step size $(h=0.005)$, ran for long enough to ensure the corresponding Markov chain is well-mixed. The initial measure is chosen as a Dirac measure at $(\bs{1}_d, \bs{0}_d)$, where $\bs{1}_d, \bs{0}_d$ are $d$-dimensional vectors filled with $1$ and $0$ respectively.  We pick threshold $\epsilon=0.1$, and for each $\alpha \in \{0, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100\}$, we try all combinations of $(\gamma, h) \in \{0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100\} \times \{5, 1, 0.5, 0.1, 0.05, 0.01, 0.005\}$ for Algorithm \ref{alg:HFHR} (which degenerate to an ULD discretization when $\alpha=0$), and empirically find the best combination that requires the fewest iterations to meet $\norm{\mathbb{E}_{\mu_k} \bs{q} - \mathbb{E}_{\mu} \bs{q} }  \le \epsilon$. We find that $h=5$ already surpasses the stability limit of ULD algorithm, hence the range of step size in the experiment covers the largest step size that can be practically used for ULD algorithm.  We repeat the experiments with 10 different random seeds for better statistical significance. The results are shown in Figure \ref{fig::verify_iteration_complexity}. When $\alpha > 0$, HFHR algorithm consistently outperforms ULD algorithm. In particular, when $\alpha=1$, which is empirically the best value we found for this experiment, HFHR algorithm achieves the specified $\epsilon$-closeness nearly 6$\times$ times faster than ULD algorithm. This empirical study corroborates that the acceleration HFHR dynamics creates also carries through to its discretization, and the acceleration of HFHR algorithm over ULD algorithm can be significant.

\subsection{Bayesian Neural Network}\label{subsec:bnn}
We now consider a Bayesian neural network (BNN) which is a compelling learning model \citep{wilson2020case}; however, the focus will not be on its learning capability, and instead we will just consider its training, which amounts to a practical, high-dimensional, and multi-modal example of sampling tasks. The specific network considered here has a fully-connected architecture with [22, 10, 2] neurons in each layer and ReLU activation. Standard Gaussian prior is used for all parameters. We compare both ULD and HFHR on Parkinson's disease data set

from UCI machine learning repository \citep{Dua:2019}.

Choices of hyper-parameter for Algorithm \ref{alg:HFHR} and ULD algorithm are systematically investigated. For each pair $(\gamma, \alpha) \in \{0.1, 0.5, 1,  5, 10,  50, 100\}^2$, we empirically tune the step size to the stability limit of ULD algorithm, simulate 10,000 independent realizations of the Markov chain, and use the ensemble to conduct Bayesian posterior prediction. HFHR will then use the same step size. For each $\gamma$, we plot the negative log likelihood of HFHR algorithm (with different $\alpha$ choices) and ULD algorithm on training and test data in Figure \ref{fig:bnn_2}. Cases where $\alpha$ is too large for numerical stability are not drawn.

From Figure \ref{fig:bnn_2}, we find that HFHR converges significantly faster than ULD in a wide range of setups. In general, the (strongly) log-concave assumption required in Theorem  \ref{thm:main} and Theorem  \ref{thm:discretization} do not hold for multimodal target distributions, and even the Poincar\'e inequality needed in the nonconvex result of Theorem \ref{thm:simple_result} may not be verified for complex models. However, this numerical result shows that HFHR still accelerates ULD for highly complex models such as BNN, even when there is no obvious theoretical guarantee.

It showcases the applicability and effectiveness of HFHR as a general sampling algorithm.

\begin{figure}[h!]
\centering
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=1e-1_h=5e-3_train_nll.pdf}
		\caption{$\gamma=0.1 \, (h=0.005)$}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=1_h=1e-2_train_nll.pdf}
		\caption{$\gamma=1 \,  (h=0.01)$}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=10_h=5e-2_train_nll.pdf}
		\caption{$\gamma=10 \, (h=0.05)$}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=100_h=1e-1_train_nll.pdf}
		\caption{$\gamma=100 \, (h=0.1)$}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=1e-1_h=1e-2_train_nll.pdf}
		\caption{$\gamma=0.1 \, (h=0.01)$} \label{fig:r1c2}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=1_h=2e-2_train_nll.pdf}
		\caption{$\gamma=1 \, (h=0.02)$} \label{fig:r2c2}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=10_h=1e-1_train_nll.pdf}
		\caption{$\gamma=10 \, (h=0.1)$}
	\end{subfigure}
    \begin{subfigure}{0.24\textwidth}
		\centering
		\includegraphics[width=\textwidth]{img/gamma=100_h=2e-1_train_nll.pdf}
		\caption{$\gamma=100 \, (h=0.2)$}
	\end{subfigure}

    \caption{ Training Negative Log-Likelihood (NLL) for various $\gamma$. Top row: step sizes are below the stability limit of ULD algorithm; Bottom row: a further increased step size would go above the stability limit of ULD algorithm }\label{fig:bnn_2}

\end{figure}

\section{Conclusion and Discussion}
This paper proposes HFHR, an accelerated gradient-based MCMC method for sampling. To demonstrate the acceleration enabled by HFHR, the geometric ergodicity of HFHR (both the continuous and the discretized versions) is quantified, and its convergence is provably faster than Underdamped Langevin Dynamics, which by itself is often already considered as an acceleration of Overdamped Langevin Dynamics. As HFHR is based on a new perspective, which is to turn NAG-SC optimizer with \textbf{finite} learning rate into a sampler, there are a number of interesting directions in which this work can be extended. Besides further theoretical investigations that aim at refining the error bounds, examples also include the followings: to scale HFHR up to large data sets, full gradient may be replaced by stochastic gradient (SG) --- how to quantify, and hence optimize the performance of SG-HFHR? Can the generalization ability of HFHR trained learning models (e.g., BNN) be quantified, and how does it compare with that by LMC and/or KLMC? We plan to study these in follow-up works.

\acks{
The authors sincerely thank Michael Tretyakov and Yian Ma for helpful discussions. MT was partially supported by NSF DMS-1847802 and ECCS-1936776.
}

\newpage
\bibliography{reference.bib}

\newpage
\appendix
\section{Notations}\label{app:notation}
We introduce a few notations that are used in the main text as well as some proof. When $\nabla f$ is $L$-Lipschitz, the drift term $\begin{bmatrix} \bs{p} - \alpha \nabla f(\bs{q}) \\ -\gamma \bs{p} - \nabla f(\bs{q})\end{bmatrix}$ in HFHR dynamics is also $L^\prime$-Lipschitz, as proved in Lemma \ref{lemma:lipschitz}, where
\[
    L^\prime = \sqrt{2} \max\left\{ \sqrt{1 + \alpha^2} \max\left\{\frac{1}{\sqrt{2}}, L \right\} , \sqrt{1 + \gamma^2}\right\}.
\]

We show in Lemma \ref{lemma:contraction} that a linear-transformed HFHR dynamics satisfies the nice contraction property, the linear transformation $P$ we use is defined as
\[
    P = \begin{bmatrix} \gamma I & I \\ 0 & \sqrt{1+\alpha\gamma}I\end{bmatrix} \in \mathbb{R}^{2d \times 2d}.
\]
Denote the largest and the smallest singular value of $P$ by
\begin{align*}
    \sigma_{\text{max}} =& \sqrt{\frac{\alpha \gamma}{2}+\frac{\gamma^{2}}{2} + \frac{\sqrt{\alpha^{2} \gamma^{2}-2 \alpha \gamma^{3}+4 \alpha \gamma+\gamma^{4}+4}}{2}+1}, \\
    \sigma_{\text{min}} =&s \sqrt{\frac{\alpha \gamma}{2}+\frac{\gamma^{2}}{2}-\frac{\sqrt{\alpha^{2} \gamma^{2}-2 \alpha \gamma^{3}+4 \alpha \gamma+\gamma^{4}+4}}{2}+1}
\end{align*}
and its condition number by
\[
    \kappa^\prime = \frac{\sigma_{\text{max}}}{\sigma_{\text{min}}} = \sqrt{\frac{\frac{\alpha \gamma}{2}+\frac{\gamma^{2}}{2} + \frac{\sqrt{\alpha^{2} \gamma^{2}-2 \alpha \gamma^{3}+4 \alpha \gamma+\gamma^{4}+4}}{2}+1}{\frac{\alpha \gamma}{2}+\frac{\gamma^{2}}{2}-\frac{\sqrt{\alpha^{2} \gamma^{2}-2 \alpha \gamma^{3}+4 \alpha \gamma+\gamma^{4}+4}}{2}+1}}.
\]
The rate $\lambda^\prime$ of exponential convergence of transformed HFHR dynamics is characterized in Lemma \ref{lemma:contraction} and is defined as
\[
    \lambda^\prime = \min\left\{ \frac{m}{\gamma} + \alpha m, \frac{\gamma^2 - L}{\gamma} \right\}
\]
given that $\gamma^2 > L$.

\section{Proof of Theorem \ref{thm:invariant_distribution}}
\begin{proof}
The Fokker-Plank equation of HFHR is given by
\[
    \partial_t \rho_t = - \nablax \cdot \left(\m{\pbold \\ -\nabla f(\qbold)} \rho_t \right)
    + \alpha\left(\nablaq \cdot ( \nabla f(\qbold) \rho_t ) + \Delta_{\qbold} \rho_t\right)
    + \gamma\left(\nablap \cdot ( \pbold \rho_t ) + \Delta_{\pbold} \rho_t\right)
\]
where $\nablax = (\nablaq, \nablap)$. For $\pi \propto e^{-f(\qbold) - \frac{1}{2} \|\pbold\|^2}$, we have
\begin{align*}
    & \nablax \cdot \left(\m{\pbold \\ -\nabla f(\qbold)} \pi \right) = \innerprod{\m{\pbold \\ -\nabla f(\qbold)}}{\nablax \pi} = 0, \\
    & \Delta_{\qbold} \pi = -\nablaq \cdot(\pi \nabla f(\qbold)) \\
    & \Delta_{\pbold} \pi = -\nablap \cdot(\pi \pbold)
\end{align*}
Therefore $\partial_t \pi = 0$ and hence $\pi$ is the invariant distribution of HFHR.
\end{proof}

\section{Proof of Theorem \ref{thm:simple_result}}
\begin{proof}
    The Fokker-Planck equation of HFHR is given by
    \begin{equation}\label{eq:FP}
         \partial_t \rho_t + \nabla \cdot (\rho_t\bs{J}) = 0, \quad \text{where }         \bs{J} = \left (
            \begin{bmatrix}
            \pbold - \alpha \nabla f(\bs{q}) - \alpha \nabla_{\bs{q}} \log \rho_t \\
            -\gamma \bs{p} - \nabla f(\bs{q}) - \gamma  \nabla_{\bs{p}} \log \rho_t
            \end{bmatrix}
                \right )
    \end{equation}
    Since
    \[
        \nabla \cdot \left(\rho_t \begin{bmatrix} -\nabla_{\bs{p}} \log \rho_t \\ \nabla_{\bs{q}} \log \rho_t \end{bmatrix}\right) = 0,
    \]
    we can then rewrite $\bs{J}$ as
    \[
    \bs{J} = \frac{\pi}{\rho_t}
             \begin{bmatrix}
            -\alpha I & I\\
            -I & -\gamma I
            \end{bmatrix}
            \begin{bmatrix}
                \nabla_{\qbold} \frac{\rho_t}{\pi} \\
                \nabla_{\pbold}  \frac{\rho_t}{\pi}
            \end{bmatrix}.
    \]

    By Lemma \ref{lemma:PI_joint}, $\pi$ satisfies PI with constant $\lambda_{\text{PI}}(\pi) = \min\{\PIconstant{\mu}, 1\}$ as it is well known that $\nu$ satisfies Poincar\'e's inequality with $\lambda_{\text{PI}}(\nu) = 1$ \cite[Theorem  3.20]{boucheron2013concentration}.
    The time derivative of $\chi^2(\rho_t, \pi)$ is
    \[
        \frac{d}{dt} \chi^2(\rho_t, \pi)
      = - \int 2(\frac{\rho_t}{\pi} - 1) \, \nabla \cdot (\rho_t \bs{J}) \,d\bs{x}
      = 2\int \langle \nabla \frac{\rho_t}{\pi}, \bs{J} \rangle \rho_t \,d\bs{x}
      \le -2\min\{\alpha, \gamma\} \int \|\nabla \frac{\rho_t}{\pi}\|^2 d\pi
    \]

    Let $g = \frac{\rho_t}{\pi} - 1$, by Poincar\'e's inequality, we have
    $
        \int \|\frac{\rho_t}{\pi} - 1\|^2 d\pi \le \frac{1}{\lambda_{\text{PI}}(\pi)} \int \|\nabla  \frac{\rho_t}{\pi}\|^2 d\pi.
    $
    Therefore we have the time derivative of $\chi^2(\rho_t, \pi)$ is bounded by
    \[
        \frac{d}{dt}\chi^2(\rho_t, \pi) \le -2\lambda_{\text{PI}}(\pi) \min\{\alpha, \gamma\}  \chi^2(\rho_t, \pi)
        = - 2\min\{\PIconstant{\mu}, 1\} \min\{\alpha, \gamma\}  \chi^2(\rho_t, \pi)
    \]
    and the desired result follows by Gronwall's inequality.
\end{proof}

\section{Proof of Theorem \ref{thm:main}}
\begin{proof}
    Two main tools we use in this proof are
    \begin{itemize}
        \item a carefully-crafted Lyapunov function, motivated by \cite{ma2019there} and
        \item the Poincar\'e's inequality of the joint invariant distribution $\pi$.
    \end{itemize}

    More specifically, we will consider the following Lypunov function
    \begin{equation}\label{eq:augmented_lyapunov}
        \mathcal{L}(\rho_t) = \chi^2(\rho_t \| \pi) + \Einfty [\innerproduct{\nablax \frac{\rho_t}{\pi}}{S \nablax \frac{\rho_t}{\pi}}]
    \end{equation}
    where $\nablax = (\nablaq, \nablap)$ and $S = \m{aI & bI \\ bI & dI} \in \mathbb{R}^{2d \times 2d}$ is a positive definite matrix to be determined later. Denote
    \begin{equation}\label{eq:cross_term}
        \crossterm = \Einfty [\innerproduct{\nablax \frac{\rho_t}{\pi}}{S \nablax \frac{\rho_t}{\pi}}].
    \end{equation}

    It is well known that the standard Gaussian measure $\nu$ satisfies PI with PI constant $\PIconstant{\nu} = 1$ \cite{boucheron2013concentration}. Therefore, by Lemma \ref{lemma:PI_joint}, the joint invariant distribution $\pi$ satisfies PI with PI constant $\PIconstant{\pi} = \min\{1, \PIconstant{\mu}\} = \min\{1, \lambda\}$. This dependence on $\lambda$, however, is undesirable when $\lambda \gg 1$ because it does not reflect fast convergence of HFHR.

    In order to unify the two cases $\lambda > $ and $\lambda < 1$, we will work with a rescaled version of HFHR. To this end, we need to introduce a larger class of dynamics parametrized by an inverse temperature parameter $\beta$
    \begin{equation}\label{eq:tempered_HFHR}
        \begin{cases}
            d\qbold = \left( \pbold - \alpha \nabla f(\qbold) \right)dt + \sqrt{2\alpha \beta^{-1}} d\bs{B}^1_t \\
            d\pbold = \left( -\gamma \pbold - \nabla f(\qbold) \right)dt + \sqrt{2\gamma \beta^{-1}} d\bs{B}^2_t
        \end{cases}.
    \end{equation}
    We refer it as tempered HFHR and denote it by tempered-HFHR($\alpha, \gamma, \beta$). It is easy to see that the plain HFHR described in Equation  \eqref{eq:HFHR} is a tempered HFHR with $\beta = 1$, i.e. tempered-HFHR$(\alpha, \beta, 1)$.

    \paragraph{Rescaling}
    Since $\mu \propto e^{-f(\qbold)}$ satisfies PI with PI constant $\lambda$, it is easy to see that $\tilde{\mu} \propto e^{-\tilde{f}(\qbold)}$ where $f(\qbold) = \lambda \tilde{f}(\qbold)$, satisfies PI with PI constant 1.

    With the rescaled potential $\tilde{f}$, HFHR process rewrites as
    \[
        \begin{cases}
            d\qbold = (\pbold - \alpha \lambda \nabla \tilde{f}(\qbold)) dt + \sqrt{2\alpha} d\bs{B}_t^1 \\
            d\pbold = (-\gamma \pbold - \lambda \nabla \tilde{f}(\qbold)) dt + \sqrt{2\gamma} d\bs{B}_t^2
        \end{cases}
    \]
    Introduce rescaled velocity $\tpbold$ via  $\qbold(t)=\sqrt{\lambda}\tpbold(t)$, then the SDE becomes
    \[
        \begin{cases}
            d\qbold = (\sqrt{\lambda}\tpbold - \alpha \lambda \nabla \tilde{f}(q)) dt + \sqrt{2\alpha} d\bs{B}_t^1 \\
            d\tpbold = (-\gamma \tpbold - \sqrt{\lambda} \nabla \tilde{f}(q)) dt + \sqrt{2\gamma/\lambda} d\bs{B}_t^2
        \end{cases}
    \]
    Introduce rescaled dissipation parameters $\tilde{\alpha},\tilde{\gamma}$ via
    $\alpha=\tilde{\alpha}\sqrt{\lambda}^{-1}$, $\gamma=\tilde{\gamma}\sqrt{\lambda}$. Then the SDE rewrites as
    \[
        \begin{cases}
            d\qbold = (\sqrt{\lambda}\tpbold - \tilde{\alpha} \sqrt{\lambda} \nabla \tilde{f}(\qbold)) dt + \sqrt{2\tilde{\alpha}\sqrt{\lambda}^{-1}} d\bs{B}_t^1 \\
            d\tpbold = (-\tilde{\gamma} \sqrt{\lambda} \tpbold - \sqrt{\lambda} \nabla \tilde{f}(\qbold)) dt + \sqrt{2\tilde{\gamma}\sqrt{\lambda}^{-1}} d\bs{B}_t^2
        \end{cases}
    \]
    Rescale time via $\tau=\sqrt{\lambda}t$, then
    \begin{equation} \label{eq:rescaled_HFHR}
        \begin{cases}
            d\qbold = (\tpbold - \tilde{\alpha} \nabla \tilde{f}(\qbold)) d\tau + \sqrt{2\tilde{\alpha}\lambda^{-1}} d\bs{B}_\tau^1 \\
            d\tpbold = (-\tilde{\gamma} \tpbold - \nabla \tilde{f}\qbold)) d\tau + \sqrt{2\tilde{\gamma}\lambda^{-1}} d\bs{B}_\tau^2
        \end{cases}
    \end{equation}

    It is easy to see that the rescaled HFHR in Equation  \eqref{eq:rescaled_HFHR} is a tempered-HFHR($\tilde{\alpha}, \tilde{\gamma}, \lambda$).

    \paragraph{Apply Lemma \ref{lemma:tempered_HFHR_convergence_result} to the Tempered HFHR} In Equation  \eqref{eq:rescaled_HFHR}, we have that $\tilde{f}\in C^2(\mathbb{R}^d)$ is convex and $\frac{L}{\lambda}$-smooth. Moreover, $\tilde{\mu} \propto e^{-\tilde{f}(\qbold)}$ satisfies PI with PI constant 1.

    Since
    \begin{align*}
        \tilde{\gamma}^2 \ge \max\{2, \frac{L}{\lambda}\} \iff \gamma^2 \ge \max\{2\lambda, L\} \\
        \tilde{\alpha} \le \tilde{\gamma} - \frac{2}{\tilde{\gamma}} \iff \alpha \le \frac{\gamma}{\lambda} - \frac{2}{\gamma},
    \end{align*}
    we can then apply Lemma \ref{lemma:tempered_HFHR_convergence_result} to the rescaled HFHR in Equation  \eqref{eq:rescaled_HFHR} and obtain the following result
    \begin{equation}\label{eq:rescaled_result}
        \chi^2(\tilde{\rho}_\tau \| \tilde{\pi}) \le e^{- (\frac{1}{2\tilde{\gamma}} + \frac{1}{16}\tilde{\alpha}) \tau} \left\{ \chi^2(\tilde{\rho}_0 \| \tilde{\pi}) + \mathbb{E}_{\tilde{\pi}} \left[ \innerprod{\nabla_{\tilde{\bs{x}}} \frac{\tilde{\rho}_0}{\tilde{\pi}}}{\tilde{S} \nabla_{\tilde{\bs{x}}} \frac{\tilde{\rho}_0}{\tilde{\pi}}} \right] \right\},
    \end{equation}
    where $\tilde{\rho}_\tau(\qbold, \tpbold)$ is the law of $(\qbold_\tau, \tpbold_\tau)$ at time $\tau$, $\tilde{\pi} \propto e^{-\lambda \tilde{H}(\qbold, \tpbold)}$ with $\tilde{H}(\qbold, \tpbold) = \tilde{f}(\qbold) + \frac{1}{2} \|\tpbold\|^2$, $\nabla_{\tilde{\bs{x}}} = (\nablaq, \nabla_{\tpbold})$ and $
    \tilde{S} = \begin{bmatrix} \tilde{a}I & \tilde{b}I \\ \tilde{b}I & \tilde{d}I \end{bmatrix} \text{ with }
    \tilde{a} = (\frac{2}{\tilde{\gamma}} + \tilde{\alpha})\tilde{b}, \tilde{d} = \tilde{\gamma} \tilde{b}, \tilde{b} = \frac{1}{\tilde{\gamma}\lambda}
    $ is a positive definite matrix.

    \paragraph{Substitute Back} We first write all measures with respect to the original position and momentum variables $(\qbold, \pbold)$. Since
    $
        \m{\qbold \\ \pbold} = \m{I & 0 \\ 0 & \sqrt{\lambda}I} \m{\qbold \\ \tpbold} \triangleq P \m{\qbold \\ \tpbold},
    $
    by change of variable formula we have
    \begin{align*}
        \rho_\tau(\qbold, \pbold) =& \tilde{\rho}_\tau(\qbold, \frac{1}{\sqrt{\lambda}}\pbold) \lambda^{-\frac{d}{2}} \\
        \pi (\qbold, \pbold) =& \tilde{\pi}(\qbold, \frac{1}{\sqrt{\lambda}}\pbold) \lambda^{-\frac{d}{2}}
        \propto \exp\{-\lambda \left[\tilde{f}(\qbold) + \frac{1}{2} \|\frac{\pbold}{\sqrt{\lambda}}\|^2 \right] \}
        = \exp\{ f(\qbold) - \frac{1}{2}\|\pbold\|^2 \}
    \end{align*}
    So $\pi(\qbold, \pbold)$ is indeed the joint invariant distribution defined in Equation  \eqref{eq:joint_invariant_distribution}. Therefore,
    \begin{align*}
        \chi^2(\tilde{\rho}_\tau \| \tilde{\pi}) =& \int (\frac{\tilde{\rho}_\tau(\qbold, \tpbold)}{\tilde{\pi}(\qbold, \tpbold)} - 1)^2 \tilde{\pi}(\qbold, \tpbold)d\qbold d\tpbold \\
        =& \int  ( \frac{\tilde{\rho}_\tau(\qbold, \frac{\pbold}{\sqrt{\lambda}})}{\tilde{\pi}(\qbold, \frac{\pbold}{\sqrt{\lambda}}) } - 1)^2 \tilde{\pi}(\qbold, \frac{\pbold}{\sqrt{\lambda}}) \lambda^{-\frac{d}{2}}d\qbold d\pbold  \qquad (\pbold = \sqrt{\lambda}\tpbold)\\
        =& \int (\frac{\rho_\tau(\qbold, \pbold)}{\pi(\qbold, \pbold)} - 1)^2 \pi(\qbold, \pbold) d\qbold d\pbold \\
        =& \chi^2(\rho_\tau \| \pi)
    \end{align*}
    Similar derivation, combined with chain rule, leads to
    \[
        \mathbb{E}_{\tilde{\pi}} \left[ \innerprod{\nabla_{\tilde{\bs{x}}} \frac{\tilde{\rho}_\tau}{\tilde{\pi}}}{\tilde{S} \nabla_{\tilde{\bs{x}}} \frac{\tilde{\rho}_\tau}{\tilde{\pi}} } \right] = \Einfty \left[ \innerprod{P \nablax \frac{\rho_\tau}{\pi}}{\tilde{S} P\nablax \frac{\rho_\tau}{\pi} } \right]
    \]
    where $\nablax = (\nablaq, \nablap)$.

    Next we substitute back the original parameters $\gamma = \tilde{\gamma} \sqrt{\lambda}$ and $\alpha = \frac{\tilde{\alpha}}{\sqrt{\lambda}}$ and rewrite the result in Equation  \eqref{eq:rescaled_result} as
    \begin{equation*}
        \chi^2(\rho_\tau \| \pi) \le e^{- (\frac{\sqrt{\lambda}}{2\gamma} + \frac{\sqrt{\lambda} \alpha}{16} )\tau} \left\{ \chi^2(\rho_0 \| \pi) + \Einfty \left[ \innerprod{\nablax P\frac{\rho_0}{\pi}}{\tilde{S} P \nablax \frac{\rho_0}{\pi}} \right] \right\},
    \end{equation*}

    Now write $S \triangleq P^T \tilde{S} P = \frac{1}{\gamma}\m{(\frac{2}{\gamma} + \alpha)I & I \\ I & \gamma I}$ write $\tau$ as $t$ we obtain
    \begin{equation}
        \chi^2(\rho_t \| \pi) \le e^{- (\frac{\sqrt{\lambda}}{2\gamma} + \frac{\sqrt{\lambda} \alpha}{16} )t} \left\{ \chi^2(\rho_0 \| \pi) + \Einfty \left[ \innerprod{\nablax P\frac{\rho_0}{\pi}}{S P \nablax \frac{\rho_0}{\pi}} \right] \right\}
    \end{equation}
\end{proof}

\subsection{Proof of Themrem \ref{thm:exp_convergence_coupling}}

\begin{proof}
Consider two copies of HFHR that are driven by the same Brownian motion
\[
    \begin{cases}
        d\qbold_t = (\pbold_t - \alpha \nabla f(\qbold_t))dt + \sqrt{2\alpha} d\bs{B}^1_t \\
        d\pbold_t = (-\gamma \pbold_t - \nabla f(\qbold_t))dt + \sqrt{2\gamma} d\bs{B}^2_t
    \end{cases},
    \quad
    \begin{cases}
        d\tilde{\qbold}_t = (\tpbold_t - \alpha \nabla f(\tilde{\qbold}_t))dt + \sqrt{2\alpha} d\bs{B}^1_t \\
        d\tilde{\pbold_t} = (-\gamma \tpbold_t - \nabla f(\tilde{\qbold}_t))dt + \sqrt{2\gamma} d\bs{B}^2_t
    \end{cases},
\]
where we set $(\tilde{\qbold}_0, \tpbold_0) \sim \pi$,  $\pbold_0 = \tpbold_0$ and $\qbold_0$ such that
\[
    W_2^2(\mu_0, \mu) = \mathbb{E} \left[ \| \qbold_0 - \tilde{\qbold}_0 \|^2_2 \right], \quad \qbold_0 \sim \mu_0
\]

Denote $\phipsi = P\deltaqp$ where $P$ is defined in Appendix \ref{app:notation}. By Lemma \ref{lemma:contraction} and the assumption on $\alpha, \gamma$, we have
\[
    \norm{\phipsi}^2 \le e^{-2(\frac{m}{\gamma} + m \alpha) t} \norm{\begin{bmatrix} \bs{\phi}_0 \\ \bs{\psi}_0\end{bmatrix}}^2.
\]
Therefore we obtain
\begin{align*}
    W_2^2(\mu_t, \mu) =& \inf_{(\bs{q}_t, \tilde{\bs{q}}_t) \sim \Pi(\mu_t, \mu)} \mathbb{E}\norm{\bs{q}_t - \tilde{\bs{q}}_t}^2 \\
    \le& \inf_{(\bs{q}_t, \tilde{\bs{q}}_t) \sim \Pi(\mu_t, \mu), (\bs{p}_t, \tilde{\bs{p}}_t) \sim \Pi(\nu_t, \nu)} \mathbb{E} \norm{ \begin{bmatrix} \bs{q}_t - \tilde{\bs{q}}_t \\ \bs{p}_t - \tilde{\bs{p}}_t\end{bmatrix} }^2 \\
    \le& \mathbb{E}\|P^{-1}\|_2^2 \norm{\phipsi}^2 \\
    \le& \mathbb{E} \|P^{-1}\|_2^2 e^{-2(\frac{m}{\gamma} + m\alpha)t} \norm{\begin{bmatrix} \bs{\phi}_0 \\ \bs{\psi}_0\end{bmatrix}}^2 \\
    \le& (\kappa^\prime)^2  e^{-2(\frac{m}{\gamma} + m\alpha)t} \norm{\begin{bmatrix} \bs{q}_0 - \tilde{\bs{q}}_0  \\ \bs{p}_0 - \tilde{\bs{p}}_0\end{bmatrix}}^2 \\
    =& (\kappa^\prime)^2  e^{-2(\frac{m}{\gamma} + m\alpha)t} W_2^2(\mu_0, \mu)
\end{align*}
Taking square root yields the desired result.

\end{proof}

\section{Proof of Theorem \ref{thm:discretization}}
\begin{proof}
Denote $t_k = kh$, the solution of the HFHR dynamics at time $t$ by $\bs{x}_{0, \bs{x}_0}(t)$, the $k$-th iterates of the Strang's splitting method of HFHR dynamics by $\bar{\bs{x}}_{0, \bs{x}_0}(kh)$. Both $\bs{x}_{0, \bs{x}_0}(t)$ and $\bar{\bs{x}}_{0, \bs{x}_0}(kh)$ start from the same initial value $\bs{x}_0$. The linear transformation $P$ defined in Appendix \ref{app:notation}, transforms the solution of HFHR dynamics into $\bs{y}_{0, P\bs{x}_0}(t) = P \bs{x}_{0, \bs{x}_0}(t)$ and the Strang's splitting discretization of HFHR into $\bar{\bs{y}}_{0, P \bs{x}_0}(t) = P\bar{\bs{x}}_{0, \bs{x}_0}(t)$.

For the ease of notation, we write $\bs{y}_{0, \bs{y}_0}(t_k)$ as $\bs{y}_k$ and $\bar{\bs{y}}_{0, \bs{y}_0}(t_k)$ as $\bar{\bs{y}}_k$. We have the following identity
\begin{align*}
    \mathbb{E} \norm{\bs{y}_{k+1} - \bar{\bs{y}}_{k+1}}^2
    =& \mathbb{E} \norm{\bs{y}_{t_k, \bs{y}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h) }^2\\
    =& \mathbb{E} \norm{\bs{y}_{t_k, \bs{y}_k}(h) - \bs{y}_{t_k, \bar{\bs{y}}_k}(h) + \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h) }^2 \\
    =& \underbrace{\mathbb{E}\norm{\bs{y}_{t_k, \bs{y}_k}(h) - \bs{y}_{t_k, \bar{\bs{y}}_k}(h)}^2}_{\circled{1}} + \underbrace{\mathbb{E}\norm{\bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)}^2}_{\circled{2}} \\
    +& 2 \underbrace{\mathbb{E} \left\langle \bs{y}_{t_k, \bs{y}_k}(h) - \bs{y}_{t_k, \bar{\bs{y}}_k}(h) , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle}_{\circled{3}}
\end{align*}
By Lemma \ref{lemma:contraction}, when $0 < h < \frac{1}{2\lambda^\prime}$, term \circled{1} can be upper bounded as
\begin{align*}
    \mathbb{E}\norm{\bs{y}_{t_k, \bs{y}_k}(h) - \bs{y}_{t_k, \bar{\bs{y}}_k}(h)}^2 \le& e^{-2\lambda^\prime h} \mathbb{E}\norm{\bs{y}_k - \bar{\bs{y}}_k}^2 \\
    \le& \left(1 -2\lambda^\prime h + 2(\lambda^\prime)^2h^2 \right)\mathbb{E}\norm{\bs{y}_k - \bar{\bs{y}}_k}^2 \\
    \le& \left(1 - \lambda^\prime h  \right)\mathbb{E}\norm{\bs{y}_k - \bar{\bs{y}}_k}^2
\end{align*}
where the second inequality is due to $e^{-x} \le 1 - x + \frac{x^2}{2}, \forall x > 0$.

For term \circled{2}, we have by Lemma \ref{lemma:local} that
\begin{align*}
    \mathbb{E}\norm{\bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)}^2 \le \sigma_{\text{max}}^2 \, \mathbb{E}\norm{\bs{x}_{t_k, \bar{\bs{x}}_k}(h) - \bar{\bs{x}}_{t_k, \bar{\bs{x}}_k}(h)}^2 \le \sigma_{\text{max}}^2 \, C_2^2 h^3
\end{align*}
where $\sigma_{\text{max}}$ is the largest singular value of matrix $P$.

For term \circled{3}, we have by Lemma \ref{lemma:Z} that
\begin{align*}
    &2\mathbb{E} \left\langle \bs{y}_{t_k, \bs{y}_k}(h) - \bs{y}_{t_k, \bar{\bs{y}}_k}(h) , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle \\
    =& 2\mathbb{E} \left\langle \bs{y}_k - \bar{\bs{y}}_k + \bs{z} , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle \\
    =& \underbrace{2\mathbb{E} \left\langle \bs{y}_k - \bar{\bs{y}}_k , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle}_{\circled{3a}} + \underbrace{2\mathbb{E} \left\langle  \bs{z} , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle}_{\circled{3b}}
\end{align*}

For term \circled{3a}, by the tower property of conditional expectation, we have
\begin{align*}
    2\mathbb{E} \left\langle \bs{y}_k - \bar{\bs{y}}_k , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle
    =& 2\mathbb{E}  \left[\mathbb{E} \left[ \left\langle \bs{y}_k - \bar{\bs{y}}_k , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle \bigg| \mathcal{F}_k \right]\right] \\
    =& 2\mathbb{E} \left\langle \bs{y}_k - \bar{\bs{y}}_k ,  \mathbb{E} \left[   \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h) \bigg| \mathcal{F}_k \right]\right\rangle \\
    \le& 2 \sqrt{\mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k}^2}  \sqrt{\mathbb{E} \norm{ \mathbb{E} \left[   \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h) \bigg| \mathcal{F}_k \right]}^2} \\
    \le& 2 \sqrt{ \mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k}^2}  \sqrt{ \sigma_{\text{max}}^2 \mathbb{E} \norm{ \mathbb{E} \left[   \bs{x}_{t_k, \bar{\bs{x}}_k}(h) - \bar{\bs{x}}_{t_k, \bar{\bs{x}}_k}(h) \bigg| \mathcal{F}_k \right]}^2} \\
    \le& 2 \sqrt{ \mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k}^2}  \sqrt{ \sigma_{\text{max}}^2 C_1^2  h^4} \\
    \le& 2\sigma_{\text{max}} C_1 \sqrt{ \mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k}^2}   h^2.
\end{align*}

For term \circled{3b}, when $0 < h < \frac{1}{4L^{\prime\prime}}$we have by Lemma \ref{lemma:Z} and Lemma \ref{lemma:local}
\begin{align*}
    2\mathbb{E} \left\langle  \bs{z} , \bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)\right\rangle
    \le& 2 \sqrt{\mathbb{E} \norm{\bs{z}}^2} \sqrt{\mathbb{E}\norm{\bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)}^2} \\
    =& 2 \sqrt{\mathbb{E} \norm{\bs{z}}^2} \sqrt{\mathbb{E}\left[\mathbb{E} \left[\norm{\bs{y}_{t_k, \bar{\bs{y}}_k}(h) - \bar{\bs{y}}_{t_k, \bar{\bs{y}}_k}(h)}^2 \bigg| \mathcal{F}_k \right] \right]} \\
    =& 2 \sqrt{\mathbb{E} \norm{\bs{z}}^2} \sqrt{\sigma_{\text{max}}^2 \mathbb{E}\left[\mathbb{E} \left[\norm{\bs{x}_{t_k, \bar{\bs{x}}_k}(h) - \bar{\bs{x}}_{t_k, \bar{\bs{x}}_k}(h)}^2 \bigg| \mathcal{F}_k \right] \right]} \\
    \le& 2 \sigma_{\text{max}} \sqrt{\tilde{C}\mathbb{E}\norm{\bs{y}_k - \bar{\bs{y}}_k}^2 h^2} \sqrt{C_2^2 h^3 } \\
    \le& 2 \sigma_{\text{max}} C_2 \sqrt{\tilde{C}} \sqrt{\mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k}^2} h^\frac{5}{2}
\end{align*}
where $\tilde{C} = 2\left( L^{\prime\prime}\right )^2 = 2(\kappa^\prime)^2 \left(L^\prime\right)^2$ is from Lemma \ref{lemma:Z} and Lemma \ref{lemma:lipschitz}.

Recall both $C_1$ and $C_2$ depend on $\norm{\bs{x}_k}$ and we would like to upper bound this term. To this end,  consider $\tilde{\bs{x}}(t)$, a solution of HFHR dynamics with initial value $\tilde{\bs{x}}_0$ that follows the invariant distribution $\tilde{\bs{x}}_0 \sim \pi$ and realizes $W_2(\pi_0, \pi)$, i.e., $\mathbb{E}\norm{\tilde{\bs{x}}_0 - \bs{x}_0}^2 = W_2^2(\pi_0, \pi)$.

Denote $\tilde{\bs{x}}_k = \tilde{\bs{x}}(kh)$ and $e_k = \left( \mathbb{E} \norm{\bs{y}_k - \bar{\bs{y}}_k }^2\right)^\frac{1}{2}$, we then have
\begin{align*}
    \mathbb{E}\norm{\bar{\bs{x}}_k}^2 =& \mathbb{E} \norm{\bs{x}_k + \bar{\bs{x}}_k - \bs{x}_k}^2 \\
    \le& 2\mathbb{E}\norm{\bs{x}_k}^2 + 2\mathbb{E}\norm{\bar{\bs{x}}_k - \bs{x}_k}^2 \\
    \le& 4\mathbb{E}\norm{\tilde{\bs{x}}_k}^2 + 4\mathbb{E}\norm{\tilde{\bs{x}}_k - \bs{x}_k}^2 + 2\mathbb{E}\norm{\bar{\bs{x}}_k - \bs{x}_k}^2 \\
    =& 4\mathbb{E}\norm{\tilde{\bs{x}}_k}^2 + 4\mathbb{E}\norm{P^{-1} P (\tilde{\bs{x}}_k - \bs{x}_k)}^2 + 2\mathbb{E}\norm{P^{-1} P(\bar{\bs{x}}_k - \bs{x}_k)}^2 \\
    \le& 4\left( \int_{\mathbb{R}^d} \norm{\bs{q}}^2 d\mu + d \right) + \frac{4}{\sigma_{\text{min}}^2}\mathbb{E}\norm{P (\tilde{\bs{x}}_k - \bs{x}_k)}^2 + \frac{2}{\sigma_{\text{min}}^2} \mathbb{E}\norm{\bar{\bs{y}}_k - \bs{y}_k}^2 \\
    \stackrel{(i)}{\le}& 4\left( \int_{\mathbb{R}^d} \norm{\bs{q}}^2 d\mu + d \right) + \frac{4}{\sigma_{\text{min}}^2} e^{-2\lambda^\prime kh}\mathbb{E}\norm{P (\tilde{\bs{x}}_0 - \bs{x}_0)}^2 + \frac{2}{\sigma_{\text{min}}^2} e_k^2 \\
    \le& 4\left( \int_{\mathbb{R}^d} \norm{\bs{q}}^2 d\mu + d \right) + 4\kappa^2 W_2^2(\pi_0, \pi) + \frac{2}{\sigma_{\text{min}}^2} e_k^2 \\
    \triangleq& F e_k^2 + G
\end{align*}
where $(i)$ is due to Lemma \ref{lemma:contraction}. Recall from Lemma \ref{lemma:local}, we have
\begin{align*}
    C_1 \le A_1 \sqrt{\mathbb{E} \norm{\bar{x}_k}^2} + B_1 \le A_1 \sqrt{F} e_k + (A_1 \sqrt{G} + B_1) \triangleq U_1 e_k + V_1\\
    C_2 \le A_2 \sqrt{\mathbb{E} \norm{\bar{x}_k}^2} + B_2 \le A_2 \sqrt{F} e_k + (A_2 \sqrt{G} + B_2) \triangleq U_2 e_k + V_2
\end{align*}
where
\begin{align*}
    A_1 =& (L+G)\max\{ \alpha + 1.25, \gamma + 1\}(1.74 + 0.71\alpha) \\
    B_1 =& (L+G)\max\{ \alpha + 1.25, \gamma +1 \}\left[ 0.5\alpha + (1.26\sqrt{\alpha} + 1.14 \alpha \sqrt{\alpha} + 2.32\sqrt{\gamma})\sqrt{hd}\right] \\
    A_2 =& L\max\{\alpha + 1.25, \gamma + 1\}(1.92 + 2.30\alpha L )\sqrt{h}\\
    B_2 =& L\max\{\alpha + 1.25, \gamma + 1\}(2.60\sqrt{\alpha} + 3.34\sqrt{\gamma}h)\sqrt{d}
\end{align*}
Combine the above and bounds for terms \circled{1}, \circled{2}, \circled{3a} and \circled{3b}, we then obtain
\begin{align*}
    e_{k+1}^2 \le& (1 -\lambda^\prime h) e_k^2 +  \sigma_{\text{max}}^2 \, C_2^2 h^3 + 2\sigma_{\text{max}} C_1 e_k   h^2 + 2\sigma_{\text{max}}C_2 \sqrt{\tilde{C}} e_k h^\frac{5}{2} \\
    \le& (1 - \lambda^\prime h) e_k^2 + \sigma_{\text{max}}^2 2(U_2^2 e_k^2 + V_2^2)h^3 + 2\sigma_{\text{max}} (U_1 e_k + V_1) e_k h^2  + 2\sigma_{\text{max}} (U_2 e_k + V_2) \sqrt{\tilde{C}} e_k h^\frac{5}{2}\\
    =& \left( 1 - \lambda^\prime h + 2\sigma_{\text{max}}^2 U_2^2 h^3 + 2\sigma_{\text{max}} U_1 h^2 + 2\sigma_{\text{max}} U_2 \sqrt{\tilde{C}} h^\frac{5}{2} \right) e_k^2 \\
    &+ \left(2\sigma_{\text{max}} V_1 + 2\sigma_{\text{max}} V_2 \sqrt{\tilde{C} h}  \right) e_k h^2 + 2\sigma_{\text{max}}^2 V_2^2 h^3 \\
    \le& \left( 1 - \lambda^\prime h + 2\sigma_{\text{max}}^2 U_2^2 h^3 + 2\sigma_{\text{max}} U_1 h^2 + 2 \sigma_{\text{max}} U_2 \sqrt{\tilde{C}} h^\frac{5}{2} \right) e_k^2 + \frac{\lambda^\prime}{8}  h e_k^2  \\
    &+ \frac{2\left(2\sigma_{\text{max}} V_1 + 2 \sigma_{\text{max}} V_2 \sqrt{\tilde{C} h}  \right)^2}{\lambda^\prime}h^3 + 2\sigma_{\text{max}}^2 V_2^2 h^3 \\
    =& \left( 1 - \frac{7}{8}\lambda^\prime h + 2\sigma_{\text{max}}^2 U_2^2 h^3 + 2\sigma_{\text{max}} U_1 h^2 + 2 \sigma_{\text{max}} U_2 \sqrt{\tilde{C}} h^\frac{5}{2} \right) e_k^2 \\
    &+ \left( \frac{2\left(2\sigma_{\text{max}} V_1 + 2 \sigma_{\text{max}} V_2 \sqrt{\tilde{C} h}  \right)^2}{\lambda^\prime} + 2\sigma_{\text{max}}^2 V_2^2 \right) h^3 \\
    \stackrel{(i)}{\le}& (1 - \frac{1}{2}\lambda^\prime h) e_k^2 + \left( \frac{2\left(2\sigma_{\text{max}} V_1 + 2 \sigma_{\text{max}} V_2 \sqrt{\tilde{C} h}  \right)^2}{\lambda^\prime} + 2\sigma_{\text{max}}^2 V_2^2 \right) h^3 \\
    \triangleq& (1 - \frac{1}{2}\lambda^\prime h) e_k^2 + K h^3
\end{align*}
where $(i)$ is due to $h < \min\{ h_1, h_2,  h_3 \}$ and
\begin{align*}
    h_1 =& \frac{\sqrt{\lambda^\prime}}{4\sqrt{2} \kappa^\prime L \max\{\alpha + 1.25, \gamma + 1\} (1.92 + 2.30\alpha L)}, \\
    h_2 =& \frac{\lambda^\prime}{16\sqrt{2} \kappa^\prime (L+G) \max\{\alpha+1.25, \gamma + 1\}(1.74 + 0.71\alpha)}, \\
    h_3 =& \frac{\lambda^\prime}{8\kappa^\prime L \max\{\alpha + 1.25, \gamma + 1\} (1.92 + 2.30\alpha L) }.
\end{align*}
Unfolding the above inequality, we arrive at
\begin{align*}
    e_k^2 \le& \left(1 - \frac{\lambda^\prime}{2}h\right)^{k} e_0^2 + \left( 1 + (1 - \frac{\lambda^\prime}{2}h) + \cdots + (1 - \frac{\lambda^\prime}{2}h)^{k-1} \right) K h^3 \\
    \stackrel{(i)}{\le}&K h^3 \sum_{i=0}^\infty \left(1 - \frac{\lambda^\prime}{2}h\right)^i  \\
    =& \frac{2K}{\lambda^\prime} h^2
\end{align*}
where $(i)$ is due to $e_k = 0$. Therefore
\[
    \left(\mathbb{E}\norm{\bs{x}_k - \bar{\bs{x}}_k}^2\right)^\frac{1}{2} = \left(\mathbb{E}\norm{P^{-1}(\bs{y}_k - \bar{\bs{y}}_k)}^2\right)^\frac{1}{2} \le \frac{1}{\sigma_{\text{min}}} e_k \le \frac{1}{\sigma_{\text{min}}} \sqrt{\frac{2K}{\lambda^\prime}}h
\]
Collecting all the constants and we have
\begin{align*}
    \frac{1}{\sigma_{\text{min}}} \sqrt{\frac{2K}{\lambda^\prime}} \le&\frac{8\kappa^\prime}{\lambda^\prime} (L+G) \max\{\alpha+1.25, \gamma+1\} (1.74+0.71\alpha) \left( \sqrt{\int_{\mathbb{R}^d} \norm{\bs{q}}^2 d\mu + d} + \kappa^\prime W_2(\pi_0, \pi)\right) \\
    +& \frac{4\kappa^\prime}{\lambda^\prime} (L+G)\max\{\alpha+1.25,\gamma+1\}\left(0.5\alpha + (1.26\sqrt{\alpha} + 1.14\alpha \sqrt{\alpha} + 2.32\sqrt{\gamma}) \sqrt{d} \right) \\
    +& \frac{8\kappa^\prime}{\sqrt{\lambda^\prime}} \left(\frac{\sqrt{\kappa^\prime L^\prime}}{\sqrt{\lambda^\prime}} + 1 \right) L\max\{\alpha+1.25,\gamma+1\}(1.92+2.30\alpha L) \left( \sqrt{\int_{\mathbb{R}^d} \norm{\bs{q}}^2 d\mu + d} + \kappa^\prime W_2(\pi_0, \pi)\right) \\
    +& \frac{4\kappa^\prime}{\sqrt{\lambda^\prime}} \left(\frac{\sqrt{\kappa^\prime L^\prime}}{\sqrt{\lambda^\prime}} + 1 \right) L\max\{\alpha+1.25,\gamma+1\}(2.60\sqrt{\alpha} + 3.34\sqrt{\gamma})\sqrt{d} \\
    \triangleq C
\end{align*}

It is clear that in terms of the dependence on dimension $d$, we have $C = \mathcal{O}(\sqrt{d})$. In the regime where $\frac{\gamma^2 - L}{\gamma} \ge \frac{m}{\gamma} + m \alpha$, then $\lambda^\prime = \frac{m}{\gamma} + m\alpha$. Recall the definition of $\kappa^\prime$ and there exist $A^\prime, B^\prime > 0$ such that $\kappa^\prime \le A^\prime \sqrt{\alpha} + B^\prime$. It follows that
\[
    C \le \frac{a_1 \alpha^3 + a_2 \alpha^\frac5{2} + a_3 \alpha^2 + a_4 \alpha^\frac{3}{2} + a_5 \alpha + a_6 \alpha^\frac{1}{2} + a_7}{\lambda^\prime} \le \frac{b_1 \alpha^3 + b_2}{\lambda^\prime} =  \frac{b_1\alpha^3 + b_2}{ \frac{m}{\gamma} + m\alpha}
\]
for some positive constants $a_1, a_2, a_3, a_4, a_5, a_6, a_7, b_1, b_2 > 0$ and independent of $\alpha$, in particular, we have $b_1 = \mathcal{O}(\sqrt{d}), b_2= \mathcal{O}(\sqrt{d})$.
\end{proof}

\subsection{Proof of Theorem \ref{thm:wasserstein}}
\begin{proof}
Denote the $k$-th iterate of the Strang's splitting method of HFHR by $\bar{\bs{x}}_k$ with time step $h$, the solution of HFHR dynamics at time $hk$ by $\bs{x}_k$. Both $\bar{\bs{x}}_k$ and $\bs{x}_k$ start from $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0\end{bmatrix}$. Also denote the solution of HFHR dynamics starting from $\tilde{\bs{x}}_0$ at time $kh$ by $\tilde{\bs{x}}_k$ where $\tilde{\bs{x}}_0 = \begin{bmatrix} \tilde{\bs{q}}_0 \\ \bs{p}_0\end{bmatrix}$, $(\tilde{\bs{q}}_0, \tilde{\bs{p}}_0) \sim \pi$ and  $\mathbb{E}\norm{\begin{bmatrix}\bs{q}_0 - \tilde{\bs{q}}_0 \\ \bs{p}_0 - \tilde{\bs{p}}_0\end{bmatrix}}^2 = W_2^2(\pi_0, \pi)$. Since $\pi$ is the invariant distribution of HFHR dynamics, it follows that $\tilde{\bs{x}}_k \sim \pi$.

By Lemma \ref{lemma:contraction} and Theorem \ref{thm:discretization}, we have
\begin{align*}
    \displaystyle W_2^2(\mu_k, \mu) =& \inf_{\xi \in \Pi(\mu_k, \mu)} \mathbb{E}_{(\bs{q}_1, \bs{q}_2) \sim \xi} \norm{\bs{q}_1 - \bs{q}_2}^2 \\
    \le& \inf_{\xi \in \Pi(\pi_k, \pi)} \mathbb{E}_{(\bs{x}_1, \bs{x}_2) \sim \xi} \norm{\bs{x}_1 - \bs{x}_2}^2 \\
    \le& \mathbb{E} \norm{ \bar{\bs{x}}_k - \tilde{\bs{x}}_k}^2 \\
    \le& 2C^2h^2 + 2\mathbb{E} \norm{P^{-1} P (\bs{x}_k - \tilde{\bs{x}}_k)}^2 \\
    \le& 2C^2h^2 + 2\|P^{-1}\|_2^2 \mathbb{E} \norm{ P (\bs{x}_k - \tilde{\bs{x}}_k)}^2 \\
    \le& 2C^2h^2 + 2\|P^{-1}\|_2^2 e^{-2\lambda^\prime kh }\mathbb{E} \norm{ P (\bs{x}_0 - \tilde{\bs{x}}_0)}^2 \\
    \le& 2C^2h^2 + 2(\kappa^\prime)^2 e^{-2\lambda^\prime kh } W_2^2(\pi_0, \pi)
\end{align*}
Take square root on both sides and apply $\sqrt{a^2 + b^2} \le a + b$, we obtain
\[
    W_2(\mu_k, \mu) \le \sqrt{2}Ch + \sqrt{2}\kappa^\prime e^{-\lambda^\prime kh }W_2(\pi_0, \pi).
\]
\end{proof}

\subsection{Proof of Corollary \ref{corollary:iteration complexity}}
\begin{proof}
By Theorem \ref{thm:wasserstein},  we have
\[
    W_2(\mu_k, \mu) \le \sqrt{2}Ch + \sqrt{2}\kappa^\prime e^{-\lambda^\prime kh }W_2(\pi_0, \pi).
\]
Given any target accuracy $\epsilon > 0$, if we run the Strang's splitting method of HFHR with $h^\star = \min\{h_0, \frac{\epsilon}{2\sqrt{2}C}\}$, then after $k^\star = \frac{1}{\lambda^\prime} \max\{\frac{1}{h_0}, \frac{2\sqrt{2}C}{\epsilon}\} \log \frac{2\sqrt{2}\kappa^\prime W_2(\pi_0, \pi)}{\epsilon}$, we have
\[
    W_2(\mu_{k^\star}, \mu) \le \sqrt{2}Ch + \sqrt{2}\kappa^\prime e^{-\lambda^\prime kh }W_2(\mu_0, \mu) \le \frac{\epsilon}{2} + \frac{\epsilon}{2} = \epsilon.
\]
Recall $C = \mathcal{O}(\sqrt{d})$, when high accuracy is needed, e.g. $\epsilon < 2\sqrt{2}Ch_0$, the iteration complexity to reach $\epsilon$-accuracy under 2-Wasserstein distance is $k^\star = \mathcal{O}(\frac{\sqrt{d}}{\epsilon} \log\frac{1}{\epsilon}) = 2\sqrt{2} \frac{C}{\lambda^\prime} \frac{1}{\epsilon} \log \frac{2\sqrt{2} \kappa^\prime W_2(\pi_0, \pi)}{\epsilon} = \tilde{\mathcal{O}}(\frac{\sqrt{d}}{\epsilon})$. Recall from Theorem \ref{thm:discretization}, $C \le \frac{b_1\alpha^3 + b_2}{\frac{m}{\gamma} + m\alpha}$, we have
\[
    \frac{C}{\lambda^\prime} \le \frac{b_1 \alpha^3 + b_2}{ (\frac{m}{\gamma} + m\alpha)^2}
\]
Denote $g(\alpha) = \frac{b_1 \alpha^3 + b_2}{ (\frac{m}{\gamma} + m\alpha)^2}$, simple calculation shows that $g^\prime(0) < 0$ and $\lim_{\alpha \to \infty }g(\alpha) = +\infty$, hence
\[
    \alpha^\prime = \argmin_{\alpha \ge 0} g(\alpha)
\]
exists and is positive.
\end{proof}

\section{Technical/Auxiliary Lemmas and Their Proof}\label{sec:auxiliary results}
\subsection{Poincar\'e's Inequalities for Product Measure}

\begin{lemma}\label{lemma:PI_joint}
    Suppose $\mathcal{X}_1 = \mathcal{X}_2 = \mathbb{R}^d$, and measures $\mu_1 \in \mathcal{P}(\mathcal{X}_1), \mu_2 \in \mathcal{P}(\mathcal{X}_2)$ satisfy Poincar\'e's inequality with constant $\lambda_{\text{PI}}(\mu_1), \lambda_{\text{PI}}(\mu_2)$. Then the product measure $\mu = \mu_1 \otimes \mu_2 \in \mathcal{P}(\mathcal{X}_1 \times \mathcal{X}_2)$ satisfies Poincar\'e's inequality with constant $\PIconstant{\mu} = \min\{\lambda_{\text{PI}}(\mu_1), \lambda_{\text{PI}}(\mu_2)\}$.
\end{lemma}
\begin{proof}
    For any smooth function $f(x_1, x_2)$, denote $g(x_1) = \int fd\mu_2$ and it follows that $\int g d\mu_1 = \int f d\mu$. We have
    \begin{align*}
        \int (f - \int f d\mu)^2 d\mu =& \int (f - g + g - \int f\mu)^2 d\mu \\
        =& \int (f - g)^2 d\mu + 2 \int (f-g)(g - \int fd\mu)d\mu + \int (g - \int f d\mu)^2d\mu \\
        =& \int(\int (f - g)^2 d\mu_2)d\mu_1  + \int (g - \int f d\mu)^2d\mu_1  \\
        \le& \frac{1}{\lambda_{\text{PI}}(\mu_2)} \int (\int \|\nabla_{x_2} f \|^2 d\mu_2) d\mu_1 + \frac{1}{\lambda_{\text{PI}}(\mu_1)} \int \|\nabla_{x_1} g\|^2 d\mu_1 \\
        =& \frac{1}{\lambda_{\text{PI}}(\mu_2)} \int \|\nabla_{x_2} f \|^2 d\mu + \frac{1}{\lambda_{\text{PI}}(\mu_1)} \int \|\int \nabla_{x_1} f d\mu_2\|^2 d\mu_1 \\
        \stackrel{(i)}{\le}& \frac{1}{\lambda_{\text{PI}}(\mu_2)} \int \|\nabla_{x_2} f \|^2 d\mu + \frac{1}{\lambda_{\text{PI}}(\mu_1)} (\int (\int \|\nabla_{x_1} f\|^2 d\mu_1)^\frac{1}{2} d\mu_2)^2 \\
        \stackrel{(ii)}{\le}& \frac{1}{\lambda_{\text{PI}}(\mu_2)} \int \|\nabla_{x_2} f \|^2 d\mu + \frac{1}{\lambda_{\text{PI}}(\mu_1)} \int  \|\nabla_{x_1} f\|^2 d\mu\\
        \le& \frac{1}{\min\{\lambda_{\text{PI}}(\mu_1), \lambda_{\text{PI}}(\mu_2)\}} \int \|\nabla f \|^2 d\mu
    \end{align*}
    where $(i)$ is due to Minkowski's inequality and $(ii)$ is due to Holder's inequality.
\end{proof}

\subsection{Tempered HFHR with Unit PI Constant}
\begin{lemma}\label{lemma:tempered_HFHR_convergence_result}
Under Conditions A\ref{asp:smooth}, A\ref{asp:convex} and suppose $\gamma^2 \ge \max\{2, L\}$ and $\alpha \le \gamma - \frac{2}{\gamma}$. Then the tempered HFHR($\alpha, \gamma, \beta$) in Equation  \eqref{eq:tempered_HFHR} converges to $\pi \propto e^{-\beta H(\qbold, \pbold)}$ where $H(\qbold, \pbold) = f(\qbold) + \frac{1}{2}\|\pbold\|^2$. Moreover, if the joint invariant distribution $\pi$ satisfies PI with PI constant $\PIconstant{\pi} = \beta$, we have the following exponential convergence
\[
    \chi^2(\rho_t \| \pi) \le e^{- (\frac{1}{2\gamma} + \frac{1}{16}\alpha) t} \left\{ \chi^2(\rho_0 \| \pi) + \Einfty \left[ \innerprod{\nablax \frac{\rho_0}{\pi}}{S \nablax \frac{\rho_0}{\pi}} \right] \right\},
\]
where $\chi^2(\mu, \nu) = \int (\frac{d\mu}{d\nu} - 1)^2 d\nu$ and $\rho_t$ is the joint law of $(\qbold_t, \pbold_t)$ of tempered HFHR($\alpha, \gamma, \beta$) at time $t$,  $\nablax = (\nablaq, \nablap)$ and $S \in \mathbb{R}^{2d \times 2d}$ is a symmetric matrix, more specifically,
\[
    S = \begin{bmatrix} aI & bI \\ bI & dI \end{bmatrix} \text{ with }
    a = (\frac{2}{\gamma} + \alpha)b, d = \gamma b, b = \frac{1}{\gamma\beta}.
\]
\end{lemma}

\begin{proof}
    Denote the eigenvalues of $\nabla^2 f$ by $\eta_i, i=1,2,\cdots,d$. By convexity assumption on $f$ and $L$-smoothness assumption on $\nabla f$, we have $0 \le \eta_i \le L, i=1,2,\cdots,d$.

    By assumption, $\kappa \propto e^{-f(\qbold)}$ satisfies PI with PI constant 1, it is easy to see $\mu \propto e^{-\beta f(\qbold)}$ satisfies PI with PI constant $\beta$.
    It is well known that standard Gaussian measure satisfies PI with constant 1, so $\nu \propto e^{-\frac{\beta}{2}\|\pbold\|^2}$ satisfies PI with PI constant $\beta$. By Lemma \ref{lemma:PI_joint}, we know $\pi$ satisfies PI with constant $\lambda_{\text{PI}}(\pi) = \min\{\PIconstant{\mu}, \PIconstant{\nu}\} = \beta$.

    Consider the following Lyapunov function, $\chi^2$ divergence augmented by the cross term
    \begin{equation*}
        \mathcal{L}(\rho_t) = \chi^2(\rho_t \| \pi) + \crossterm
    \end{equation*}
    where $\crossterm$ is defined in Equation  \eqref{eq:cross_term} with
    \begin{equation}
        S = \begin{bmatrix} aI & bI \\ bI & dI \end{bmatrix} \text{ with }
        a = (\frac{2}{\gamma} + \alpha)b, d = \gamma b, b = \frac{1}{\gamma\beta}
    \end{equation}

    By direct computation and Lemma \ref{lemma:cross_term}, we have
    \begin{equation}\label{eq:time_derivative_augmented_lyapunov_function}
        \frac{d}{dt} \mathcal{L}(\rho_t) \le - \Einfty
        \left[
            \innerproduct{\nablax \frac{\rho_t}{\pi}}{(\beta^{-1}D + M_{\text{cross}}) \nablax \frac{\rho_t}{\pi}}
        \right]
    \end{equation}
    where $D \triangleq \begin{bmatrix} 2\alpha I & 0 \\ 0 & 2\gamma I  \end{bmatrix}$ and $M_{\text{cross}}$ is the matrix from Equation  \eqref{eq:M_cross} in Lemma \ref{lemma:cross_term}.

    Denote
    \begin{align*}
        A \triangleq& \beta^{-1}D + M_{\text{cross}} \\
        =& \beta^{-1}\begin{bmatrix} 2\alpha I & 0 \\0 & 2\gamma I \end{bmatrix} +
        \begin{bmatrix}  2bI + 2a\alpha \nabla^2 f(\qbold) & -(a-b\alpha)\nabla^2 f(\qbold) + b\gamma I + dI \\
        -(a-b\alpha)\nabla^2 f(\qbold) + b\gamma I + dI & -2b\nabla^2 f(\qbold) + 2d\gamma I + 2\gamma I
        \end{bmatrix} \\
        =& 2\beta^{-1}\begin{bmatrix}\alpha I & 0\\0 & \gamma I\end{bmatrix}
        + 2a\alpha \begin{bmatrix} \nabla^2 f(\qbold) & 0 \\ 0 & 0 \end{bmatrix}
        + \frac{2}{\gamma\beta} \begin{bmatrix} I & \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) \\  \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) & \gamma^2 I - \nabla^2 f(\qbold)\end{bmatrix}
    \end{align*}

    Denote
    \[
    E =\frac{1}{\PIconstant{\pi}} I + S = \frac{1}{\beta} I + S = \frac{1}{\beta} I
    + \frac{1}{\gamma \beta}\m{(\frac{2}{\gamma} + \alpha) I & I \\ I & \gamma I}.
    \]
    The rest of the proof is dedicated to matrix analysis of $A$ and $E$. Denote $P \succeq Q$ if $P - Q$ is a positive semi-definite matrix and we will frequently use the following property of block matrix
    \begin{equation}\label{eq:block_matrix_det}
        \det\left( \m{A_{11} & A_{12} \\ A_{21} & A_{22}} \right) = \det \left( A_{11}A_{22} -A_{12}A_{21} \right) \mbox{ if } A_{21}A_{22} = A_{22}A_{21}.
    \end{equation}
    First since $\alpha \le \gamma - \frac{2}{\gamma}$, we have
    \begin{equation}\label{eq:small_result1}
        E = \frac{1}{\beta} I + \frac{1}{\gamma\beta}\begin{bmatrix} (\alpha +\frac{2}{\gamma})I & I \\I & \gamma I \end{bmatrix} \preceq \frac{1}{\beta} I +  \frac{1}{\gamma\beta}\begin{bmatrix} \gamma I & I \\I & \gamma I \end{bmatrix}
        = \frac{1}{\beta}\begin{bmatrix} 2I & \frac{1}{\gamma}I \\ \frac{1}{\gamma}I & 2I \end{bmatrix}
    \end{equation}
    Since $2a\alpha \m{\nabla^2 f(\qbold) & 0 \\ 0 & 0} \succeq 0$, we have
    \begin{equation}\label{eq:small_result2}
        A \succeq 2\beta^{-1}\m{\alpha I & 0\\0 & \gamma I}
        + \frac{2}{\gamma\beta} \m{ I & \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) \\  \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) & \gamma^2 I - \nabla^2 f(\qbold)}
    \end{equation}
    Now we consider the following difference
    \begin{align*}
        & \, A - (\frac{1}{2\gamma} + \z \alpha) E \\
        \succeq& \underbrace{\frac{2}{\gamma\beta} \m{I & \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) \\ \gamma I - \frac{1}{\gamma}\nabla^2 f(\qbold) & 2\gamma^2 I - \nabla^2 f(\qbold)}
        - \frac{1}{\gamma\beta}\m{I & \frac{1}{2\gamma}I \\ \frac{1}{2\gamma}I & I}}_{B} \\
        +& 2\beta^{-1} \m{\alpha I & 0 \\ 0 & 0}  - \frac{\z \alpha}{2\beta}\m{I & \frac{1}{2\gamma}I \\ \frac{1}{2\gamma}I & I}
    \end{align*}
    where the matrix inequality is due to assumption $\gamma^2 \ge L$,  Equation  \eqref{eq:small_result1} and \eqref{eq:small_result2}.

    For $B$, we have
    \[
        B = \frac{1}{\gamma\beta} \m{I & 2\gamma I - \frac{2}{\gamma}\nabla^2 f(\qbold) - \frac{1}{2\gamma}I \\ 2\gamma I - \frac{2}{\gamma}\nabla^2 f(\qbold) - \frac{1}{2\gamma}I & 4\gamma^2 I -2\nabla^2 f(\qbold) - I}
    \]
    Using Equation  \eqref{eq:block_matrix_det} and diagonalization of $\nabla^2 f(\qbold)$, we know the eigenvalues of $B$ are the collections of eigenvalues of the following $2 \times 2$ matrices
    \[
        B_i = \frac{1}{\gamma\beta} \m{1 & 2\gamma I - \frac{2\eta_i}{\gamma} - \frac{1}{2\gamma}\\ 2\gamma - \frac{2\eta_i}{\gamma} - \frac{1}{2\gamma}I & 4\gamma^2 -2\eta_i - 1}
    \]
    Notice
    \[
        \det(B_i) =  \frac{1}{\gamma^2\beta^2} (1 - \frac{1}{4\gamma^2} + \frac{6\eta_i \gamma^2 - 2\eta_i - 4\eta_i^2}{\gamma^2}) \ge \frac{1}{\gamma^2\beta^2} (1 - \frac{1}{4\gamma^2}) > 0
    \]
    hence each $B_i$ is positive definite and the smaller eigenvalue of each $B_i$ is
    \begin{align*}
        \lambda^-(B_i)
        =& \frac{1}{2\gamma\beta}\left[
        4\gamma^2 - 2\eta_i - \sqrt{1 - \frac{1}{4\gamma^2} + \frac{6\eta_i \gamma^2 - 2\eta_i - 4\eta_i^2}{\gamma^2}}
        \right] \\
        =& \frac{1}{2\gamma\beta}\left[
        4\gamma^2 - 2\eta_i - \sqrt{(1 - \frac{1}{4\gamma^2}) + \frac{(2\eta_i \gamma^2 - 2\eta_i) + (4\eta_i \gamma^2 -4\eta_i^2)}{\gamma^2}}
        \right]\\
        \stackrel{(i)}{\ge}& \frac{1}{2\gamma\beta}\left[
        2\gamma^2 - \sqrt{1 + 2\gamma^2 - 2 + 4L}
        \right] \\
        \stackrel{(ii)}{\ge}& \frac{1}{2\gamma\beta} (2\gamma^2 - \sqrt{6}\gamma)\\
        \ge& \frac{\gamma}{8\beta}
    \end{align*}
    where $(i), (ii)$ follow by $0\le\eta_i\le L$ and the assumption $\gamma^2 \ge \max\{2, L\}$.

    Therefore, the smallest eigenvalue of $B$ is also lower bounded
    \[
        \lambda^{\min}(B) = \min_{i=1,2,\cdots, d} \lambda^-(B_i) \ge \frac{\gamma}{8\beta}.
    \]
    and this is equivalent to $B \succeq \frac{8\gamma}{\beta} I$.

    We now obtain
    \begin{align*}
         A - (\frac{1}{2\gamma} + \z \alpha)E
        \succeq& 2\beta^{-1} \m{\alpha I & 0\\0&0} + \frac{\gamma}{8\beta} I - \frac{\z \alpha}{2\beta} \m{I&\frac{1}{2\gamma}I \\\frac{1}{2\gamma}I & I}\\
        =& \underbrace{\frac{1}{8\beta} \m{(4\alpha - \frac{1}{4}\alpha + \gamma)I & -\frac{\alpha}{8\gamma}I \\
        -\frac{\alpha}{8\gamma}I & (\gamma - \frac{\alpha}{4})I }}_{F}
    \end{align*}
    By Equation  \eqref{eq:block_matrix_det}, it is easy to see that the eigenvalues of $F$
    are identical to the eigenvalues (ignoring multiplicity) of the following $2 \times 2$ matrix
    $
        \tilde{F} = \frac{1}{8\beta} \m{4\alpha - \frac{1}{4}\alpha + \gamma & -\frac{\alpha}{8\gamma} \\
        -\frac{\alpha}{8\gamma} & \gamma - \frac{\alpha}{4} }
    $ and we have
    \[
        \det(\tilde{F}) = \frac{1}{(8\beta)^2} \left( -\frac{15\alpha^2}{16} - \frac{\alpha^2}{64\gamma^2} + \frac{7\alpha\gamma}{2} + \gamma^2 \right)
        \stackrel{(i)}{\ge} \frac{1}{(8\beta)^2} \left( \frac{1}{16}\gamma^2 - \frac{1}{64} + \frac{7\alpha\gamma}{2}  \right)
        \stackrel{(ii)}{>} 0
    \]
    where $(i)$ and $(ii)$ follow by the assumption $\alpha \le \gamma - \frac{2}{\gamma} < \gamma$ and $\gamma^2 > 2$. Therefore $\tilde{F}$ is a positive definite matrix and all of its eigenvalues are positive, and hence the eigenvalues of $F$ are also positive, equivalently $F \succeq 0$.

    We now have established the relation
    $
        A - (\frac{1}{2\gamma} + \frac{1}{16}\alpha)E \succeq 0
    $
    and return to the time derivative in Equation  \eqref{eq:time_derivative_augmented_lyapunov_function}
    \begin{align*}
        \frac{d}{dt}\mathcal{L}(\rho_t) \le& -\Einfty \left[ \innerproduct{\nablax \frac{\rho_t}{\pi}}{A \nablax \frac{\rho_t}{\pi}}  \right] \\
        \le& -(\frac{1}{2\gamma} + \frac{1}{16}\alpha) \Einfty \left[ \innerproduct{\nablax \frac{\rho_t}{\pi}}{E \nablax \frac{\rho_t}{\pi}}  \right] \\
        =& -(\frac{1}{2\gamma} + \frac{1}{16}\alpha) \left[ \frac{1}{\PIconstant{\pi}} \Einfty \left[ \|\nablax \frac{\rho_t}{\pi}\|^2  \right] + \crossterm \right] \\
        \stackrel{(i)}{\le}& -(\frac{1}{2\gamma} + \frac{1}{16}\alpha) \left[ \Einfty [ (\frac{\rho_t}{\pi} - 1)^2 ] + \crossterm \right]\\
        =& -(\frac{1}{2\gamma} + \frac{1}{16}\alpha) \mathcal{L}(\rho_t)
    \end{align*}
    where $(i)$ is due to Poincar\'e's inequality.

    By Gronwall's inequality, $\mathcal{L}$ has exponential decay
    $
        \mathcal{L}(\rho_t) \le e^{- (\frac{1}{2\gamma} + \frac{1}{16}\alpha) t} \mathcal{L}(\rho_0)
    $
    and since $S$ is positive definite, we further have
    \[
        \chi^2(\rho_t \| \pi) \le \mathcal{L}(\rho_t) \le e^{- (\frac{1}{2\gamma} + \frac{1}{16}\alpha) t} \mathcal{L}(\rho_0),
    \]
\end{proof}

\subsection{Time Derivative of $M_{\text{cross}}$}
\begin{lemma}\label{lemma:cross_term}
\begin{equation*}
    \frac{d}{dt} \crossterm \le - \beta^{-1} \int \innerproduct{\nablax \frac{\rho_t}{\pi}}{ M_{\text{cross}} \nablax \frac{\rho_t}{\pi}} d\pi
\end{equation*}
where $\rho_t$ is the law of $(\qbold_t, \pbold_t)$ of tempered-HFHR($\alpha, \gamma, \beta$) and
\begin{equation}\label{eq:M_cross}
    M_{\text{cross}} =
        \begin{bmatrix}
            2a\alpha \nabla^2 f(\bs{\qbold}) + 2b I  &
            -a \nabla^2 f(\qbold) + b\alpha \nabla^2 f(\qbold) +b\gamma I + d I \\
            -a \nabla^2 f(\qbold) + b\alpha \nabla^2 f(\qbold) +b\gamma I + d I &
            -2b \nabla^2 f(\qbold) + 2d\gamma I
        \end{bmatrix}
\end{equation}
\end{lemma}

\begin{proof}
    For better readability, we collect some notations used in the proof here
    \[
        \mu(\qbold) \propto e^{-\beta f(\qbold)}, \quad \nu(\pbold) \propto e^{-\frac{\beta}{2}\|\pbold\|^2}, \quad \pi(\qbold, \pbold) \propto e^{-\beta H(\qbold, \pbold)}
    \]
    where $H(\qbold, \pbold) = f(\qbold) + \frac{1}{2}\|\pbold\|^2$ and write tempered HFHR again for reference
    \begin{equation*}
        \begin{cases}
            d\qbold = \left( \pbold - \alpha \nabla f(\qbold) \right)dt + \sqrt{2\alpha \beta^{-1}} d\bs{B}^1_t \\
            d\pbold = \left( -\gamma \pbold - \nabla f(\qbold) \right)dt + \sqrt{2\gamma \beta^{-1}} d\bs{B}^2_t
        \end{cases}.
    \end{equation*}
    The Fokker-Planck equation of tempered HFHR in Equation  \eqref{eq:tempered_HFHR} is given by
    \begin{equation}
         \partial_t \rho_t + \nabla \cdot (\rho_t\bs{J}) = 0, \quad \text{where }         \bs{J} = \left (
            \begin{bmatrix}
            \pbold - \alpha \nabla f(\bs{q}) - \alpha \beta^{-1} \nabla_{\bs{q}} \log \rho_t \\
            -\gamma \bs{p} - \nabla f(\bs{q}) - \gamma \beta^{-1}  \nabla_{\bs{p}} \log \rho_t
            \end{bmatrix}
                \right )
    \end{equation}
    Since
    $
        \nabla \cdot \left(\rho_t \begin{bmatrix} -\nabla_{\bs{p}} \log \rho_t \\ \nabla_{\bs{q}} \log \rho_t \end{bmatrix}\right) = 0,
    $
    we can then further simplify $\bs{J}$ to
    \begin{equation}\label{eq:vector_field}
        \bs{J} = -\beta^{-1} \frac{\pi}{\rho_t} A
            \begin{bmatrix}
                \nabla_{\qbold} \frac{\rho_t}{\pi} \\
                \nabla_{\pbold}  \frac{\rho_t}{\pi}
            \end{bmatrix},
    \end{equation}
    where $A = \A$.

    The functional derivative w.r.t. $\rho_t$ is
    \[
        \frac{\delta \crossterm}{\delta \rho_t} = 2 (\nablax)^\star (S \nablax \frac{\rho_t}{\pi})
    \]
    where $(\nablax)^\star$ is the adjoint operator with respect to $\Einfty[\innerproduct{\cdot}{\cdot}]$ and
    \[
        (\nablax)^\star = -\nablax^T - \nablax^T \log \pi = ((\nablaq)^\star, (\nablap)^\star) = (-\nablaq^T + \beta (\bs{\nabla} f(\qbold))^T, -\nablap^T +  \beta \pbold^T).
    \]
    The time derivative of the $\crossterm$ is
    \begin{align}
        \frac{d}{dt} \crossterm(\rho_t) =& \int \frac{\delta \crossterm}{\delta \rho_t} \partial_t \rho_t d\bs{x} \nonumber\\
        =& -\int \frac{\delta \crossterm}{\delta \rho_t} \nablax \cdot (\rho_t \bs{J}) d\bs{x} \nonumber\\
        =& \int \innerprod{\nablax \frac{\delta \crossterm}{\delta \rho_t}}{\bs{J}}\rho_t d\bs{x} \nonumber\\
        =& -2 \beta^{-1} \int \innerprod{\nablax (\nablax)^\star(S\nablax \frac{\rho_t}{\pi})}{\A \nablax \frac{\rho_t}{\pi}} d\pi \nonumber\\
        =& -2 \beta^{-1} \int \innerprod{\nablax (\nablax)^\star(S\nablax h)}{\A \nablax h} d\pi \qquad (h \triangleq \frac{\rho_t}{\pi}) \label{term:2}
    \end{align}
    For the term in Equation  \eqref{term:2}, we have
    \begin{align}
        & -2 \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablax)^\star S \nablax h)}{\A\nablax h}] \nonumber \\
        =& -2a \beta^{-1} \Einfty [\innerproduct{\nablax (\nablaq)^\star\nablaq h}{\A \nablax h}] \label{term:2a} \\
        & - 2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h) }{\A \nablax h}] \label{term:2b} \\
        & -2d \beta^{-1} \Einfty[\innerproduct{\nablax (\nablap)^\star \nablap h}{\A \nablax h}] \label{term:2c}
    \end{align}
    For the cross term in Equation  \eqref{term:2b}, we have
    \begin{align}
        & - 2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h) }{\A \nablax h}] \nonumber \\
        =& -2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\D \nablax h}] \nonumber \\
        & -2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\Q \nablax h}] \nonumber \\
        =& -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablaq ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\nablaq h}] \label{term:2b1}\\
        & -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablap ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\nablap h}] \label{term:2b2}\\
        & -2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\Q \nablax h}] \label{term:2b3}
    \end{align}
    For the term in Equation  \eqref{term:2b1}, we have
    \begin{align*}
        & -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablaq ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\nablaq h}] \nonumber \\
        =& -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablaq (\nablaq)^\star \nablap h }{\nablaq h} + \innerproduct{\nablaq (\nablap)^\star \nablaq h}{\nablaq h}] \nonumber \\
        =& -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablap h }{ \nablaq (\nablaq)^\star \nablaq h} + \innerproduct{(\nablaq)^\star \nablaq  \nablaq h}{\nablap h}] \nonumber \quad (\nablap \mbox{ commutes with } \nablaq, (\nablaq)^\star) \\
        =& -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablap h }{ ((\nablaq)^\star \nablaq + \nablaq (\nablaq)^\star) \nablaq h}] \nonumber \\
        =& -2b\alpha \beta^{-1} \Einfty[\innerproduct{\nablap h }{ (2(\nablaq)^\star \nablaq + [\nablaq, (\nablaq)^\star]) \nablaq h}] \nonumber \\
        =& -4b\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablaq h}]
        - 2b\alpha \Einfty[\innerproduct{\nablap h}{\nabla^2 f(\qbold) \nablaq h}]
    \end{align*}
    where we make use of the commutator $[\nablaq, (\nablaq)^\star]$ of $\nablaq$ and $(\nablaq)^\star$
    \[
        [\nablaq, (\nablaq)^\star] = \nablaq(\nablaq)^\star - (\nablaq)^\star \nablaq = \beta \nabla^2 f(\qbold) + \nablaq^T \nablaq - \nablaq \nablaq^T.
    \]
    For the term in Equation  \eqref{term:2b2}, we have
    \begin{align*}
        & -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablap ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\nablap h}] \nonumber \\
        =& -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablap (\nablaq)^\star \nablap h }{\nablap h} + \innerproduct{\nablap (\nablap)^\star \nablaq h}{\nablap h}] \nonumber \\
        =& -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablaq h }{ (\nablap)^\star \nablap  \nablap h} + \innerproduct{\nablaq h}{ \nablap (\nablap)^\star \nablap h}] \nonumber \quad (\nablaq \mbox{ commutes with } \nablap, (\nablap)^\star) \\
        =& -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablaq h }{ ((\nablap)^\star \nablap + \nablap (\nablap)^\star) \nablap h}] \nonumber \\
        =& -2b\gamma \beta^{-1} \Einfty[\innerproduct{\nablaq h }{ (2(\nablap)^\star \nablap + [\nablap, (\nablap)^\star]) \nablap h}] \nonumber \\
        =& -4b\gamma \beta^{-1}\Einfty[\innerproductF{\nablap \nablaq h}{\nablap \nablap h}]
        - 2b\gamma \Einfty[\innerproduct{\nablap h}{ \nablaq h}]
    \end{align*}
    where we make use of the commutator $[\nablap, (\nablap)^\star]$ of $\nablap$ and $(\nablap)^\star$
    \[
        [\nablap, (\nablap)^\star] = \nablap(\nablap)^\star - (\nablap)^\star \nablap = \beta I + \nablap^T \nablap - \nablap \nablap^T.
    \]
    For the term in Equation  \eqref{term:2b3}, we have
    \begin{align*}
        & -2b \beta^{-1} \Einfty[\innerproduct{\nablax ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)}{\Q \nablax h}] \\
        =& -2b \beta^{-1} \Einfty[ ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h) \cdot (\nablax h)^\star(\Q \nablax h)] \\
        =& -2b \beta^{-1} \Einfty[ ((\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h) \cdot (-(\nablaq)^\star \nablap h + (\nablap)^\star \nablaq h)] \\
        =& -2b \beta^{-1} \Einfty[ ((\nablap)^\star \nablaq h)^2 - ((\nablaq)^\star \nablap h)^2] \\
        =& -2b \beta^{-1} \Einfty[ \innerproduct{\nablaq h}{\nablap (\nablap)^\star \nablaq h} - \innerproduct{\nablap h}{\nablaq (\nablaq)^\star \nablap h}] \\
        =& -2b \beta^{-1} \Einfty[ \innerproduct{\nablaq h}{\nablap (\nablap)^\star \nablaq h} - \innerproduct{\nablaq h}{(\nablap)^\star \nablap \nablaq h} + \innerproduct{\nablaq h}{(\nablap)^\star \nablap \nablaq h}  - \innerproduct{\nablap h}{\nablaq (\nablaq)^\star \nablap h}] \\
        =& -2b \beta^{-1} \Einfty[ \innerproduct{\nablaq h}{ [\nablap, (\nablap)^\star] \nablaq h} + \innerproduct{\nablaq h}{(\nablap)^\star \nablap \nablaq h}  - \innerproduct{\nablap h}{\nablaq (\nablaq)^\star \nablap h}] \\
        =& -2b \beta^{-1} \Einfty[ \innerproduct{\nablaq h}{ [\nablap, (\nablap)^\star] \nablaq h} + \innerproduct{\nablap h}{(\nablaq)^\star \nablaq \nablap h}  - \innerproduct{\nablap h}{\nablaq (\nablaq)^\star \nablap h}] \\
        =& -2b \beta^{-1} \Einfty[ \innerproduct{\nablaq h}{ [\nablap, (\nablap)^\star] \nablaq h} - \innerproduct{\nablap h}{[\nablaq, (\nablaq)^\star] \nablap h}] \\
        =& -2b \Einfty[ \innerproduct{\nablaq h}{  \nablaq h} - \innerproduct{\nablap h}{\nabla^2 f(\qbold) \nablap h}]
    \end{align*}
    For the quadratic term in Equation  \eqref{term:2a}, we have
    \begin{align*}
        & -2a \beta^{-1} \Einfty [\innerproduct{\nablax (\nablaq)^\star\nablaq h}{\A \nablax h}] \\
        =& -2a \beta^{-1} \Einfty [ (\nablaq)^\star\nablaq h \cdot (\nablax)^\star\begin{bmatrix} \alpha \nablaq h - \nablap h \\ \nablaq h + \gamma \nablap h \end{bmatrix}] \\
        =& -2a \beta^{-1} \Einfty [ (\nablaq)^\star\nablaq h \cdot (\alpha (\nablaq)^\star\nablaq h - (\nablaq)^\star\nablap h + (\nablap)^\star \nablaq h + \gamma(\nablap)^\star \nablap h)] \\
        =& -2a\alpha \beta^{-1} \Einfty[\innerproduct{\nablaq h}{\nablaq (\nablaq)^\star \nablaq h}] - 2a \beta^{-1} \Einfty[-(\nablaq)^\star\nablaq h (\nablaq)^\star \nablap h + (\nablaq)^\star \nablaq h (\nablap)^\star \nablaq h] \\
        & - 2a\gamma \beta^{-1} \Einfty[(\nablaq)^\star \nablaq (\nablap)^\star \nablap h] \\
        =& -2a\alpha \beta^{-1} \Einfty[\innerproduct{\nablaq h}{((\nablaq)^\star \nablaq + \beta \nabla^2 f(\qbold)) \nablaq h}] + 2a \beta^{-1} \Einfty[\innerproduct{\nablap h}{[\nablaq, (\nablaq)^\star] \nablaq h}] \\
        &- 2a \beta^{-1} \Einfty[\innerproductF{\nablap \nablaq h}{\nablaq \nablaq h}] \\
        =& -2a\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablaq h}{\nablaq \nablaq h}] -2a\alpha \Einfty[\innerproduct{\nablaq h}{\nabla^2 f(\qbold) \nablaq h}] + 2a  \Einfty[\innerproduct{\nablap h}{\nabla^2 f(\qbold) \nablaq h}] \\
        & - 2a\gamma \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablap h}]
    \end{align*}
    Similarly, for the term in Equation  \eqref{term:2c}, we have
    \begin{align*}
        & -2d \beta^{-1} \Einfty [\innerproduct{\nablax (\nablap)^\star\nablap h}{\A \nablax h}] \\
        =& -2d \beta^{-1} \Einfty [ (\nablap)^\star\nablap h \cdot (\nablax)^\star\begin{bmatrix} \alpha \nablaq h - \nablap h \\ \nablaq h + \gamma \nablap h \end{bmatrix}] \\
        =& -2d \beta^{-1} \Einfty [ (\nablap)^\star\nablap h \cdot (\alpha (\nablaq)^\star\nablaq h - (\nablaq)^\star\nablap h + (\nablap)^\star \nablaq h + \gamma(\nablap)^\star \nablap h)] \\
        =& -2d\alpha \beta^{-1} \Einfty[\innerproduct{\nablap h}{\nablap (\nablaq)^\star \nablaq h}] - 2d \beta^{-1} \Einfty[-(\nablap)^\star\nablap h (\nablaq)^\star \nablap h + (\nablap)^\star \nablap h (\nablap)^\star \nablaq h] \\
        & - 2d\gamma \beta^{-1} \Einfty[\innerproduct{\nablap}{\nablap(\nablap)^\star \nablap h}] \\
        =& -2d\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablap h}] - 2d \beta^{-1} \Einfty[\innerproduct{\nablaq h}{[\nablap, (\nablap)^\star]\nablap h}] \\
        & - 2d\gamma \beta^{-1} \Einfty[\innerproduct{\nablap}{((\nablap)^\star \nablap + \beta I) \nablap h}] \\
        =& -2d\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablap h}] - 2d \Einfty[\innerproduct{\nablaq h}{\nablap h}] \\
        & - 2d\gamma \beta^{-1} \Einfty[\innerproductF{\nablap \nablap h}{\nablap \nablap h}] -2d\gamma \Einfty[\innerproduct{\nablap h}{ \nablap h}]
    \end{align*}
    We now collect all terms in regular Euclidean inner product, i.e., $\innerproduct{}{}$, we have
    \[
        \begin{cases}
        \innerproduct{\nablaq h}{\nablaq h}: - 2a\alpha \nabla^2 f(\qbold) - 2b \\
        \innerproduct{\nablaq h}{\nablap h}: 2a\nabla^2  f(\qbold) - 2b\alpha  \nabla^2 f(\qbold) -2b\gamma  - 2d  \\
        \innerproduct{\nablap h}{\nablap h}: 2b\nabla^2 f(\qbold) - 2d\gamma
        \end{cases}
    \]
    Therefore, if we denote
    \begin{equation}
        M_{\text{cross}} =
        \begin{bmatrix}
            2a\alpha \nabla^2 f(\qbold) + 2b I  & -a \nabla^2 f(\qbold) + b\alpha \nabla^2 f(\qbold) +b\gamma I + d I \\
           -a \nabla^2 f(\qbold) + b\alpha \nabla^2 f(\qbold) +b\gamma I + d I & - 2b \nabla^2 f(\qbold) + 2d\gamma I
        \end{bmatrix}
    \end{equation}
    then the component containing regular Euclidean inner product can be written in a compact form
    \[
        - \Einfty[\innerproduct{\nablax h}{ M_{\text{cross}} \nablax h}].
    \]
    Next, we collect all terms in Frobenius inner product, i.e., $\innerproductF{}{}$, we have
    \begin{align*}
        & \left\{
        -2a\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablaq h}{\nablaq \nablaq h}]
        -4b\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablaq h}]
        -2d\alpha \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablap h}]
        \right\} \\
        +& \left\{
        - 2a\gamma \beta^{-1} \Einfty[\innerproductF{\nablaq \nablap h}{\nablaq \nablap h}]
        -4b\gamma  \beta^{-1} \Einfty[\innerproductF{\nablap \nablaq h}{\nablap \nablap h}]
        - 2d\gamma \beta^{-1} \Einfty[\innerproductF{\nablap \nablap h}{\nablap \nablap h}]
        \right\} \\
        =& -2\alpha \beta^{-1} \Einfty[\innerproductF{\nablax \nablaq h}{S \nablax \nablaq h}] -2\gamma \beta^{-1} \Einfty[\innerproductF{\nablax \nablap h}{S \nablax \nablap h}]
    \end{align*}
    Now we sum up all terms and obtain
    \begin{align*}
        \frac{d}{dt} \crossterm(\rho_t) =& - \Einfty [\innerproduct{\nablax \frac{\rho_t}{\pi}}{ M_{\text{cross}} \nablax  \frac{\rho_t}{\pi}}] \\
        & -2\alpha \beta^{-1} \Einfty[\innerproductF{\nablax \nablaq h}{S \nablax \nablaq h}] -2\gamma \beta^{-1} \Einfty[\innerproductF{\nablax \nablap h}{S \nablax \nablap h}]\\
        \le& - \Einfty [\innerproduct{\nablax \frac{\rho_t}{\pi}}{ M_{\text{cross}} \nablax \frac{\rho_t}{\pi}}]
    \end{align*}
\end{proof}

\subsection{Dependence of error of SDE on initial values}
\begin{lemma}\label{lemma:Z}
    Consider the following two SDE with different initial condition
    \[
        \begin{cases}
        d\bs{x}_t = \bs{a}(\bs{x}_t) dt + \bs{\sigma} d \bs{W}_t,\\
        \bs{x}(0) = \bs{x}_0
        \end{cases}\quad
        \begin{cases}
        d\bs{y}_t = \bs{a}(\bs{y}_t) dt + \bs{\sigma} d \bs{W}_t,\\
        \bs{y}(0) = \bs{y}_0
        \end{cases}
    \]
    where $\bs{a}(\bs{u}) \in \mathbb{R}^d$ is $L$-Lipschitz, and $\bs{\sigma} \in \mathbb{R}^{n \times n}$ is a constant matrix. For $0 < h < \frac{1}{4L}$, we have the following representation
    \[
        \bs{x}_h - \bs{y}_h = \bs{x}_0 - \bs{y}_0 + \bs{z}
    \]
    with
    \[
        E\norm{\bs{z}}^2 \le 2L^2 \norm{\bs{x}_0 - \bs{y}_0}^2 h^2
    \]
\end{lemma}
\begin{proof}
Let $\bs{z} = (\bs{x}_h - \bs{y}_h) - (\bs{x}_0 - \bs{y}_0) = \int_0^h \bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s) ds$.
Ito's lemma readily implies that
\begin{align*}
    \mathbb{E}\norm{\bs{x}_h - \bs{y}_h}^2 =& \norm{\bs{x}_0 - \bs{y}_0}^2 + 2 \mathbb{E} \int_0^h \innerprod{\bs{x}_s - \bs{y}_s}{ \bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s)} ds \\
    \le& \norm{\bs{x}_0 - \bs{y}_0}^2 + 2L \int_0^h \mathbb{E} \norm{\bs{x}_s - \bs{y}_s}^2 ds
\end{align*}
By Gronwall's inequality, it follows that
\[
    \mathbb{E} \norm{\bs{x}_h - \bs{y}_h}^2 \le \norm{\bs{x}_0 - \bs{y}_0}^2 e^{2Lh} \le 2\norm{\bs{x}_0 - \bs{y}_0}^2, \mbox{ for } 0 < h < \frac{1}{4L}
\]
and
\begin{align*}
    \mathbb{E}\norm{\bs{z}}^2 =  \norm{ \mathbb{E} \left[\int_0^h \bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s) ds \right] }^2
    \le&  \left( \int_0^h \norm{\mathbb{E} \left[\bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s) \right]} ds \right) ^2  \\
    \le&  \int_0^h 1^2 ds \int_0^h \norm{\mathbb{E} \left[\bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s)\right]}^2 ds \\
    \le& h  \int_0^h \mathbb{E}  \norm{\bs{a}(\bs{x}_s) - \bs{a}(\bs{y}_s)}^2 ds   \\
    \le& L^2 h \int_0^h \mathbb{E} \norm{\bs{x}_s - \bs{y}_s}^2 ds \\
    \le& 2L^2 \norm{\bs{x}_0 - \bs{y}_0}^2 h^2
\end{align*}
\end{proof}

\subsection{Growth bound of SDE with additive noise}
\begin{lemma}\label{lemma:growth-bound}
Consider the following SDE with constant diffusion
\[
    \begin{cases}
    d\bs{x}_t = \bs{a}(\bs{x}_t) dt + \bs{\sigma} d \bs{W}_t,\\
    \bs{x}(0) = \bs{x}_0
    \end{cases}
\]
where $\bs{a}(\bs{x}) \in \mathbb{R}^d$ is $L$-smooth, i.e., $|\bs{a}(\bs{y} ) - \bs{a}(\bs{x})| \le L |\bs{y} - \bs{x}|$, $\bs{a}(\bs{0}) = \bs{0}$  and $\bs{\sigma} \in \mathbb{R}^{d \times d}$ is a constant matrix independent of time $t$ and $\bs{x}_t$. Then for $0 < h < \frac{1}{4L}$, we have
\[
     \mathbb{E} \norm{\bs{x}_h - \bs{x}_0}^2 \le 2.57 \left( \|\bs{\sigma}\|_F^2 + 2hL^2\norm{\bs{x}_0}^2 \right) h.
\]
\end{lemma}
\begin{proof}
We have
\begin{align*}
    \mathbb{E}\norm{\bs{x}_h - \bs{x}_0}^2 =& \mathbb{E} \norm{ \int_0^h \bs{a}(\bs{x}_t)dt + \int_0^h \bs{\sigma} d\bs{W}_t }^2\\
    \le& 2\mathbb{E}\norm{\int_0^h \bs{a}(\bs{x}_t)dt}^2 + 2\mathbb{E} \norm{\int_0^h \bs{\sigma}d\bs{W}_t}^2 \\
    \stackrel{(i)}{=}& 2\mathbb{E}\norm{\int_0^h \bs{a}(\bs{x}_t)dt}^2 + 2 \int_0^h \|\bs{\sigma}\|_F^2 dt \\
    \le& 2\mathbb{E}\left[ \left(\int_0^h \norm{\bs{a}(\bs{x}_t)} dt \right)^2 \right] + 2h \|\bs{\sigma}\|_F^2  \\
    \le& 2\mathbb{E}\left[ \left(\int_0^h \norm{\bs{a}(\bs{x}_t) - \bs{a}(\bs{x}_0)} dt + \int_0^h \norm{\bs{a}(\bs{x}_0)} dt \right)^2 \right] + 2h \|\bs{\sigma}\|_F^2 \\
    \le& 2\mathbb{E}\left[ \left(L \int_0^h \norm{\bs{x}_t - \bs{x}_0} dt + h\norm{\bs{a}(\bs{x}_0)}\right)^2 \right] + 2h \|\bs{\sigma}\|_F^2 \\
    \le& 4\mathbb{E}\left[ L^2\left( \int_0^h \norm{\bs{x}_t - \bs{x}_0} dt \right)^2 + h^2\norm{\bs{a}(\bs{x}_0)}^2 \right] + 2h \|\bs{\sigma}\|_F^2 \\
    \stackrel{(ii)}{\le}& 2h \|\bs{\sigma}\|_F^2 + 4h^2\norm{\bs{a}(\bs{x}_0)}^2 + 4L^2 h \int_0^h \mathbb{E}\norm{\bs{x}_t - \bs{x}_0 }^2dt
\end{align*}
where $(i)$ is due to Ito's isometry, $(ii)$ is due to Cauchy-Schwarz inequality and $\|\bs{\sigma}\|_F$ is the Frobenius norm of $\bs{\sigma}$. By Gronwall's inequality, we obtain
\[
    \mathbb{E} \norm{\bs{x}_h - \bs{x}_0}^2 \le \left( 2h \|\bs{\sigma}\|_F^2 + 4h^2\norm{\bs{a}(\bs{x}_0)}^2 \right) \exp\left\{ 4L^2h^2 \right\}.
\]
Since $\norm{\bs{a}(\bs{x}_0)} = \norm{\bs{a}(\bs{x}_0) - \bs{a}(\bs{0})} \le L \norm{\bs{x}_0}$, when $0 < h < \frac{1}{4L}$, we finally reach at
\[
    \mathbb{E} \norm{\bs{x}_h - \bs{x}_0}^2 \le 2\left( \|\bs{\sigma}\|_F^2 + 2hL^2\norm{\bs{x}_0}^2 \right) e^\frac{1}{4} h \le 2.57 \left( \|\bs{\sigma}\|_F^2 + 2hL^2\norm{\bs{x}_0}^2 \right) h.
\]
\end{proof}

\subsection{Lipschitz continuity of the drift of HFHR dynamics}
\begin{lemma}\label{lemma:lipschitz}
Assume $\nabla f$ is $L$-Lipschitz, i.e. $\norm{\nabla f(\bs{x}) - \nabla f(\bs{y})} \le L \norm{\bs{x} - \bs{y}}$, then the drift term of HFHR dynamics
\[
    \begin{bmatrix} \bs{p} - \alpha \nabla f(\bs{q}) \\ -\gamma \bs{p} - \nabla f(\bs{q}) \end{bmatrix}
\]
is $L^\prime$-Lipschitz, where $L^\prime \triangleq \sqrt{2} \max\{ \sqrt{1 + \alpha^2} \max\{\frac{1}{\sqrt{2}}, L\}, \sqrt{1 + \gamma^2}\}$. Let $P$ be defined in Appendix \ref{app:notation} and $\begin{bmatrix} \bs{\phi} \\ \bs{\psi} \end{bmatrix} = P \begin{bmatrix} \bs{q} \\ \bs{p} \end{bmatrix}$, then $\begin{bmatrix} \bs{\phi} \\ \bs{\psi} \end{bmatrix}$ satisfies the following SDE
\[
    \begin{bmatrix} d\bs{\phi} \\ d\bs{\psi} \end{bmatrix} = P \begin{bmatrix} \bs{p}(\bs{\phi}, \bs{\psi}) - \alpha \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \\ -\gamma \bs{p}(\bs{\phi}, \bs{\psi}) - \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \end{bmatrix} dt + P \begin{bmatrix} \sqrt{2\alpha} I & 0 \\ 0 & \sqrt{2\gamma} I \end{bmatrix} \begin{bmatrix} d\bs{W} \\ d\bs{B} \end{bmatrix}
\]
and the drift term
\[
    P \begin{bmatrix} \bs{p}(\bs{\phi}, \bs{\psi}) - \alpha \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \\ -\gamma \bs{p}(\bs{\phi}, \bs{\psi}) - \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \end{bmatrix}
\]
is $L^{\prime \prime}$-Lipschitz, where $L^{\prime\prime} = \kappa^\prime L^\prime$ and $\kappa^\prime$ is the condition number of $P$.
\end{lemma}
\begin{proof}
By direct computation and Cauchy-Schwarz inequality, we have
\begin{align*}
    &\norm{\begin{bmatrix} \bs{p}_1 - \alpha \nabla f(\bs{q}_1) \\ -\gamma \bs{p}_1 - \nabla f(\bs{q}_1) \end{bmatrix} - \begin{bmatrix} \bs{p}_2 - \alpha \nabla f(\bs{q}_2) \\ -\gamma \bs{p}_2 - \nabla f(\bs{q}_2) \end{bmatrix}} \\

    =& \sqrt{\norm{-\alpha \left( \nabla f(\bs{q}_1) - \nabla f(\bs{q}_2) \right) + (\bs{p}_1 - \bs{p}_2) }^2 + \norm{-\left( \nabla f(\bs{q}_1) - \nabla f(\bs{q}_2) \right) -\gamma (\bs{p}_1 - \bs{p}_2)}^2} \\
    \le& \sqrt{2\alpha^2\norm{\nabla f(\bs{q}_1) - \nabla f(\bs{q}_2)} + 2\norm{\bs{p}_1 - \bs{p}_2 }^2 + 2\norm{ \nabla f(\bs{q}_1) - \nabla f(\bs{q}_2)} + 2\gamma^2 \norm{\bs{p}_1 - \bs{p}_2}^2} \\
    \le& \sqrt{(2\alpha^2L^2 + 2L^2)\norm{\bs{q}_1 -\bs{q}_2} + (2 + 2 \gamma^2)\norm{\bs{p}_1 - \bs{p}_2 }^2} \\
    \le& \sqrt{2} \max\{L \sqrt{1 + \alpha^2}, \sqrt{1 + \gamma^2}\} \norm{\begin{bmatrix} \bs{q}_1 - \bs{q}_2 \\ \bs{p}_1\ - \bs{p}_2 \end{bmatrix}}  \\
    \le& \sqrt{2} \max\{ \sqrt{1 + \alpha^2} \max\{\frac{1}{\sqrt{2}}, L \} , \sqrt{1 + \gamma^2}\} \norm{\begin{bmatrix} \bs{q}_1 - \bs{q}_2 \\ \bs{p}_1\ - \bs{p}_2 \end{bmatrix}} \\
    \triangleq & L^\prime \norm{\begin{bmatrix} \bs{q}_1 - \bs{q}_2 \\ \bs{p}_1\ - \bs{p}_2 \end{bmatrix}}
\end{align*}
By Ito's lemma, we have
\[
    \begin{bmatrix} d\bs{\phi} \\ d\bs{\psi} \end{bmatrix} = P \begin{bmatrix} \bs{p}(\bs{\phi}, \bs{\psi}) - \alpha \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \\ -\gamma \bs{p}(\bs{\phi}, \bs{\psi}) - \nabla f(\bs{q}(\bs{\phi}, \bs{\psi})) \end{bmatrix} dt + P \begin{bmatrix} \sqrt{2\alpha} I & 0 \\ 0 & \sqrt{2\gamma} I \end{bmatrix} \begin{bmatrix} d\bs{W} \\ d\bs{B} \end{bmatrix}
\]
Using the Lipschitz constant obtained for the drift of HFHR, we further have
\begin{align*}
    & \norm{    P \begin{bmatrix} \bs{p}(\bs{\phi}_1, \bs{\psi}_1) - \alpha \nabla f(\bs{q}(\bs{\phi}_1, \bs{\psi}_1)) \\ -\gamma \bs{p}(\bs{\phi}_1, \bs{\psi}_1) - \nabla f(\bs{q}(\bs{\phi}_1, \bs{\psi}_1)) \end{bmatrix} -     P \begin{bmatrix} \bs{p}(\bs{\phi}_2, \bs{\psi}_2) - \alpha \nabla f(\bs{q}(\bs{\phi}_2, \bs{\psi}_2)) \\ -\gamma \bs{p}(\bs{\phi}_2, \bs{\psi}_2) - \nabla f(\bs{q}(\bs{\phi}_2, \bs{\psi}_2)) \end{bmatrix}} \\
    \le& \sigma_{\text{max}} \norm{\begin{bmatrix} \bs{p}_1 - \alpha \nabla f(\bs{q}_1) \\ -\gamma \bs{p}_1 - \nabla f(\bs{q}_1) \end{bmatrix} - \begin{bmatrix} \bs{p}_2 - \alpha \nabla f(\bs{q}_2) \\ -\gamma \bs{p}_2 - \nabla f(\bs{q}_2) \end{bmatrix}} \\
    \le& \sigma_{\text{max}} L^\prime \norm{\begin{bmatrix} \bs{q}_1 - \bs{q}_2 \\ \bs{p}_1\ - \bs{p}_2 \end{bmatrix}} \\
    \le& \sigma_{\text{max}} L^\prime \norm{P^{-1}\begin{bmatrix} \bs{\phi}_1 - \bs{\phi}_2 \\ \bs{\psi}_1\ - \bs{\psi}_2 \end{bmatrix}} \\
    \le& \sigma_{\text{max}} L^\prime \frac{1}{\sigma_{\text{min}}} \norm{\begin{bmatrix} \bs{\phi}_1 - \bs{\phi}_2 \\ \bs{\psi}_1\ - \bs{\psi}_2 \end{bmatrix}} \\
    =& \kappa^\prime L^\prime \norm{\begin{bmatrix} \bs{\phi}_1 - \bs{\phi}_2 \\ \bs{\psi}_1\ - \bs{\psi}_2 \end{bmatrix}}
\end{align*}
where $\sigma_{\text{max}}, \sigma_{\text{min}}$ and $\kappa^\prime$ are the largest, smallest singular values and the condition number (w.r.t. 2-norm) of matrix $P$.
\end{proof}
\begin{remark}
The following inequalities associated with $L^\prime$ will turn out to be useful in many proofs
\[
    L^\prime \ge 1, \, L^\prime \ge \sqrt{2}\gamma, L^\prime \ge \sqrt{2} \alpha, L \ge \sqrt{2} L \mbox{ and } L^\prime \ge \sqrt{2} \alpha L.
\]
\end{remark}

\subsection{Contraction of (Transformed) HFHR Dynamics}
\begin{lemma}\label{lemma:contraction}
Suppose $f$ is $L$-smooth, $m$-strongly convex and $\gamma^2 > L$. Consider two copies of HFHR dynamics $\begin{bmatrix}\qbold_t \\ \pbold_t \end{bmatrix}$, $\begin{bmatrix}\tqbold_t \\ \tpbold_t \end{bmatrix}$ (driven by the same Brownian motion) with initialization $\begin{bmatrix}\qbold_0 \\ \pbold_0 \end{bmatrix}$, $\begin{bmatrix}\tqbold_0 \\ \tpbold_0 \end{bmatrix}$ respectively, then we have
\[
   \norm{P\begin{bmatrix}
   \qbold_t - \tqbold_t \\
   \pbold_t - \tpbold_t
   \end{bmatrix}} \le e^{-\lambda^\prime t} \norm{P\begin{bmatrix}
   \qbold_0 - \tqbold_0 \\
   \pbold_0 - \tpbold_0
   \end{bmatrix}}
\]
where $P = \begin{bmatrix} \gamma I & I \\ 0 & \sqrt{1+\alpha\gamma}I\end{bmatrix}$ and $\lambda^\prime = \min\{ \frac{m}{\gamma} + \alpha m, \frac{\gamma^2 - L}{\gamma} \}$.
\end{lemma}
\begin{proof}
Consider two copies of HFHR that are driven by the same Brownian motion
\[
    \begin{cases}
        d\qbold_t = (\pbold_t - \alpha \nabla f(\qbold_t))dt + \sqrt{2\alpha} d\bs{B}^1_t \\
        d\pbold_t = (-\gamma \pbold_t - \nabla f(\qbold_t))dt + \sqrt{2\gamma} d\bs{B}^2_t
    \end{cases},
    \quad
    \begin{cases}
        d\tilde{\qbold}_t = (\tpbold_t - \alpha \nabla f(\tilde{\qbold}_t))dt + \sqrt{2\alpha} d\bs{B}^1_t \\
        d\tilde{\pbold_t} = (-\gamma \tpbold_t - \nabla f(\tilde{\qbold}_t))dt + \sqrt{2\gamma} d\bs{B}^2_t
    \end{cases}.
\]
The difference of the two copies satisfies the following equation
\begin{align*}
    \frac{d}{dt}\begin{bmatrix}\qbold_t - \tilde{\qbold}_t \\ \pbold_t - \tilde{\pbold_t}\end{bmatrix}
    =& -
    \begin{bmatrix}
    \alpha H_t & -I \\
    H_t & \gamma I
    \end{bmatrix}
    \begin{bmatrix}\qbold_t - \tilde{\qbold}_t \\ \pbold_t - \tilde{\pbold_t}\end{bmatrix} \triangleq - A \begin{bmatrix}\qbold_t - \tilde{\qbold}_t \\ \pbold_t - \tilde{\pbold_t}\end{bmatrix}
\end{align*}
where $H_t = \int_0^1 \nabla^2 f(\tilde{\qbold}_t + s(\qbold - \tilde{\qbold}_t)) ds$. Denote the eigenvalues of $H_t$ by $\eta_i, 1\le i \le d$, by strong convexity and smoothness assumption on $f$, we have $m \le \eta_i \le L, 1\le i \le d$.

Denote $\phipsi = P\deltaqp$ and consider $\mathcal{L}_t = \frac{1}{2}\norm{\phipsi}^2$, we have
\begin{align*}
    \frac{d}{dt}\mathcal{L}_t =& - \phipsi^T P A P^{-1} \phipsi \\
    =& - \phipsi^T \frac{1}{2}(PAP^{-1} + (P^{-1})^T A^T P^T) \phipsi \\
    =& - \phipsi^T \frac{1}{\gamma} \begin{bmatrix} (1+\alpha\gamma)H_t & 0_{d\times d} \\ 0_{d\times d} & \gamma^2I - H_t \end{bmatrix} \phipsi \\
    \triangleq& -\phipsi^T B(\alpha) \phipsi
\end{align*}
It is easy to see that
\begin{align*}
    \lambda_{\min}(B(\alpha)) = \min_{i=1,2,\cdots,d} \{ \min\{ \frac{\eta_i}{\gamma} + \alpha \eta_i, \gamma - \frac{\eta_i}{\gamma}\} \} \ge \min\{ \frac{m}{\gamma} + \alpha m, \frac{\gamma^2 - L}{\gamma} \}  \triangleq \lambda^\prime.
\end{align*}
Therefore we have
$
    \frac{d}{dt}\mathcal{L}_t \le - 2\lambda_{\min}{B(\alpha)} \mathcal{L}_t \le - 2 \lambda^\prime \mathcal{L}_t
$. By Gronwall's inequality, we obtain
\[
    \norm{\phipsi}^2 \le e^{-2\lambda^\prime t} \norm{\begin{bmatrix} \bs{\phi}_0 \\ \bs{\psi}_0 \end{bmatrix}}^2.
\]
and the desired inequality follows by taking square root.
\end{proof}

\subsection{Local error between the exact Strang's splitting method and HFHR dynamics}
\begin{lemma}\label{lemma:local-exact}
Assume $f$ is $L$-smooth and $\bs{0} \in \argmin_{\bs{x} \in \mathbb{R}^d} f(\bs{x})$, i.e. $\nabla f(\bs{0}) = \bs{0}$. If $0 < h \le  \frac{1}{4 L^\prime} $, then compared with the HFHR dynamics, the exact Strang's splitting method has local mathematical expectation of deviation of order $p_1 = 2$ and local mean-squared error of order $p_2 = 2$, i.e. there exist constants $\widehat{C}_1, \widehat{C}_2 > 0$ such that
\begin{equation*}
    \norm{\mathbb{E} \bs{x}(h)  - \mathbb{E} \hat{\bs{x}}(h) } \le \widehat{C}_1  h^{p_1}
\end{equation*}
\begin{equation*}
    \left(\mathbb{E} \left[ \norm{ \bs{x}(h) - \hat{\bs{x}}(h) }^2 \right]\right)^\frac{1}{2} \le \widehat{C}_2 h^{p_2}
\end{equation*}
where $\bs{x}(h) = \begin{bmatrix} \bs{q}(h) \\ \bs{p}(h) \end{bmatrix}$ is the solution of the HFHR dynamics with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$ and $\hat{\bs{x}}(h) = \begin{bmatrix} \hat{\bs{q}}(h) \\ \hat{\bs{p}}(h) \end{bmatrix}$ is the solution of the implementable Strang's splitting with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$, $p_1 = 2$ and $p_2 = 2$. More concretely, we have
\[
  \widehat{C}_1 = L \max\{\alpha + 1.25, \gamma + 1\} \left(1.74 \norm{\bs{x}_0} + ( 1.26\sqrt{\alpha} + 2.84\sqrt{\gamma})\sqrt{hd}\right),
\]
\[
    \widehat{C}_2 =  L \max\{ \alpha + 1.25, \gamma + 1\} \left( 1.92 \norm{\bs{x}_0} + (1.30\sqrt{\alpha} + 3.22\sqrt{\gamma})\sqrt{hd} \right).
\]
\end{lemma}
\begin{proof}
The exact Strang's splitting integrator with step size $h$ reads as $\phi^{\frac{h}{2}} \circ \psi^{h} \circ \phi^{\frac{h}{2}}$ where
\begin{equation*}
    \phi:
    \begin{cases}
        d\bs{q} = \bs{p} dt \\
        d\bs{p} = -\gamma \bs{p} dt + \sqrt{2\gamma} d\bs{B}
    \end{cases}
    \,
    \psi:
    \begin{cases}
        d\bs{q} = -\alpha \nabla f(\bs{q}) dt + \sqrt{2\alpha} d\bs{W} \\
        d\bs{p} = -\nabla f(\bs{q}) dt
    \end{cases}.
\end{equation*}
The $\phi$ flow can be explicitly solved and the solution is
\begin{equation*}
    \begin{cases}
        \bs{q}(t) = \bs{q}_0 + \frac{1 - e^{-\gamma t}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^t \frac{1 - e^{-\gamma(t - s)}}{\gamma} d\bs{B}(s)\\
        \bs{p}(t) = e^{-\gamma t}\bs{p}_0 + \sqrt{2\gamma} \int_0^t e^{-\gamma (t - s)} d\bs{B}(s)
    \end{cases}.
\end{equation*}
The $\psi$ flow can be written as
\begin{equation*}
    \begin{cases}
        \bs{q}(t) = \bs{q}_0 - \int_0^t \alpha \nabla f(\bs{q}(s)) ds + \sqrt{2\alpha} \int_0^t d\bs{W}(s) \\
        \bs{p}(t) = \bs{p}_0 - \int_0^t \nabla f(\bs{q}(s)) ds
    \end{cases}.
\end{equation*}

The solution of one-step exact Strang's  splitting integrator with step size $h$ can be written as
\[
    \begin{cases}
    \bs{q}_3 = \bs{q}_2(h) + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_2(h) + \sqrt{2\gamma} \int_\frac{h}{2}^h \frac{1 - e^{-\gamma(h - s)}}{\gamma} d\bs{B}(s)\\
    \bs{p}_3 = e^{-\gamma \frac{h}{2}}\bs{p}_2(h) + \sqrt{2\gamma} \int_\frac{h}{2}^h e^{-\gamma (h - s)} d\bs{B}(s) \\
    \bs{q}_2(r) = \bs{q}_1 - \int_0^r \alpha \nabla f(\bs{q}_2(s)) ds + \sqrt{2\alpha}\int_0^r d\bs{W}(s)  \qquad (0 \le r \le h)\\
    \bs{p}_2(r) = \bs{p}_1 - \int_0^r \nabla f(\bs{q}_2(s)) ds \\
    \bs{q}_1 = \bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(s)\\
    \bs{p}_1 = e^{-\gamma \frac{h}{2}}\bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)
    \end{cases}
\]
Therefore, we have $\hat{\bs{q}}(h) =  \bs{q}_3, \hat{\bs{p}}(h) = \bs{p}_3$ and
\begin{align*}
    \hat{\bs{q}}(h) =& \sqrt{2\gamma} \int_\frac{h}{2}^h \frac{1 - e^{-\gamma(h - s)}}{\gamma} d\bs{B}(s)
    + \underbrace{\bs{q}_1 - \int_0^h \alpha \nabla f(\bs{q}_2(s)) ds + \sqrt{2\alpha}\int_0^h d\bs{W}(s)}_{\bs{q}_2(h)} \\
    & + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \left[ \underbrace{\bs{p}_1 - \int_0^h \nabla f(\bs{q}_2(s)) ds}_{\bs{p}_2(h)} \right] \\
    =&  \sqrt{2\gamma} \int_\frac{h}{2}^h \frac{1 - e^{-\gamma(h - s)}}{\gamma} d\bs{B}(s)  - \int_0^h \alpha \nabla f(\bs{q}_2(s)) ds + \sqrt{2\alpha}\int_0^h d\bs{W}(s) - \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \int_0^h \nabla f(\bs{q}_2(s)) ds \\
    & +\underbrace{\bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(s)}_{\bs{q}_1} + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \left[ \underbrace{e^{-\gamma \frac{h}{2}}\bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)}_{\bs{p}_1} \right] \\
    =& \bs{q}_0  + \frac{1 - e^{-\gamma h}}{\gamma}  \bs{p}_0 - \left(\alpha +  \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma}\right) \int_0^h \nabla f(\bs{q}_2(s))  ds \\
    &+ \sqrt{2\alpha}\int_0^h d\bs{W}(s) + \sqrt{2\gamma} \int_{\frac{h}{2}}^h \frac{1 - e^{-\gamma(h - s)}}{\gamma} d\bs{B}(s) + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma (\frac{h}{2} -s)}}{\gamma}d\bs{B}(s) \\
    & + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma}  \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)  \\
    \hat{\bs{p}}(h) =& e^{-\gamma \frac{h}{2}} \left[ \underbrace{\bs{p}_1 - \int_0^h \nabla f(\bs{q}_2(s)) ds}_{\bs{p}_2(h)} \right] + \sqrt{2\gamma} \int_{\frac{h}{2}}^h e^{-\gamma (h - s)} d\bs{B}(s)  \\
    =& e^{-\gamma \frac{h}{2}} \left[ \underbrace{e^{-\gamma \frac{h}{2}}\bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)}_{\bs{p}_1} \right] - e^{-\gamma \frac{h}{2}} \int_0^h \nabla f(\bs{q}_2(s)) ds + \sqrt{2\gamma} \int_{\frac{h}{2}}^h e^{-\gamma (h - s)} d\bs{B}(s) \\
    =& e^{-\gamma h} \bs{p}_0 - e^{-\gamma \frac{h}{2}} \int_0^h \nabla f(\bs{q}_2(s)) ds + e^{-\gamma \frac{h}{2}} \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s) + \sqrt{2\gamma} \int_{\frac{h}{2}}^h e^{-\gamma (h - s)} d\bs{B}(s)
\end{align*}
It is clear that $\hat{\bs{q}}(h), \hat{\bs{p}}(h)$ should be compared with the exact solution of HFHR at time $h$, which can be written as
\begin{align*}
    \bs{q}(h) =& \bs{q}_0
    + \frac{1 - e^{-\gamma h}}{\gamma} \bs{p}_0
    - \int_0^h \left(\frac{1 - e^{-\gamma (h - s)}}{\gamma} + \alpha \right) \nabla f(\bs{q}(s)) ds
    + \sqrt{2\alpha} \int_0^h d\bs{W}_s
    + \sqrt{2\gamma} \int_0^h \frac{1 - e^{-\gamma (h-s)}}{\gamma} d\bs{B}_s\\
    \bs{p}(h) =& e^{-\gamma h} \bs{p}_0 - \int_0^h e^{-\gamma(h - s)} \nabla f(\bs{q}(s)) ds + \sqrt{2\gamma} \int_0^h e^{-\gamma(h -s )} d\bs{B}(s)
\end{align*}
Subtracting $\bs{q}(h), \bs{p}(h)$ from $\hat{\bs{q}}(h), \hat{\bs{p}}(h)$ respectively, we obtain
\begin{align*}
    \hat{\bs{q}}(h) - \bs{q}(h) =& -\left( \alpha +  \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \right) \int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds  \\
    & + \int_0^h \left( \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right) \nabla f(\bs{q}(s)) ds \\
    \hat{\bs{p}}(h) - \bs{p}(h) =& -e^{-\gamma \frac{h}{2}}\int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds
    + \int_0^h \left( e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right) \nabla f(\bs{q}(s))ds
\end{align*}
It should be clear now that we will need to bound the term $\nabla f(\bs{q}_2) - \nabla f(\bs{q})$ and $\nabla f(\bs{q})$. Since
\begin{align*}
    \bs{q}_2(r) =& \bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma}\bs{p}_0 + \sqrt{2\gamma}\int_0^\frac{h}{2} \frac{1 - e^{-\gamma (\frac{h}{2} - s)}}{\gamma}d\bs{B}(s) - \alpha \int_0^r \nabla f(\bs{q}_2(s))ds + \sqrt{2\alpha} \int_0^r d\bs{W}(s) \\
    \bs{q}(r) =& \bs{q}_0
    + \frac{1 - e^{-\gamma r}}{\gamma} \bs{p}_0
    - \int_0^r \left(\frac{1 - e^{-\gamma (r - s)}}{\gamma} +\alpha \right) \nabla f(\bs{q}(s)) ds
    + \sqrt{2\alpha} \int_0^r d\bs{W}(s) \\
    &+ \sqrt{2\gamma} \int_0^r \frac{1 - e^{-\gamma (r-s)}}{\gamma} d\bs{B}(s) \nonumber,
\end{align*}
we then have
\begin{align*}
    \bs{q}_2(r) - \bs{q}(r) =& \frac{e^{-\gamma r} - e^{-\gamma \frac{h}{2}}}{\gamma}\bs{p}_0 - \alpha \int_0^r \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds + \int_0^r \frac{1 - e^{-\gamma (r - s)}}{\gamma} \nabla f(\bs{q}(s)) ds \\
    +& \sqrt{2\gamma}\int_0^\frac{h}{2} \frac{1 - e^{-\gamma (\frac{h}{2} - s)}}{\gamma}d\bs{B}(s) - \sqrt{2\gamma} \int_0^r \frac{1 - e^{-\gamma (r - s)}}{\gamma} d\bs{B}(s) \nonumber
\end{align*}
By Lemma \ref{lemma:lipschitz} and \ref{lemma:growth-bound}, when $0 < h < \frac{1}{4L^\prime}$, we have the following for the solution of HFHR dynamics
\[
    \mathbb{E}[\norm{\bs{x}_{0, \bs{x}_0}(h) - \bs{x}_0}^2] \le \widehat{C}_0 h
\]
where $\widehat{C}_0 = 5.14 \left\{(\alpha + \gamma)d +  h\left( L^\prime \right)^2 \norm{\bs{x}_0}^2 \right\}$ and hence
\begin{align}
    \mathbb{E}\left[ \int_0^r \norm{\nabla f(\bs{q}(s))}^2 ds \right]
    \le& \mathbb{E} \left[ 2 \int_0^r \norm{\nabla f (\bs{q}(0))}^2 ds + 2\int_0^r \norm{\nabla f(\bs{q}(s)) - \nabla f(\bs{q}(0))}^2ds  \right] \nonumber \\
    \le& \mathbb{E}\left[ 2L^2r \norm{\bs{q}(0)}^2  + 2L^2 \int_0^r \norm{\bs{q}(s) - \bs{q}(0)}^2ds \right] \nonumber\\
    \le& 2L^2 r \norm{\bs{x}_0}^2 + 2L^2 \mathbb{E}\left[ \int_0^r \norm{\bs{q}(s) - \bs{q}(0)}^2ds \right] \nonumber\\
    \le& 2L^2 r \norm{\bs{x}_0}^2 + 2L^2\widehat{C}_0 \int_0^r sds \nonumber\\
    \le& L^2 r \left( 2\norm{\bs{x}_0}^2 + h \widehat{C}_0 \right) \nonumber\\
    \le& L^2 r \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right) \label{eq:tool1}
\end{align}
Now $\mathbb{E}\left[ \norm{\bs{q}_2 - \bs{q} }^2 \right]$ can be bounded as follow
\begin{align*}
    & \mathbb{E} \left[\norm{ \bs{q}_2(r) - \bs{q}(r) }^2 \right] \\
    \le& 5 \left\{ \left(\frac{e^{-\gamma r} - e^{-\gamma \frac{h}{2}}}{\gamma} \right)^2\norm{\bs{p}_0}^2  + \alpha^2 \mathbb{E} \norm{\int_0^r \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds }^2 + \mathbb{E} \norm{ \int_0^r \frac{1 - e^{-\gamma(r-s)}}{\gamma} \nabla f(\bs{q}(s)) ds }^2 \right\} \\
    & + 5\left\{ 2\gamma \mathbb{E} \norm{ \int_0^\frac{h}{2} \frac{1 - e^{-\gamma (\frac{h}{2} - s)}}{\gamma} d\bs{B}(s) }^2 + 2\gamma \mathbb{E} \norm{ \int_0^r \frac{1 - e^{-\gamma (r - s)}}{\gamma} d\bs{B}(s) }^2\right\} \qquad \mbox{(Cauchy-Schwartz Inequality)}\nonumber \\
    \le& 5 \left\{\frac{h^2}{4}\norm{\bs{x}_0}^2  + \alpha^2 L^2 r  \int_0^r \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) }^2 ds +  \int_0^r \left(\frac{1 - e^{-\gamma(r-s)}}{\gamma} \right)^2 ds \int_0^r \mathbb{E}  \norm{\nabla f(\bs{q}(s)) }^2 ds\right\} \\
    &+ 5\left\{ \frac{\gamma d h^3}{12} + \frac{2\gamma d}{3} r^3  \right\} \\
    \le& 5 \left\{ \frac{h^2}{4} \norm{ \bs{x}_0 }^2 + \alpha^2 L^2 r \int_0^r \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) }^2 ds + \frac{h^3}{3} \mathbb{E} \left[ \int_0^r \norm{ \nabla f(\bs{q}(s)) }^2 \right] + \frac{3\gamma d}{4} h^3
    \right\} \\
    \le& 5 \left\{
    \frac{h^2}{4} \norm{ \bs{x}_0 }^2
    + \frac{3\gamma d}{4} h^3 + \frac{h^3}{3} L^2 \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right) r
    + \alpha^2 L^2 r \int_0^r \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) }^2 ds
    \right\} \\
    \le& 5 h^2 \left\{
    \frac{1}{4} \norm{ \bs{x}_0 }^2
    + \frac{3\gamma d}{4} h
    + \frac{h^2}{3} L^2 \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right) \right\}
    + 5 \alpha^2 L^2 h \int_0^r \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) }^2 ds
\end{align*}
By Gronwall's inequality and $0 < h \le \frac{1}{4 L^\prime}$, we have
\begin{align}
    \mathbb{E} \left[  \norm{ \bs{q}_2(r) - \bs{q}(r) }^2 \right]
    \le& 5 h^2 \left\{ \frac{1}{4} \norm{\bs{x}_0}^2 + \frac{3\gamma d}{4} h + \frac{h^2}{3} L^2 \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right) \right\}  \exp\{5\alpha^2 L^2 h^2\} \nonumber\\
    \le& 5 h^2 \left\{ \frac{1}{4} \norm{\bs{x}_0}^2 + \frac{3\gamma d}{4} h + \frac{h^2}{3} L^2 \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right) \right\}  e^\frac{5}{32} \nonumber\\
    \le& 5.85 h^2 \left\{ 0.28 \norm{\bs{x}_0}^2 + (0.06\alpha + 0.81\gamma)h d \right\} \nonumber\\
    \le& h^2 \left\{ 1.64 \norm{\bs{x}_0}^2 + (0.36\alpha + 4.74\gamma)h d \right\}. \label{eq:tool2}
\end{align}
With bounds in Equation \eqref{eq:tool1} and \eqref{eq:tool2},  we are now ready to show $p_1$ and $p_2$. For $p_1$, i.e. the order of the mathematical expectation of deviation, we have
\begin{align*}
    &\norm{\mathbb{E}\left[ \begin{bmatrix}\hat{\bs{q}}(h) \\ \hat{\bs{p}}(h) \end{bmatrix} - \begin{bmatrix}\bs{q}(h) \\ \bs{p}(h) \end{bmatrix} \right] } \\
    \le& \norm{\mathbb{E}\left[ \hat{\bs{q}}(h)  -  \bs{q}(h) \right] } +  \norm{\mathbb{E}\left[ \hat{\bs{p}}(h)  -  \bs{p}(h) \right] } \\
    \le& \left( \alpha + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma}  \right) \norm{ \int_0^h \mathbb{E} \left[ \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) \right]ds } + \norm{ \int_0^h \left( \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right) \mathbb{E}\left[\nabla f(\bs{q}(s))\right] ds }\\
    & + e^{-\gamma\frac{ h}{2}} \norm{\int_0^\frac{h}{2} \mathbb{E} \left[ \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) \right] ds} + \norm{\int_0^h \left( e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right)  \mathbb{E} \left[\nabla f(\bs{q}(s)) \right] ds } \\
    \le& \left( \alpha + 1 + \frac{h}{2}  \right) L \int_0^h \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s)} ds \\
    & + \int_0^h \left(\left| \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right| + \left|  e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right| \right) \norm{ \mathbb{E}\left[\nabla f(\bs{q}(s))\right]} ds \\
    \le& L\left( \alpha + 1 + \frac{h}{2}  \right) \int_0^h \mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) } ds \\
    & + \left\{ \left(\int_0^h \left| \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right|^2 ds\right)^\frac{1}{2} + \left( \int_0^h \left|  e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right|^2 ds \right)^\frac{1}{2} \right\} \left( \int_0^h \norm{ \mathbb{E}\left[\nabla f(\bs{q}(s))\right]}^2 ds \right)^\frac{1}{2} \\
    \le& L\left( \alpha + 1 + \frac{h}{2}  \right) \int_0^h \left(\mathbb{E} \norm{ \bs{q}_2(s) - \bs{q}(s) }^2 \right)^\frac{1}{2} ds + \frac{1 + \gamma}{2\sqrt{3}} h^\frac{3}{2} \left( \mathbb{E} \int_0^h \norm{ \left[\nabla f(\bs{q}(s))\right]}^2 ds \right)^\frac{1}{2} \\
    \le& L\left( \alpha + 1 + \frac{h}{2}  \right) h^2 \left\{ 1.64 \norm{\bs{x}_0}^2 + (0.36\alpha + 4.74\gamma)h d \right\}^\frac{1}{2} + \frac{1 + \gamma}{2\sqrt{3}} h^2 L \left( 2.33\norm{\bs{x}_0}^2 + 5.14 (\alpha + \gamma) dh \right)^\frac{1}{2} \\
    \le& L\left( \alpha + 1.25 \right) h^2 \left( 1.29 \norm{\bs{x}_0} + \sqrt{0.36\alpha + 4.74\gamma} \sqrt{h d} \right) + (1 + \gamma) h^2 L \left( 0.45\norm{\bs{x}_0} + \sqrt{0.43\alpha + 0.43\gamma} \sqrt{dh} \right) \\
    \le& L h^2 \max\{\alpha + 1.25, \gamma + 1\} \left(1.74 \norm{\bs{x}_0} + ( 1.26\sqrt{\alpha} + 2.84\sqrt{\gamma})\sqrt{hd}\right)
\end{align*}
The above derivation proves $p_1 = 2$ with
\[
    \widehat{C}_1 = L \max\{\alpha + 1.25, \gamma + 1\} \left(1.74 \norm{\bs{x}_0} + ( 1.26\sqrt{\alpha} + 2.84\sqrt{\gamma})\sqrt{hd}\right).
\]

We now proceed with $p_2$, i.e. mean-square error
\begin{align*}
    &\mathbb{E} \norm{\begin{bmatrix}\hat{\bs{q}}(h) \\ \hat{\bs{p}}(h) \end{bmatrix} -  \begin{bmatrix}\bs{q}(h) \\ \bs{p}(h) \end{bmatrix}}^2 \\
    \le& 2\left( \alpha +  \frac{h}{2} \right)^2 \mathbb{E} \norm{\int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds}^2 + 2\mathbb{E} \norm{ \int_0^h \left( \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right) \nabla f(\bs{q}(s)) ds }^2 \\
    & + 2 \mathbb{E} \norm{ \int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}(s)) ds }^2 + 2 \mathbb{E} \norm{ \int_0^h \left( e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right) \nabla f(\bs{q}(s))ds }^2 \\
    \le& 2\left( (\alpha + \frac{h}{2})^2 + 1 \right) L^2 \mathbb{E} \left(\int_0^h |\bs{q}_2(s) - \bs{q}(s)| ds\right)^2 +  2\int_0^h \left| \frac{1 - e^{-\gamma (h - s)}}{\gamma} - \frac{1 - e^{-\gamma\frac{h}{2}}}{\gamma} \right|^2 ds \int_0^h \mathbb{E} \norm{\nabla f(\bs{q}(s))}^2 ds \\
    & + 2\int_0^h \left| e^{-\gamma(h - s)} - e^{-\gamma \frac{h}{2}} \right|^2 ds  \, \int_0^h \mathbb{E} \norm{ \nabla f(\bs{q}(s)) }^2 ds \\
    \le& 2\left( (\alpha + \frac{h}{2})^2 + 1 \right) L^2 h \int_0^h \mathbb{E}|\bs{q}_2(s) - \bs{q}(s)|^2 ds
    + \frac{1 + \gamma^2}{6} h^3 \int_0^h \mathbb{E} |\nabla f(\bs{q}(s))|^2 ds\\
    \le& 2\left( (\alpha + \frac{h}{2})^2 + 1 \right) L^2  \left\{ 1.64\norm{\bs{x}_0}^2 + (0.36\alpha + 4.74\gamma) hd\right\} h^4 + \frac{1 + \gamma^2}{6} L^2 \left\{ 2.33\norm{\bs{x}_0}^2 + 5.14(\alpha + \gamma) hd \right\} h^4 \\
    \le& L^2 \max\{ (\alpha + 1.25)^2, 1 + \gamma^2\} \left( 3.67 \norm{\bs{x}_0}^2 + (1.68\alpha + 10.34\gamma)hd \right) h^4
\end{align*}
The above derivation implies $p_2 = 2$ with
\[
    \widehat{C}_2 =  L \max\{ \alpha + 1.25, 1 + \gamma\} \left( 1.92 \norm{\bs{x}_0} + (1.30\sqrt{\alpha} + 3.22\sqrt{\gamma})\sqrt{hd} \right).
\]
\end{proof}

\subsection{Local error between Algorithm \ref{alg:HFHR} and the exact Strang's splitting method}
\begin{lemma}\label{lemma:local-implementable}
Assume $f$ is $L$-smooth, $ \bs{0} \in \argmin_{\bs{x} \in \mathbb{R}^d} f(\bs{x})$, i.e. $\nabla f(\bs{0}) = \bs{0}$ and the operator $\nabla \Delta f$ grows at most linearly, i.e. $\norm{\nabla \Delta f(\bs{q})} \le G\sqrt{1 + \norm{\bs{q}}^2 }$. If $0 < h \le \frac{1}{4 L^\prime}$, then compared with the exact Strang's splitting method of HFHR dynamics, the implementable Strang's splitting method has local mathematical expectation of deviation of order $p_1 = 2$ and local mean-squared error of order $p_2 = 1.5$, i.e. there exist constants $\bar{C}_1, \bar{C}_2 > 0$ such that
\[
    \norm{ \mathbb{E} \hat{\bs{x}}(h)  - \mathbb{E} \bar{\bs{x}}(h) } \le \bar{C}_1 h^{p_1}
\]
\[
    \left(\mathbb{E} \left[ \norm{ \hat{\bs{x}}(h) - \bar{\bs{x}}(h) }^2 \right]\right)^\frac{1}{2} \le \bar{C}_2 h^{p_2}
\]
where $\hat{\bs{x}}(h) = \begin{bmatrix} \hat{\bs{q}}(h) \\ \hat{\bs{p}}(h) \end{bmatrix}$ is the solution of the exact Strang's splitting method for HFHR with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$ and $\bar{\bs{x}}(h) = \begin{bmatrix} \bar{\bs{q}}(h) \\ \bar{\bs{p}}(h) \end{bmatrix}$ is the one-step result of Algorithm \ref{alg:HFHR} with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$, $p_1 = 2$ and $p_2 = 1.5$. More concretely, we have
\[
    \bar{C}_1 = \alpha (\alpha + 1.125) (L+G) \left[ 0.5 + 0.71\norm{\bs{x}_0} + (1.14\sqrt{\alpha} + 0.21 \sqrt{\gamma} h) \sqrt{hd} \right]
\]
and
\[
    \bar{C}_2 = L(\alpha + 0.73) \left( 2.30\sqrt{h} \alpha L \norm{\bs{x}_0} + (2.27\sqrt{\alpha} + 0.12 \sqrt{\gamma} h )\sqrt{d} \right).
\]
\end{lemma}

\begin{proof}
The solution of one-step exact Strang's  splitting integrator with step size $h$ can be written as
\begin{equation*}
    \begin{cases}
    \bs{q}_3 = \bs{q}_2(h) + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_2(h) + \sqrt{2\gamma} \int_\frac{h}{2}^h \frac{1 - e^{-\gamma(h - s)}}{\gamma} d\bs{B}(s)\\
    \bs{p}_3 = e^{-\gamma \frac{h}{2}}\bs{p}_2(h) + \sqrt{2\gamma} \int_\frac{h}{2}^h e^{-\gamma (h - s)} d\bs{B}(s) \\
    \bs{q}_2(r) = \bs{q}_1 - \int_0^r \alpha \nabla f(\bs{q}_2(s)) ds + \sqrt{2\alpha}\int_0^r d\bs{W}(s)  \qquad (0 \le r \le h)\\
    \bs{p}_2(r) = \bs{p}_1 - \int_0^r \nabla f(\bs{q}_2(s)) ds \\
    \bs{q}_1 = \bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(s)\\
    \bs{p}_1 = e^{-\gamma \frac{h}{2}}\bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)
    \end{cases}
\end{equation*}
and the solution of one-step implementable Strang's splitting integrator with step size $h$ can be written as
\begin{equation*}
    \begin{cases}
    \bar{\bs{q}}_3 = \bar{\bs{q}}_2(h) + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bar{\bs{p}}_2(h) + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(\frac{h}{2} + s)\\
    \bar{\bs{p}}_3 = e^{-\gamma \frac{h}{2}}\bar{\bs{p}}_2(h) + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(\frac{h}{2} + s) \\
    \bar{\bs{q}}_2(r) = \bs{q}_1 - \int_0^r \alpha \nabla f(\bs{q}_1) ds + \sqrt{2\alpha}\int_0^r d\bs{W}(s)  \qquad (0 \le r \le h)\\
    \bar{\bs{p}}_2(r) = \bs{p}_1 - \int_0^r \nabla f(\bs{q}_1) ds \\
    \bs{q}_1 = \bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(s)\\
    \bs{p}_1 = e^{-\gamma \frac{h}{2}}\bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} e^{-\gamma (\frac{h}{2} - s)} d\bs{B}(s)
    \end{cases}
\end{equation*}
Note that in the implementable Strang's splitting method, $\phi$ flow can be explicitly integrated and hence $\bs{q}_1, \bs{p}_1$ are the same as that in the exact Strang's splitting method.

First, we will bound the deviation of mathematical expectation and mean squared error of $\bs{q}_2(h) - \bar{\bs{q}}_2(h)$ and $\bs{p}_2(h) - \bar{\bs{p}}_2(h)$. We have
\begin{equation}\label{eq: tmp-1}
    \begin{cases}
        \bs{q}_2(h) - \bar{\bs{q}}_2(h) =& -\alpha \int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1) ds \\
        \bs{p}_2(h) - \bar{\bs{p}}_2(h) =& -\int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1) ds
    \end{cases}
\end{equation}

Square both sides of the first equation in \eqref{eq: tmp-1} and take expectation, we obtain
\begin{align*}
    \mathbb{E}\norm{\bs{q}_2(h) - \bar{\bs{q}}_2(h)}^2 =& \alpha^2 \mathbb{E}\norm{\int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1) ds}^2 \\
    \le& \alpha^2 \mathbb{E}\left(\int_0^h \norm{\nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1)} ds\right)^2 \\
    \le& \alpha^2 L^2 \mathbb{E}\left(\int_0^h \norm{\bs{q}_2(s) - \bs{q}_1} ds\right)^2 \\
    \le& \alpha^2 L^2 h \int_0^h \mathbb{E}\norm{\bs{q}_2(s) - \bs{q}_1}^2 ds
\end{align*}
Note that $\bs{q}_2$ is the solution of a rescaled overdamped Langevin dynamics whose drift vector field is $\alpha L$-Lipschitz, by conditional expectation version of Lemma \ref{lemma:growth-bound}, for $0 < h < \frac{1}{4L^\prime} < \frac{1}{4\alpha L}$, we have
$
    \mathbb{E}\norm{\bs{q}_2(h) - \bs{q}_1}^2 \le \bar{C}_0 h
$ with $\bar{C}_0 = 5.14\left\{ \alpha d +  h(\alpha L)^2 \mathbb{E}\norm{\bs{q}_1}^2 \right\}$ and it follows that
\[
    \begin{cases}
        \mathbb{E}\norm{\bs{q}_2(h) - \bar{\bs{q}}_2(h)}^2 \le& \alpha^2 L^2 \bar{C}_0  h^3 \\
        \mathbb{E}\norm{\bs{p}_2(h) - \bar{\bs{p}}_2(h)}^2 \le& L^2 \bar{C}_0  h^3.
    \end{cases}
\]
Now consider $p_1$, i.e., the deviation of mathematical expectation. By Ito's lemma, we have
\begin{align}
    & \bs{q}_2(h) - \bar{\bs{q}}_2(h) \nonumber \\
    =& -\alpha \int_0^h \nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1) ds \nonumber\\
    =& -\alpha \int_0^h \left[\int_0^s -\alpha \nabla^2 f(\bs{q}_2(r)) \nabla f(\bs{q}_2(r)) dr  +  \alpha \int_0^s \nabla \Delta f(\bs{q}_2(r)) dr + \rho \right]\, ds  \label{eq:tmp-2}
\end{align}
where $\rho$ is a stochastic integral term. Take expectation and norm for Equation \eqref{eq:tmp-2}, we have
\begin{align*}
    & \norm{\mathbb{E} \left[ \bs{q}_2(h) - \bar{\bs{q}}_2(h)   \right]} \\
    =& \alpha^2 \norm{\int_0^h \mathbb{E}\left[\int_0^s  \nabla^2 f(\bs{q}_2(r)) \nabla f(\bs{q}_2(r)) dr  - \int_0^s \nabla \Delta f(\bs{q}_2(r)) dr\right]\, ds} \\
    \le& \alpha^2 \int_0^h \mathbb{E}\left[\int_0^s  \|\nabla^2 f(\bs{q}_2(r))\|_2 \norm{\nabla f(\bs{q}_2(r))} dr  +  \int_0^s \norm{\nabla \Delta f(\bs{q}_2(r))} dr\right]\, ds \\
    \le& \alpha^2 \int_0^h \mathbb{E}\left[L  \int_0^s  \norm{\bs{q}_2(r)} dr  +  \int_0^s G (1 + \norm{\bs{q}_2(r)}) dr\right]\, ds \\
    =& \alpha^2 (L + G) \int_0^h  \int_0^s \mathbb{E}\norm{\bs{q}_2(r)}dr + \alpha^2 G\frac{h^2}{2} \\
    \le& \alpha^2 (L + G) \int_0^h  \int_0^s \mathbb{E}\norm{\bs{q}_2(r) - \bs{q}_1} + \mathbb{E} \norm{\bs{q}_1}dr + \alpha^2 G\frac{h^2}{2} \\
    \le& \alpha^2 (L + G) \int_0^h  \int_0^s \sqrt{\mathbb{E}\norm{\bs{q}_2(r) - \bs{q}_1}^2} + \mathbb{E} \norm{\bs{q}_1}dr + \alpha^2 G\frac{h^2}{2} \\
    \le& \alpha^2 (L + G) \sqrt{\bar{C}_0 h}  \frac{h^2}{2} + \alpha^2(L + G) \frac{h^2}{2} \mathbb{E}\norm{\bs{q}_1} + \alpha^2 G\frac{h^2}{2} \\
    \le& \alpha^2 \left\{ \frac{\sqrt{\bar{C}_0 h} + \mathbb{E}\norm{\bs{q}_1}}{2} (L+G) + \frac{G}{2} \right\} h^2 \\
    \le&  \frac{1}{2}\alpha^2   (L+G)  \left\{\sqrt{\bar{C}_0 h} + \mathbb{E}\norm{\bs{q}_1} + 1\right\} h^2
\end{align*}

Similarly, we have $\norm{\mathbb{E} \left[ \bs{p}_2(h) - \bar{\bs{p}}_2(h)   \right]} \le \frac{1}{2}\alpha   (L+G)  \left\{\sqrt{\bar{C}_0 h} + \mathbb{E}\norm{\bs{q}_1} + 1\right\} h^2$.

For $p_2$, i.e., mean-square error,  we have
\begin{align*}
    \mathbb{E} \norm{ \bs{q}_2(h) - \bar{\bs{q}}_2(h)}^2 \le& \alpha^2 \mathbb{E}\left\{\int_0^h \norm{\nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1)} ds\right\}^2 \\
    \le& \alpha^2 \mathbb{E} \left\{\int_0^h 1 ds \int_0^h  \norm{\nabla f(\bs{q}_2(s)) - \nabla f(\bs{q}_1)}^2 ds \right\} \\
    \le& \alpha^2 L^2 h \int_0^h  \mathbb{E}\norm{\bs{q}_2(s) - \bs{q}_1}^2 ds \\
    \le& \frac{\alpha^2 L^2 \bar{C}_0}{2}  h^3
\end{align*}
Similarly we obtain $\mathbb{E} \norm{ \bs{p}_2(h) - \bar{\bs{p}}_2(h)}^2 \le \frac{L^2 \bar{C}_0}{2} h^3$. Recall
\[
    \begin{cases}
        \bs{q}_3 - \bar{\bs{q}}_3 = \bs{q}_2(h) - \bar{\bs{q}}_2(h) + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} (\bs{p}_2(h) -\bar{\bs{p}}_2(h) ) \\
        \bs{p}_3 - \bar{\bs{p}}_3 = e^{-\gamma \frac{h}{2}} (\bs{p}_2(h) - \bar{\bs{p}}_2(h))
    \end{cases}.
\]
and it follows that when $0 < h \le \frac{1}{4L^\prime} < 1$
\begin{align}
    \norm{ \mathbb{E}\begin{bmatrix} \bs{q}_3 - \bar{\bs{q}}_3 \\ \bs{p}_3 - \bar{\bs{p}}_3 \end{bmatrix}} \le& \alpha (\alpha + 1 + \frac{h}{2}) (L+G) \frac{\sqrt{\bar{C}_0 h} + \mathbb{E}\norm{\bs{q}_1} + 1 }{2}  h^2 \label{eq:tool3}\\
    \mathbb{E} \norm{ \begin{bmatrix} \bs{q}_3 - \bar{\bs{q}}_3 \\ \bs{p}_3 - \bar{\bs{p}}_3 \end{bmatrix}}^2 \le& L^2 \bar{C}_0 \left( \alpha^2 + \frac{1}{2} + \frac{h^2}{4} \right) h^3  \label{eq:tool4}.
\end{align}
Finally we need to bound $\mathbb{E}\norm{\bs{q}_1}^2$ by $\mathbb{E}\norm{\bs{x}_0}^2$, to this end, we have
\begin{align}
    \mathbb{E}\norm{\bs{q}_1}^2 =& \mathbb{E}\norm{\bs{q}_0 + \frac{1 - e^{-\gamma \frac{h}{2}}}{\gamma} \bs{p}_0 + \sqrt{2\gamma} \int_0^\frac{h}{2} \frac{1 - e^{-\gamma(\frac{h}{2} - s)}}{\gamma} d\bs{B}(s)}^2 \nonumber\\
    \le& (1 + \frac{h^2}{4}) \mathbb{E} \norm{\bs{q}_0}^2 + (1 + \frac{h^2}{4}) \mathbb{E} \norm{\bs{p}_0}^2 +  2\gamma d \int_0^\frac{h}{2} \left(\frac{1 - e^{-\gamma (\frac{h}{2} - s)}}{\gamma} \right)^2 ds \nonumber\\
    \le& (1 + \frac{h^2}{4}) \mathbb{E} \norm{\bs{x}_0}^2 + \frac{\gamma d}{12}h^3 \\
    =& (1 + \frac{h^2}{4}) \norm{\bs{x}_0}^2 + \frac{\gamma d}{12}h^3  \label{eq:tool5}
\end{align}
Collecting all pieces together, including \eqref{eq:tool3}, \eqref{eq:tool4}, \eqref{eq:tool5}, the definition of $\bar{C}_0$ and $0 < h < \frac{1}{4L^\prime}$, it is not difficult to obtain the following
\begin{align*}
    \norm{ \mathbb{E}\begin{bmatrix} \bs{q}_3 - \bar{\bs{q}}_3 \\ \bs{p}_3 - \bar{\bs{p}}_3 \end{bmatrix}} \le& \bar{C}_1  h^2 \\
    \left(\mathbb{E} \norm{ \begin{bmatrix} \bs{q}_3 - \bar{\bs{q}}_3 \\ \bs{p}_3 - \bar{\bs{p}}_3 \end{bmatrix}}^2\right)^\frac{1}{2} \le& \bar{C}_2 h^\frac{3}{2}
\end{align*}
with
\[
    \bar{C}_1 = \alpha (\alpha + 1.125) (L+G) \left[ 0.5 + 0.71\norm{\bs{x}_0} + (1.14\sqrt{\alpha} + 0.21 \sqrt{\gamma} h) \sqrt{hd} \right]
\]
and
\[
\bar{C}_2 = L(\alpha + 0.73) \left( 2.30\sqrt{h} \alpha L \norm{\bs{x}_0} + (2.27\sqrt{\alpha} + 0.12 \sqrt{\gamma} h )\sqrt{d} \right)
\]

\end{proof}

\subsection{Local error between Algorithm \ref{alg:HFHR} and HFHR dynamics}
\begin{lemma}\label{lemma:local}
Assume $f$ is $L$-smooth, $\bs{0} \in \argmin_{\bs{x} \in \mathbb{R}^d} f(\bs{x})$, i.e. $\nabla f(\bs{0}) = \bs{0}$ and the operator $\nabla \Delta f$ grows at most linearly, i.e. $\norm{\nabla \Delta f(\bs{q}) }\le G\sqrt{1 + \norm{\bs{q}}^2 }$. If $0 < h \le \frac{1}{4 L^\prime}$, then compared with the HFHR dynamics, the implementable Strang's splitting method has local weak error  of order $p_1 = 2$ and local mean-squared error of order $p_2 = 1.5$, i.e. there exist constants $C_1, C_2 > 0$ such that
\[
    \norm{ \mathbb{E} \bs{x}(h)  - \mathbb{E} \bar{\bs{x}}(h) } \le C_1 h^{p_1}
\]
\[
    \left(\mathbb{E} \left[ \norm{ \bs{x}(h) - \bar{\bs{x}}(h) }^2 \right]\right)^\frac{1}{2} \le C_2 h^{p_2}
\]
where $\bs{x}(h) = \begin{bmatrix} \bs{q}(h) \\ \bs{p}(h) \end{bmatrix}$ is the solution of HFHR with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$ and $\bar{\bs{x}}(h) = \begin{bmatrix} \bar{\bs{q}}(h) \\ \bar{\bs{p}}(h) \end{bmatrix}$ is the solution of the implementable Strang's splitting with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$, $p_1 = 2$ and $p_2 = 1.5$. More concretely, we have
\[
    C_1 = (L+G) \max\{\alpha + 1.25, \gamma + 1\} \left[ 0.5\alpha + (1.74 + 0.71\alpha)\norm{\bs{x}_0} + \left( 1.26\sqrt{\alpha} + 1.14\alpha \sqrt{\alpha} + 2.32\sqrt{\gamma}\right) \sqrt{hd} \right]
\]
and
\[
    C_2 = L \max\{\alpha + 1.25, \gamma + 1\} \left[ (1.92 + 2.30\alpha L) \sqrt{h} \norm{\bs{x}_0} + (2.60\sqrt{\alpha} + 3.34\sqrt{\gamma} h) \sqrt{d} \right]
\]
\end{lemma}
\begin{proof}
Denote by $\hat{\bs{x}}(h) = \begin{bmatrix} \hat{\bs{q}}(h) \\ \hat{\bs{p}}(h) \end{bmatrix}$ the solution of the exact Strang's splitting method with initial value $\bs{x}_0 = \begin{bmatrix} \bs{q}_0 \\ \bs{p}_0 \end{bmatrix}$. By triangle inequality and Minkowski's inequality, we have
\begin{align*}
    \norm{\mathbb{E} \bs{x}(h) - \mathbb{E}\bar{\bs{x}}(h)} \le& \norm{\mathbb{E} \bs{x}(h) - \mathbb{E}\hat{\bs{x}}(h)} + \norm{\mathbb{E} \hat{\bs{x}}(h) - \mathbb{E}\bar{\bs{x}}(h)}, \\
    \left(  \mathbb{E} \norm{\bs{x}(h) - \bar{\bs{x}}(h)}^2 \right)^\frac{1}{2} \le& \left(  \mathbb{E} \norm{\bs{x}(h) - \hat{\bs{x}}(h)}^2 \right)^\frac{1}{2} + \left(  \mathbb{E} \norm{\hat{\bs{x}}(h) - \bar{\bs{x}}(h)}^2 \right)^\frac{1}{2}.
\end{align*}
By Lemma \ref{lemma:local-exact} and \ref{lemma:local-implementable}, we have
\begin{align*}
     \norm{\mathbb{E} \bs{x}(h) - \mathbb{E}\hat{\bs{x}}(h)} \le \widehat{C}_1 h^2, &\quad \norm{\mathbb{E}  \hat{\bs{x}}(h) - \mathbb{E}\bar{\bs{x}}(h)} \le \bar{C}_1 h^2 \\
     \left(  \mathbb{E} \norm{\bs{x}(h) - \hat{\bs{x}}(h)}^2 \right)^\frac{1}{2} \le \widehat{C}_2 h^\frac{3}{2}, &\quad \left(  \mathbb{E} \norm{\hat{\bs{x}}(h) - \bar{\bs{x}}(h)}^2 \right)^\frac{1}{2} \le \bar{C}_2 h^\frac{3}{2}
\end{align*}
and hence
\begin{align*}
    \norm{\mathbb{E} \bs{x}(h) - \mathbb{E}\bar{\bs{x}}(h)} \le& (\widehat{C}_1 + \bar{C}_1) h^2 \\
    \left(  \mathbb{E} \norm{\bs{x}(h) - \bar{\bs{x}}(h)}^2 \right)^\frac{1}{2} \le& (\widehat{C}_2 + \bar{C}_2) h^\frac{3}{2}
\end{align*}
with
\begin{align*}
    \widehat{C}_1 + \bar{C}_1 \le&  C_1  \\
    \triangleq& (L+G) \max\{\alpha + 1.25, \gamma + 1\} \left[ 0.5\alpha + (1.74 + 0.71\alpha)\norm{\bs{x}_0} + \left( 1.26\sqrt{\alpha} + 1.14\alpha \sqrt{\alpha} + 2.32\sqrt{\gamma}\right) \sqrt{hd} \right] \\
    \widehat{C}_2 + \bar{C}_2 \le& C_2 \triangleq L \max\{\alpha + 1.25, \gamma + 1\} \left[ (1.92 + 2.30\alpha L) \sqrt{h} \norm{\bs{x}_0} + (2.60\sqrt{\alpha} + 3.34\sqrt{\gamma} h) \sqrt{d} \right]
\end{align*}

\end{proof}

\section{$\alpha$ does create acceleration even after discretization: an analytical demonstration}
\label{sec_onAlpha}
If $\alpha \to \infty$ while $\gamma$ remains fixed, then $dq=-\alpha \nabla f(q)+\sqrt{2\alpha}dW$ is the dominant part of the dynamics, and in this case the role of $\alpha$ could be intuitively understood as to simply rescale the time of gradient flow, which does not create any algorithmic advantage, as the timestep of discretization has to scale like $1/\alpha$ in this case. However, finite $\alpha$ no longer corresponds to solely a time-scaling, but closely couples with the dynamics and creates acceleration. This is true even after the continuous dynamics is discretized by an algorithm .

We will analytically illustrate this point by considering quadratic $f$. In this case, the diffusion process remains Gaussian, and it suffices to quantify the convergence of its mean and covariance. In fact, it can be shown that both have the same speed of convergence, and therefore for simplicity we will only consider the mean process. Two demonstrations (with different focuses) will be provided.

\paragraph{Demonstration 1 (1D, $\gamma$ given; infinite acceleration).}
Consider $f(x)=x^2/2$, $\gamma$ fixed. The mean process is
\[ \begin{cases}
    \dot{q} &= p-\alpha q \\
    \dot{p} &= -q-\gamma p
\end{cases} \]
Consider, for simplicity, an Euler-Maruyama discretization of the HFHR dynamics, which coressponds to a Forward Euler discretization of the mean process (other numerical methods can be analyzed analogously):
\[ \begin{bmatrix} q_{k+1} \\ p_{k+1} \end{bmatrix} =
A \begin{bmatrix} q_k \\ p_k \end{bmatrix}, \qquad
A=\begin{bmatrix} 1-\alpha h & h \\ -h & 1-\gamma h \end{bmatrix}.
\]
We will show that, unless $\gamma=2$, an appropriately chosen $\alpha$ will converge infinitely faster than the case with $\alpha=0$, if both cases use the optimal $h$.

To do so, let us compute $A$'s eigenvalues, which are
\[
\frac{1}{2} \left( 2 - (\alpha+\gamma)h \pm h \sqrt{-4+(\alpha-\gamma)^2} \right)
\]
Consider the case where $|\alpha-\gamma|\leq 2$, then the eigenvalues are a pair of complex conjugates. Their modulus
determines the speed of convergence, and it can be computed to be
\[
    \frac{1}{2}\sqrt{(2-(\alpha+\gamma)h)^2+h^2(4-(\alpha-\gamma)^2)} = \sqrt{1-(\alpha+\gamma)h+(1+\alpha\gamma)h^2}
\]
Minimizing the quadratic function gives the optimal $h$ that ensures the fastest speed of convergence, and the optimal $h$ is
\[
  h=\frac{\alpha+\gamma}{2(1+\alpha\gamma)}
\]
and the optimal spectral radius is
\[
    \sqrt{1-\frac{(\alpha+\gamma)^2}{4(1+\alpha\gamma)}}.
\]
When one uses low-resolution ODE, in which $\alpha=0$, the optimal rate is $1-\gamma^2/4$ (note it is not surprising that the critically damped case, i.e., $\gamma=2$, will give the fastest convergence).

If $\gamma\neq 2$, the additional introduction of $\alpha$ can accelerate the convergence by reducing the spectral radius. For instance, if $\alpha=\gamma+2$, upon choosing the optimal $h=\frac{1}{1+\gamma}$, the optimal spectral radius is 0 (note in this case $A$ actually has Jordan canonical form of $\begin{bmatrix} 0 & 1 \\ 0 & 0 \end{bmatrix}$ and thus the discretization converges in 2 steps instead of 1, irrespective of the initial condition).

\paragraph{Demonstration 2 (multi-dim, $\gamma$, $\alpha$ and $h$ all to be chosen; acceleration quantified in terms of condition number).}

Consider quadratic $f$ with positive definite Hessian, whose eigenvalues are $1=\lambda_1 < \cdots < \lambda_n = \epsilon^{-1}$ for some $0 < \epsilon \ll 1$. Assume without loss of generality that $f=q_1^2/2+\epsilon^{-1} q_2^2/2$. Similar to Demonstration 1, the forward Euler discretization of the mean process is
\begin{equation} \begin{bmatrix} q_{1,k+1} \\ p_{1, k+1} \\ q_{2,k+1} \\ p_{2, k+1} \end{bmatrix} =
\begin{bmatrix} A_1 & 0 \\ 0 & A_2 \end{bmatrix} \begin{bmatrix} q_{1,k} \\ p_{1, k} \\ q_{2,k} \\ p_{2, k} \end{bmatrix}, \quad
A_1=\begin{bmatrix} 1-\alpha h & h \\ -h & 1-\gamma h \end{bmatrix}, \quad
A_2=\begin{bmatrix} 1-\alpha \epsilon^{-1} h & h \\ -\epsilon^{-1} h & 1-\gamma h \end{bmatrix}
\label{eq_conditionNumberSec_iter}
\end{equation}
We will (i) find $h$ and $\gamma$ that lead to fastest convergence of the ULD discretization, i.e. the above iteration with $\alpha=0$, and then (ii) constructively show the existence of $h$, $\gamma$ and $\alpha$ that lead to faster convergence than the optimal one in (i) --- note these may not even be the optimal choices for HFHR, but they already lead to significant acceleration. More specifically,

\textbf{\underline{(i)}} In a ULD setup, $\alpha=0$. It can be computed that the eigenvalues of $A_1$ and $A_2$ are respectively
\[
\frac{1}{2} \left( 2-h\gamma\pm h\sqrt{-4+\gamma^2} \right)   \qquad\text{and}\qquad
\frac{1}{2} \left( 2-h\gamma\pm h\sqrt{-4\epsilon^{-1}+\gamma^2} \right)
\]
We now seek $\gamma>0,h>0$ to minimize the maximum of their norms for obtaining the optimal convergence rate. This is done in cases.

\noindent \underline{Case (i1)} When $\gamma \leq 2$, both $A_1$ and $A_2$ eigenvalues are complex conjugate pairs. To minimize the maximum of their norms, let's first see if their norms could be made equal.

$A_1$ eigenvalue's norm squared $\times 4$ is
\begin{equation}
    (2-h\gamma)^2-h^2(-4+\gamma^2) = 4(h-\gamma/2)^2+4-\gamma^2
    \label{ugv9317ughoiuhoig1314}
\end{equation}
$A_2$ eigenvalue's norm squared $\times 4$ is
\begin{equation}
    (2-h\gamma)^2-h^2(-4\epsilon^{-1}+\gamma^2) = 4\epsilon^{-1}(h-\epsilon\gamma/2)^2+4-\epsilon\gamma^2
    \label{qbogoyiuhb31ogyiubo1}
\end{equation}
It can be seen that for \eqref{ugv9317ughoiuhoig1314} is always strictly smaller than \eqref{qbogoyiuhb31ogyiubo1} for any $h>0$. Therefore, the max of the two is minimized when $h=\epsilon \gamma/2$, and the corresponding max value is $4-\epsilon \gamma^2$. $\gamma$ that minimizes this max value is $\gamma=2$. Corresponding rate of convergence is
\[
    \sqrt{1-\epsilon}.
\]

\noindent \underline{Case (i2)} When $\gamma \geq 2\epsilon^{-1/2}$, both $A_1$ and $A_2$ eigenvalues are real. Since $\epsilon\ll 1$, we can order them$\times 2$ as
\[
    2-h\gamma-h\sqrt{-4+\gamma^2} <
    2-h\gamma-h\sqrt{-4\epsilon^{-1}+\gamma^2} <
    2-h\gamma+h\sqrt{-4\epsilon^{-1}+\gamma^2} <
    2-h\gamma+h\sqrt{-4+\gamma^2} < 2.
\]
To minimize the max of their norms, consider cases in which the smallest of four is negative, in which case at optimum one should have
\[
    -(2-h\gamma-h\sqrt{-4+\gamma^2}) = 2-h\gamma+h\sqrt{-4+\gamma^2}.
\]
This gives $h=2/\gamma$ (which does verify the assumption that the smallest of four is negative). Corresponding max of their norms is thus $\sqrt{1-4/\gamma^2}$. $\gamma$ that minimizes this max value is $\gamma=2\epsilon^{-1/2}$, which gives rate of convergence of
\[
    \sqrt{1-\epsilon}.
\]

\noindent \underline{Case (i3)} When $2 \leq \gamma \leq 2\epsilon^{-1/2}$, $A_1$ eigenvalues are real and $A_2$ eigenvalues are complex conjugates. Again, the max of their norms is minimized if the norms can be made all equal.

Note $A_1$ eigenvalues cannot be of the same sign, because otherwise $2-h\gamma-h\sqrt{-4+\gamma^2}=2-h\gamma+h\sqrt{-4+\gamma^2}$, which means either $h=0$ or $\gamma=2$, but if $\gamma=2$ then $2-h\gamma+h\sqrt{-4+\gamma^2}$ being equal to 2*norm of $A_2$ eigenvalue, which is $\sqrt{4\epsilon^{-1}(h-\epsilon\gamma/2)^2+4-\epsilon\gamma^2}$, leads to $h=0$ again.

Therefore, the equality of norms of $A_1$, $A_2$ eigenvalues means
\[
  -(2-h\gamma-h\sqrt{-4+\gamma^2})
  =2-h\gamma+h\sqrt{-4+\gamma^2}
  =\sqrt{4\epsilon^{-1}(h-\epsilon\gamma/2)^2+4-\epsilon\gamma^2}.
\]
The first equality gives $h\gamma=2$, which, together with the second equality, gives $h=\pm\sqrt{\frac{2\epsilon}{1+\epsilon}}$. Selecting the positive value of optimal $h$, we also obtain optimal $\gamma=\sqrt{2(1+\epsilon)}\epsilon^{-1/2}$, which is $\leq 2 \epsilon^{-1/2}$ and thus satisfying our assumption ($2 \leq \gamma \leq 2\epsilon^{-1/2}$). The corresponding rate of convergence is thus
\[
    \frac{1}{2} \left(2-h\gamma+h\sqrt{-4+\gamma^2}\right) = \sqrt{\frac{1-\epsilon}{1+\epsilon}}.
\]

\noindent \underline{Summary of (i)} Since $\sqrt{\frac{1-\epsilon}{1+\epsilon}} < \sqrt{1-\epsilon}$, the ULD Euler-Maruyama discretization converges the fastest when
\[
    h=\sqrt{\frac{2\epsilon}{1+\epsilon}}, \qquad
    \gamma=\sqrt{2(1+\epsilon)}\epsilon^{-1/2},
\]
and the corresponding discount factor of convergence is
\begin{equation}
    \sqrt{\frac{1-\epsilon}{1+\epsilon}},   \qquad \text{where $\epsilon=1/\kappa$ with $\kappa$ being Hessian's condition number.}
    \label{eq_ULDoptimalRate_conditionNum}
\end{equation}

\textbf{\underline{(ii)}} Now consider the HFHR setup. Let's first state a result: when
\begin{align}
    & \gamma=\frac{\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+\epsilon +3}{2 c \epsilon ^2+2 c \epsilon} > 0, \\
    & \alpha=\frac{-\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+3 \epsilon +1}{2 c \epsilon ^2+2 c \epsilon} > 0, \qquad
    h=c \epsilon
    \label{eq_HFHRnonoptimalParameters}
\end{align}
for any $c>0$ independent of $\epsilon$, the iteration \eqref{eq_conditionNumberSec_iter} converges with discount factor
\begin{equation}
    \frac{1}{\sqrt{2} (1+\epsilon)} \sqrt{(1-\epsilon) \left(1-\epsilon+\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+\left(4 c^2+1\right) \epsilon ^2-2 \epsilon +1}\right)}.
    \label{eq_HFHRnonoptimalRate_conditionNum}
\end{equation}
While the exact expression is lengthy, it can proved that the HFHR non-optimal discount factor \eqref{eq_HFHRnonoptimalRate_conditionNum} is strictly smaller than the ULD optimal discount factor \ref{eq_ULDoptimalRate_conditionNum} for not only small but also large $\epsilon$'s.

For some quantitative intuition, discount factors Taylor expanded in $\epsilon$ are respectively
\begin{align}
    & \text{HFHR non-optimal:}   &&
    1-2 \epsilon +\left(\frac{c^2}{2}+2\right) \epsilon ^2 +\mathcal{O}\left(\epsilon ^3\right) \qquad\qquad \\
    & \text{ULD optimal:}     && 1-\epsilon+\frac{\epsilon^2}{2} +\mathcal{O}\left(\epsilon ^3\right) \qquad\qquad
\end{align}

\begin{figure}[h]
    \vspace{-120pt}
    \centering
    \includegraphics[width=0.8\textwidth]{img/discountFactors.pdf}
    \vspace{-120pt}
    \caption{Acceleration of HFHR algorithm over ULD algorithm (despite of an additional constraint $\alpha$ may place on $h$) for multi-dimensional quadratic objectives. $1/\epsilon$ is the condition number.}
    \label{fig_discountFactors}
\end{figure}
The exact expressions of discount factors are also plotted  in Fig.\ref{fig_discountFactors} ($c=1$ was arbitrarily chosen) and one can see acceleration for any (not necessarily small) $\epsilon$.

\textbf{\underline{(ii details)}} How were values in \eqref{eq_HFHRnonoptimalParameters} chosen? Following the idea detailed in (i), we consider a case where $A_1$ eigenvalues are both real, $A_2$ eigenvalues are complex conjugates, and all their norms are equal. Note there are 3 more cases, namely real/real, complex/real, and complex/complex, but we do not optimize over all cases for simplicity --- the real/complex case is enough for outperforming the optimal ULD.

This case leads to at least the following equations
\begin{equation} \begin{cases}
    \text{tr} A_1 &= 0 \\
    \det A_1 + \det A_2 &= 0
\end{cases}
\label{agijqhrogubvqorehgb97413hbt8gob}
\end{equation}
One can solve this system of equations to obtain $\alpha$ and $\gamma$ as functions of $h$. Following the idea of choosing $h$ small enough to resolve the stiffness of the ODE
\[\begin{cases}
    \dot{q}_2 &= p_2 - \alpha \epsilon^{-1} q_2 \\
    \dot{p}_2 &= -\epsilon^{-1} q_2 - \gamma p_2
\end{cases},\]
pick $h=c \epsilon$. Then \eqref{agijqhrogubvqorehgb97413hbt8gob} gives
\begin{align*}
    & \gamma=\frac{\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+\epsilon +3}{2 c \epsilon ^2+2 c \epsilon} \\
    & \alpha=\frac{-\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+3 \epsilon +1}{2 c \epsilon ^2+2 c \epsilon}
\end{align*}
or
\begin{align*}
    & \gamma = \frac{-\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+\epsilon +3}{2 c \epsilon ^2+2 c \epsilon } \\
    & \alpha = \frac{\sqrt{4 c^2 \epsilon ^4+8 c^2 \epsilon ^3+4 c^2 \epsilon ^2+\epsilon ^2-2 \epsilon +1}+3 \epsilon +1}{2 c \epsilon ^2+2 c \epsilon }
\end{align*}
The former is our choice \eqref{eq_HFHRnonoptimalParameters} because it can be checked that the latter leads to $\det A_1 > 0$ which violates the assumption of a pair of plus and minus real eigenvalues.

It is possible to find optimal $\alpha,\gamma,h$ for HFHR for the Gaussian cases. One has to minimize $\det A_2$ under the constraint $\det A_2 > 0$ in addition to \eqref{agijqhrogubvqorehgb97413hbt8gob}. And then do similar calculations for the other 3 cases, and then finally the best among the 4 cases. Doing so however does not give enough insights to determine optimal hyperparameters for sampling general distributions.

\end{document}

config: arxiv

file: 2006.09230/(JMLR) Hessian-Free High-Resolution Nesterov Acceleration for Sampling/sample.tex