\section{Means and covariance matrices of $\pi_{\scriptscriptstyle{\operatorname{LMC}}}, \pi_{\scriptscriptstyle{\operatorname{FP}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ in the Bayesian linear regression}
\label{sec:mean_var_gaussian}
In this Section, we provide explicit expressions of the covariance matrices of $\pi_{\scriptscriptstyle{\operatorname{LMC}}},\pi_{\scriptscriptstyle{\operatorname{FP}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ in the context of the Bayesian linear regression. In this setting, the algorithms are without bias, \textit{i.e.}
\begin{equation}\label{eq:regression-without-bias}
  \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{LMC}}}(\mathrm{d} \theta) = \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{FP}}}(\mathrm{d} \theta) = \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \theta) = \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{SGD}}}(\mathrm{d} \theta) = \int_{\mathbb{R}^d} \theta \pi(\mathrm{d} \theta) = \theta^{\star} \;.
\end{equation}
Before giving the expressions of the variances in \Cref{thm:variance-piula-pisgld-pisgd-pifp}, we define $\operatorname{T}:\mathbb{R}^{d\times d}\to \mathbb{R}^{d\times d}$ for all $A\in\mathbb{R}^{d\times d}$ by
\begin{equation}\label{eq:def-operator-T}
  \operatorname{T} (A) = \expe{\parenthese{\frac{\operatorname{Id}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS} x_i x_i^{\operatorname{T}} - \Sigma}^{\otimes 2} A} = \frac{N}{p}\sum_{i=1}^{N} \parenthese{\frac{x_i x_i^{\operatorname{T}}}{\sigma^2_y} + \frac{\operatorname{Id}}{N\sigma^2_{\tp}} - \frac{\Sigma}{N}}^{\otimes 2} A \;,
\end{equation}
where $S$ is a random subsample of $\defEns{1,\ldots,N}$ with replacement of size $p\in\mathbb{N}^*$.
Note that, in this setting, $\tilde{L} = \max_{i\in\defEns{1,\ldots,N}} \norm[2]{x_i}$ and $m$ is the smallest eigenvalue of $\Sigma$.
There exists $r \in\ccintLigne{0,L/(\sqrt{p}m)}$ such that
\begin{equation}\label{eq:def-rSig}
  \operatorname{T} \preceq r^2 \Sigma^{\otimes 2}
\end{equation}
\textit{i.e.}~for all $A\in\mathbb{R}^{d\times d}$, $\operatorname{Tr}(A^{\operatorname{T}} \operatorname{T} \cdot A) \leq r^2 \operatorname{Tr}(A^{\operatorname{T}} \Sigma^{\otimes 2} A)$. Assuming that $\liminf_{N\to+\infty} N^{-1} m >0$, $r$ can be chosen independently of $N$.

\begin{theorem}\label{thm:variance-piula-pisgld-pisgd-pifp}
  Consider the case of the Bayesian linear regression.
  We have for all $\gamma\in\ooint{0,2/L}$
  \begin{equation*}
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{LMC}}}(\mathrm{d} \theta) = \parenthese{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma \Sigma \otimes \Sigma}^{-1} (2\operatorname{Id})\;,
  \end{equation*}
  and for all $\gamma\in\ooint{0,2/\{(1+r^2)L\}}$,
  \begin{align*}
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{FP}}}(\mathrm{d} \theta) &= \defEns{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})}^{-1} (2\operatorname{Id})\;, \\
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \theta) &= \defEns{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})}^{-1} \\
    &\phantom{-----}\cdot \defEns{2\operatorname{Id} + \frac{\gamma N}{p}\sum_{i=1}^{N} \parenthese{\frac{(x_i^{\operatorname{T}} \theta^{\star} - y_i)x_i}{\sigma^2_y} + \frac{\theta^{\star}}{\sigma^2_{\tp}}}^{\otimes 2}} \;, \\
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{SGD}}}(\mathrm{d} \theta) &= \defEns{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})}^{-1} \\
    &\phantom{---------}\cdot \frac{\gamma N}{p}\sum_{i=1}^{N} \parenthese{\frac{(x_i^{\operatorname{T}} \theta^{\star} - y_i)x_i}{\sigma^2_y} + \frac{\theta^{\star}}{\sigma^2_{\tp}}}^{\otimes 2} \;.
  \end{align*}
\end{theorem}

\begin{proof}
We prove the result for SGLD, the adaptation to the other algorithms is immediate. Let $\gamma\in\ooint{0,2/\{(1+r^2)L\}}$, $\theta_0$ be distributed according to $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\theta_1$ be given by \eqref{eq:SGLD}. By definition of $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$, $\theta_1$ is distributed according to $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$. We have
\begin{multline*}
  \expe{(\theta_1 - \theta^{\star})^{\otimes 2}} = \mathbb{E} \Bigg[ \Bigg[ \defEns{\operatorname{Id} - \gamma\parenthese{\frac{\operatorname{Id}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS_1} x_i x_i^{\operatorname{T}}}} (\theta_0 - \theta^{\star}) \\
  - \gamma\parenthese{\frac{\theta^{\star}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS_1} (x_i^{\operatorname{T}} \theta^{\star} - y_i) x_i} + \sqrt{2\gamma}Z_1 \Bigg]^{\otimes 2} \Bigg] \;.
\end{multline*}
Using that $\theta_0,S_1,Z_1$ are mutually independent, we obtain
\begin{multline*}
  \defEns{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma \expe{\parenthese{\frac{\operatorname{Id}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS_1} x_i x_i^{\operatorname{T}}}^{\otimes 2}}} \expe{(\theta_0 - \theta^{\star})^{\otimes 2}} \\
  = 2\operatorname{Id} + \gamma \expe{\parenthese{\frac{\theta^{\star}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS_1} (x_i^{\operatorname{T}} \theta^{\star} - y_i) x_i}^{\otimes 2}}
\end{multline*}
and
\begin{multline*}
  \defEns{\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})} \expe{(\theta_0 - \theta^{\star})^{\otimes 2}} \\
  = 2\operatorname{Id} + \frac{\gamma N}{p}\sum_{i=1}^{N} \parenthese{\frac{(x_i^{\operatorname{T}} \theta^{\star} - y_i)x_i}{\sigma^2_y} + \frac{\theta^{\star}}{\sigma^2_{\tp}}}^{\otimes 2} \;.
\end{multline*}
On $\mathbb{R}^{d\times d}$ equipped with the Hilbert-Schmidt inner product, $\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})$ is a positive definite operator. Indeed, by \eqref{eq:def-rSig},
\begin{align*}
  \operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T}) & \succeq \operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (1+r^2) \Sigma^{\otimes 2} \\
  & = \parenthese{\operatorname{Id} - \gamma \frac{1+r^2}{2} \Sigma} \otimes \Sigma + \Sigma \otimes \parenthese{\operatorname{Id} - \gamma \frac{1+r^2}{2} \Sigma} \succ 0
\end{align*}
for $\gamma\in\ooint{0,2/\{(1+r^2)L\}}$. $\operatorname{Id} \otimes \Sigma + \Sigma \otimes \operatorname{Id} - \gamma (\Sigma^{\otimes 2} + \operatorname{T})$ is thus invertible, which concludes the proof.
\end{proof}

The covariance matrices make clearly visible the different origins of the noise. The Gaussian noise is responsible of the term $2\operatorname{Id}$, while the multiplicative and additive parts of the stochastic gradient (see \eqref{eq:def-partm-parta}) are related to the operator $\operatorname{T}$ and to the term
\begin{equation}\label{eq:noise-additive-part}
  \frac{\gamma N}{p}\sum_{i=1}^{N} \parenthese{\frac{(x_i^{\operatorname{T}} \theta^{\star} - y_i)x_i}{\sigma^2_y} + \frac{\theta^{\star}}{\sigma^2_{\tp}}}^{\otimes 2}
\end{equation}
respectively.

Denote by
\begin{equation}\label{eq:def-eta-1}
  \eta_1 = \inf_{N \geq 1} \defEns{\frac{2N}{L} \wedge \frac{2N}{(1+r^2) L}} > 0 \;.
\end{equation}

\begin{corollary}\label{corollary:variances-piula-pisgld-pisgd-pifp}
  Consider the case of the Bayesian linear regression. Set $\gamma = \eta/N$ with $\eta\in\oointLigne{0,\eta_1}$ and assume that $\liminf_{N\to+\infty} N^{-1} m >0$.
  \begin{align*}
    & \int_{\mathbb{R}^d} \norm[2]{\theta-\theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{LMC}}}(\mathrm{d} \theta) = d\Theta_{N\to+\infty}(N^{-1}) \;, \;
    \int_{\mathbb{R}^d} \norm[2]{\theta-\theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{FP}}}(\mathrm{d} \theta) = d\Theta_{N\to+\infty}(N^{-1}) \;, \\
    & \int_{\mathbb{R}^d} \norm[2]{\theta-\theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \theta) = \eta d\Theta_{N\to+\infty}(1) \;,\;
    \int_{\mathbb{R}^d} \norm[2]{\theta-\theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{SGD}}}(\mathrm{d} \theta) = \eta d\Theta_{N\to+\infty}(1) \;.
  \end{align*}
\end{corollary}

Recall that, according to the Bernstein-von Mises theorem, the variance of $\pi$ is of the order $d/N$ when $N$ is large. The corollary confirms that $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ is very far from $\pi$ when the constant step size $\gamma$ is chosen proportional to $1/N$.

\section{Illustration of \Cref{lemma:mean-pi-tpstar}, \Cref{thm:bias-variance-tpstar-ula-fp} and \Cref{thm:bias-variance-tpstar-sgld-sgd}}
\label{sec:illustration-m1-m2}

We provide in \Cref{fig:illustration-m1-m2} an illustration of the results of \Cref{subsec:variance-piula-pifp-puslgd-pisgd} as the number of data items $N$ goes to infinity.

\begin{figure}
\centering
  \begin{tikzpicture}[scale=1.5]
    \draw (0,0) node[below]{$\bar{\theta}$} ;
    \draw (0,0) node{$\bullet$} ;
    \coordinate (LMC) at (0.5, 0.25);
    \draw [dashed, color=red] (LMC) -- ++(90:1) node[midway, right]{$1/N$};
    \draw [color=red] (LMC) node[right]{$\bar{\theta}_{\scriptscriptstyle{\operatorname{LMC}}}$} ;
    \draw [color=red] (LMC) node{$\bullet$} ;
    \draw [color=red] (LMC) circle (1) ;
    \coordinate (FP) at (-0.7, 0);
    \draw [dashed, color=blue] (FP) -- ++(90:1) node[midway, left]{$1/N$};
    \draw [color=blue] (FP) node[right]{$\bar{\theta}_{\scriptscriptstyle{\operatorname{FP}}}$} ;
    \draw [color=blue] (FP) node{$\bullet$} ;
    \draw [color=blue] (FP) circle (1) ;
    \coordinate (SGD) at (3, 0);
    \draw [color=green] (SGD) node{$\bullet$} ;
    \draw [color=green] (SGD) node[right]{$\bar{\theta}_{\scriptscriptstyle{\operatorname{SGD}}}$} ;
    \draw [dashed, color=green] (0,0) -- (SGD) node[midway, below]{$\sim 1$};
    \draw [color=green] (SGD) circle (3) ;
    \coordinate (SGLD) at (3.2, 0.5);
    \draw [color=magenta] (SGLD) node{$\bullet$} ;
    \draw [color=magenta] (SGLD) node[above]{$\bar{\theta}_{\scriptscriptstyle{\operatorname{SGLD}}}$} ;
    \draw [dashed] (SGLD) -- (SGD) node[midway, right]{$1/N$};
    \draw [color=magenta] (SGLD) circle (3) ;
  \end{tikzpicture}
\caption{\label{fig:illustration-m1-m2} Illustration of \Cref{lemma:mean-pi-tpstar}, \Cref{thm:bias-variance-tpstar-ula-fp} and \Cref{thm:bias-variance-tpstar-sgld-sgd} in the asymptotic $N\to+\infty$. $\bar{\theta}$, $\bar{\theta}_{\scriptscriptstyle{\operatorname{SGD}}}$, $\bar{\theta}_{\scriptscriptstyle{\operatorname{LMC}}}$, $\bar{\theta}_{\scriptscriptstyle{\operatorname{FP}}}$ and $\bar{\theta}_{\scriptscriptstyle{\operatorname{SGLD}}}$ are the means under the stationary distributions $\pi$, $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$, $\pi_{\scriptscriptstyle{\operatorname{LMC}}}$, $\pi_{\scriptscriptstyle{\operatorname{FP}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$, respectively. The associated circles indicate the order of magnitude of the covariance matrix. While LMC and SGLDFP concentrate to the posterior mean $\bar{\theta}$ with a covariance matrix of the order $1/N$, SGLD and SGD are at a distance of order $\sim 1$ of $\bar{\theta}$ and do not concentrate as $N\to+\infty$.}
\end{figure}

\section{Numerical experiments}
\label{sec:numerical-experiments}

\paragraph*{Simulated data} For illustrative purposes, we consider a Bayesian logistic regression in dimension $d=2$. We simulate $N=10^5$ covariates $\{x_i\}_{i=1}^{N}$ drawn from a standard $2$-dimensional Gaussian distribution and we denote by $\mathbf{X} \in\mathbb{R}^{N\times d}$ the matrix of covariates such that the $i^{\text{th}}$ row of $\mathbf{X}$ is $x_i$.
Our Bayesian regression model is specified by a Gaussian prior of mean $0$ and covariance matrix the identity, and a likelihood given for $y_i\in\{0,1\}$ by $p(y_i |x_i, \theta) = (1 + \mathrm{e}^{-x_i^{\operatorname{T}} \theta})^{-y_i} (1 + \mathrm{e}^{x_i^{\operatorname{T}} \theta})^{y_i - 1}$. We simulate $N$ observations $\{y_i\}_{i=1}^{N}$ under this model. In this setting, \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex} are satisfied, and \Cref{assum:strongly-convex} holds if the state space is compact.

To illustrate the results of \Cref{subsec:variance-piula-pifp-puslgd-pisgd}, we consider $10$ regularly spaced values of $N$ between $10^2$ and $10^5$ and we truncate the dataset accordingly. We compute an estimator $\hat{\theta}$ of $\theta^{\star}$ using SGD \cite{scikit-learn} combined with the BFGS algorithm \cite{scipy}. For the LMC, SGLDFP, SGLD and SGD algorithms, the step size $\gamma$ is set equal to $(1+\delta/4)^{-1}$ where $\delta$ is the largest eigenvalue of $\mathbf{X}^{\operatorname{T}} \mathbf{X}$. We start the algorithms at $\theta_0 = \hat{\theta}$ and run $n=1/\gamma$ iterations where the first $10\%$ samples are discarded as a burn-in period.

We estimate the means and covariance matrices of $\pi_{\scriptscriptstyle{\operatorname{LMC}}}, \pi_{\scriptscriptstyle{\operatorname{FP}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ by their empirical averages $\bar{\theta}_n = (1/n)\sum_{k=0}^{n-1} \theta_k$ and $\{1/(n-1)\}\sum_{k=0}^{n-1} (\theta_k - \bar{\theta}_n)^{\otimes 2}$. We plot the mean and the trace of the covariance matrices for the different algorithms, averaged over $100$ independent trajectories, in \Cref{figure:log-biases} and \Cref{figure:log-covariance} in logarithmic scale.

The slope for LMC and SGLDFP is $-1$ which confirms the convergence of $\norm{\bar{\theta}_n - \theta^{\star}}$ to $0$ at a rate $N^{-1}$. On the other hand, we can observe that $\norm{\bar{\theta}_n - \theta^{\star}}$ converges to a constant for SGD and SGLD.

\begin{figure}
\centering
\includegraphics[scale=0.4]{loglog-biases-3.pdf}
\caption{\label{figure:log-biases} Distance to $\theta^{\star}$, $\norm{\bar{\theta}_n - \theta^{\star}}$ for LMC, SGLDFP, SGLD and SGD, function of $N$, in logarithmic scale.}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.4]{loglog-covariance-3.pdf}
\caption{\label{figure:log-covariance} Trace of the covariance matrices for LMC, SGLDFP, SGLD and SGD, function of $N$, in logarithmic scale.}
\end{figure}

\paragraph*{Covertype dataset} We then illustrate our results on the covertype dataset\footnote{\url{https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/binary/covtype.libsvm.binary.scale.bz2}} with a Bayesian logistic regression model. The prior is a standard multivariate Gaussian distribution. Given the size of the dataset and the dimension of the problem, LMC requires high computational resources and is not included in the simulations. We truncate the training dataset at $N\in\defEns{10^3, 10^4, 10^5}$. For all algorithms, the step size $\gamma$ is set equal to $1/N$ and the trajectories are started at $\hat{\theta}$, an estimator of $\theta^{\star}$, computed using SGD combined with the BFGS algorithm.

We empirically check that the variance of the stochastic gradients scale as $N^2$ for SGD and SGLD, and as $N$ for SGLDFP. We compute the empirical variance estimator of the gradients, take the mean over the dimension and display the result in a logarithmic plot in \Cref{figure:log-var-grad-cover}. The slopes are $2$ for SGD and SGLD, and $1$ for SGLDFP.

On the test dataset, we also evaluate the negative loglikelihood of the three algorithms for different values of $N\in\defEns{10^3,10^4,10^5}$, as a function of the number of iterations. The plots are shown in \Cref{figure:negll-cover}. We note that for large $N$, SGLD and SGD give very similar results that are below the performance of SGLDFP.

\begin{figure}
\centering
\includegraphics[scale=0.4]{loglog-var-gradients-cover.pdf}
\caption{\label{figure:log-var-grad-cover} Variance of the stochastic gradients of SGLD, SGLDFP and SGD function of $N$, in logarithmic scale.}
\end{figure}

\begin{figure}
\centering
\includegraphics[scale=0.4]{negll-cover.pdf}
\caption{\label{figure:negll-cover} Negative loglikelihood on the test dataset for SGLD, SGLDFP and SGD function of the number of iterations for different values of $N\in\defEns{10^3,10^4,10^5}$.}
\end{figure}

\section{Introduction}
\label{sec:introduction}

Most MCMC algorithms have not been designed to process huge sample sizes, a typical setting in machine learning.  As a result, many classical MCMC methods fail in this context, because the mixing time becomes prohibitively long and the cost per iteration increases proportionally to the number of training samples $N$. The computational cost in standard Metropolis-Hastings algorithm comes from 1) the computation of the proposals, 2) the acceptance/rejection step. Several approaches to solve these issues  have been recently proposed in machine learning and computational statistics.

Among them, the stochastic gradient langevin dynamics (SGLD) algorithm, introduced in \cite{Welling:Teh:2011}, is a popular choice. This method is based on the Langevin Monte Carlo (LMC) algorithm proposed in \cite{grenander:1983,grenander:miller:1994}. Standard versions of LMC require to compute the gradient of the log-posterior at the current fit of the parameter, but avoid the accept/reject step.  The LMC algorithm is a discretization of a continuous-time process, the overdamped Langevin diffusion, which leaves invariant the target distribution $\pi$.  To further reduce the computational cost, SGLD uses unbiased estimators of the gradient of the log-posterior based on subsampling. This method has triggered a huge number of works among others \cite{ahn:welling:balan:2012,Korattikara:Chen:Welling:2014,Ahn:Shahababa:Welling:2014,chen:ding:carin:2015,
chen2014stochastic,Ding:2014,Ma:Chen:Fox:2015,Dubey:2016,bardenet:doucet:holmes:2017} and have been successfully applied to a range of state of the art machine learning problems \cite{Teh:Patterson:2013,li2016scalable}.

The properties of SGLD with decreasing step sizes have been studied in \cite{teh2016consistency}. The two key findings in this work are that 1) the SGLD algorithm converges weakly to the target distribution $\pi$, 2) the optimal rate of convergence to equilibrium scales as $n^{-1/3}$ where $n$ is the number of iterations, see \cite[Section 5]{teh2016consistency}. However, in most of the applications, constant rather than decreasing step sizes are used, see \cite{ahn:welling:balan:2012,chen2014stochastic,JMLR:v18:16-478,
Li:2016:PSG:3016100.3016149,pmlr-v32-satoa14,JMLR:v17:15-494}.
A natural question for the practical design of SGLD is the choice of the minibatch size. This size controls on the one hand the computational complexity of the algorithm per iteration and on the other hand the variance of the gradient estimator.
Non-asymptotic bounds in Wasserstein distance between the marginal distribution of the SGLD iterates and the target distribution $\pi$ have been established in \cite{dalalyan:colt:2017,dalalyan:karagulyan:2017}.  These results highlight the cost of using stochastic gradients and show that, for a given precision $\epsilon$ in Wasserstein distance, the computational cost of the plain SGLD algorithm does not improve over the LMC algorithm; \citet{vollmer:nagapetyan:2017} reports also similar results on the mean square error.

It has been suggested to use control variates to reduce the high variance of the stochastic gradients. For strongly log-concave models, \citet{vollmer:nagapetyan:2017,baker:fearnhead:fox:nemeth:2017} use the mode of the posterior distribution as a reference point and introduce the SGLDFP (Stochastic Gradient Langevin Dynamics Fixed Point) algorithm. \citet{vollmer:nagapetyan:2017,baker:fearnhead:fox:nemeth:2017} provide upper bounds on the mean square error and the Wasserstein distance between the marginal distribution of the iterates of SGLDFP and the posterior distribution. In addition, \citet{vollmer:nagapetyan:2017,baker:fearnhead:fox:nemeth:2017} show that the overall cost remains sublinear in the number of individual data points, up to a preprocessing step. Other control variates methodologies are provided for non-concave models in the form of SAGA-Langevin Dynamics and SVRG-Langevin Dynamics \cite{Dubey:2016,Chen:Wang:Zhang:2017}, albeit a detailed analysis in Wasserstein distance of these algorithms is only available for strongly log-concave models \cite{flammarion:jordan:2018}.

In this paper, we provide further insights on the links between SGLD, SGLDFP, LMC and SGD (Stochastic Gradient Descent). In our analysis, the algorithms are used with a constant step size and the parameters are set to the standard values used in practice \cite{ahn:welling:balan:2012,chen2014stochastic,JMLR:v18:16-478,
Li:2016:PSG:3016100.3016149,pmlr-v32-satoa14,JMLR:v17:15-494}. The LMC, SGLD and SGLDFP algorithms define homogeneous Markov chains, each of which admits a unique stationary distribution used as a hopefully close proxy of $\pi$. The main contribution of this paper is to show that, while the invariant distributions of LMC and SGLDFP become closer to $\pi$ as the number of data points increases, on the opposite, the invariant measure of SGLD never comes close to the target distribution $\pi$ and is in fact very similar to the invariant measure of SGD.

In \Cref{subsec:wasserstein-distance}, we give an upper bound in Wasserstein distance of order $2$ between the marginal distribution of the iterates of LMC and the Langevin diffusion, SGLDFP and LMC, and SGLD and SGD. We provide a lower bound on the Wasserstein distance between the marginal distribution of the iterates of SGLDFP and SGLD. In \Cref{subsec:variance-piula-pifp-puslgd-pisgd}, we give a comparison of the means and covariance matrices of the invariant distributions of LMC, SGLDFP and SGLD with those of the target distribution $\pi$. Our claims are supported by numerical experiments in \Cref{sec:numerical-experiments}.

\subsubsection*{Acknowledgments}

\bibliographystyle{abbrvnat}

\section{Preliminaries}
\label{sec:preliminaries}

Denote by  $\mathbf{z} = \{z_i\}_{i=1}^{N}$ the observations.  We are interested in situations where the target distribution $\pi$ arises as the posterior in a Bayesian inference problem with prior density $\pi_0(\theta)$ and a large number $N \gg 1$ of i.i.d.\ observations $z_i$ with likelihoods $p(z_i|\theta)$. In this case, $\pi(\theta)= \pi_0(\theta) \prod_{i=1}^N p(z_i|\theta)$. We denote $U_i(\theta) = -\log(p(z_i | \theta))$ for $i\in\{1,\ldots,N\}$, $U_0(\theta) = -\log(\pi_0(\theta))$, $U = \sum_{i=0}^{N} U_i$.

Under mild conditions, $\pi$ is the unique invariant probability measure of the Langevin Stochastic Differential Equation (SDE):
\begin{equation}\label{eq:langevin-SDE}
  \mathrm{d} \theta_t = - \nabla U(\theta_t) \mathrm{d} t + \sqrt{2} \mathrm{d} B_t \;,
\end{equation}
where $(B_t)_{t\geq 0}$ is a $d$-dimensional Brownian motion.
Based on this observation, Langevin Monte Carlo (LMC) is an MCMC algorithm that enables to sample (approximately) from $\pi$ using an Euler discretization of the Langevin SDE:
\begin{equation}\label{eq:LMC}
  \theta_{k+1} = \theta_k - \gamma \nabla U(\theta_k) + \sqrt{2\gamma} Z_{k+1} \;,
\end{equation}
where $\gamma>0$ is a constant step size and $(Z_k)_{k\geq 1}$ is a sequence of i.i.d.~standard $d$-dimensional Gaussian vectors. Discovered and popularised in the seminal works \cite{grenander:1983,grenander:miller:1994,roberts:tweedie-Langevin:1996}, LMC has recently received renewed attention \cite{dalalyan:2014,durmus:moulines:2015,durmus:moulines:highdimULA,dalalyan:karagulyan:2017}. However, the cost of one iteration is $N d$ which is prohibitively large for massive datasets. In order to scale up to the big data setting, \citet{Welling:Teh:2011} suggested to replace $\nabla U$ with an unbiased estimate $\nabla U_0 + (N/p)\sum_{i\inS} \nabla U_i$ where $S$ is a minibatch of $\{1,\ldots,N\}$ with replacement of size $p$. A single update of SGLD is then given for $k\in\mathbb{N}$ by
\begin{equation}\label{eq:SGLD}
  \theta_{k+1} = \theta_k - \gamma \parenthese{ \nabla U_0(\theta_k) + \frac{N}{p} \sum_{i\inS_{k+1}} \nabla U_i(\theta_k)} + \sqrt{2\gamma} Z_{k+1} \;.
\end{equation}
The idea of using only a fraction of data points to compute an unbiased estimate of the gradient at each iteration comes from Stochastic Gradient Descent (SGD) which is a popular algorithm to minimize the potential $U$. SGD is very similar to SGLD because it is characterised by the same recursion as SGLD but without Gaussian noise:
\begin{equation}\label{eq:SGD}
  \theta_{k+1} = \theta_k - \gamma \parenthese{\nabla U_0(\theta_k) + \frac{N}{p}\sum_{i\inS_{k+1}} \nabla U_i(\theta_k)} \;.
\end{equation}
Assuming for simplicity that $U$ has a minimizer $\theta^{\star}$, we can define a control variates version of SGLD, SGLDFP, see \cite{Dubey:2016,Chen:Wang:Zhang:2017}, given for $k\in\mathbb{N}$ by
\begin{equation}\label{eq:SGLDFP}
  \theta_{k+1} = \theta_k - \gamma \parenthese{ \nabla U_0(\theta_k) - \nabla U_0(\theta^{\star}) +  \frac{N}{p} \sum_{i\inS_{k+1}} \defEns{\nabla U_i(\theta_k) - \nabla U_i(\theta^{\star})}} + \sqrt{2\gamma} Z_{k+1} \;.
\end{equation}

It is worth mentioning that the objectives of the different algorithms presented so far are distinct. On the one hand, LMC, SGLD and SGDLFP are MCMC methods used to obtain approximate samples from the posterior distribution $\pi$. On the other hand, SGD is a stochastic optimization algorithm used to find an estimate of the mode $\theta^{\star}$ of the posterior distribution. In this paper, we focus on the fixed step-size SGLD algorithm and assess its ability to reliably sample from $\pi$. For that purpose and to quantify precisely the relation between LMC, SGLD, SGDFP and SGD, we make for simplicity the following assumptions on $U$.

\begin{assumption}\label{assum:gradient-Lipschitz}
  For all $i\in\defEns{0,\ldots,N}$, $U_i$ is four times continuously differentiable and for all $j \in \{2,3,4\}$, $\sup_{\theta\in \mathbb{R}^d}\norm{\operatorname{D}^j U_i(\theta)} \leq \tilde{L}$. In particular for all $i \in \{0,\ldots,N\}$, $U_i$ is $\tilde{L}$-gradient Lipschitz, \textit{i.e.}~for all $\theta_1,\theta_2\in\mathbb{R}^d$, $\norm{\nabla U_i(\theta_1) - \nabla U_i(\theta_2)} \leq \tilde{L} \norm{\theta_1 - \theta_2}$.
\end{assumption}

\begin{assumption}\label{assum:strongly-convex}
  $U$ is $m$-strongly convex, \textit{i.e.}~for all $\theta_1,\theta_2\in\mathbb{R}^d$, $\ps{\nabla U(\theta_1) - \nabla U (\theta_2)}{\theta_1 - \theta_2} \geq m\norm[2]{\theta_1 - \theta_2}$.
\end{assumption}

\begin{assumption}\label{assum:Ui-convex}
  For all $i\in\defEns{0,\ldots,N}$, $U_i$ is convex.
\end{assumption}

Note that under \Cref{assum:gradient-Lipschitz}, $U$ is four times continuously differentiable and for $j \in \{2,3,4\}$, $\sup_{\theta\in \mathbb{R}^d}\norm{\operatorname{D}^j U(\theta)} \leq L$, with $L = (N+1)\tilde{L}$ and where $\norm{\operatorname{D}^j U(\theta)} = \sup_{\norm{u_1} \leq 1,\ldots,\norm{u_j} \leq 1} \operatorname{D}^j U(\theta)[u_1,\ldots,u_j]$. In particular, $U$ is  $L$-gradient Lipschitz.  Furthermore, under \Cref{assum:strongly-convex}, $U$ has a unique minimizer $\theta^{\star}$. In this paper, we focus on the asymptotic $N \to +\infty$,. We assume that $\liminf_{N\to+\infty} N^{-1}m >0$, which is a common assumption for the analysis of SGLD and SGLDFP \cite{baker:fearnhead:fox:nemeth:2017,flammarion:jordan:2018}. In practice \cite{ahn:welling:balan:2012,chen2014stochastic,JMLR:v18:16-478,
Li:2016:PSG:3016100.3016149,pmlr-v32-satoa14,JMLR:v17:15-494}, $\gamma$ is of order $1/N$ and we adopt this convention in this article.

For a practical implementation of SGLDFP, an estimator $\hat{\theta}$ of $\theta^{\star}$ is necessary. The theoretical analysis and the bounds remain unchanged if, instead of considering SGLDFP centered w.r.t.~$\theta^{\star}$, we study SGLDFP centered w.r.t.~$\hat{\theta}$ satisfying $\mathbb{E}[\Vert \hat{\theta} - \theta^{\star} \Vert^2] = O(1/N)$. Such an estimator $\hat{\theta}$ can be computed using for example SGD with decreasing step sizes, see \cite[eq.(2.8)]{nemirovski:juditsky:lan:shapiro:2009} and \cite[Section 3.4]{baker:fearnhead:fox:nemeth:2017}, for a computational cost linear in $N$.

\section{Proofs of \Cref{subsec:wasserstein-distance}}
\label{subsec:proof-section-wasserstein-distance}

\subsection{Proof of \Cref{lemma:existence-pi-gamma}}
\label{subsec:proof-lemma-existence-pi-gamma}

The convergence in Wasserstein distance is classically done via a standard synchronous coupling \cite[Proposition 2]{dieuleveut:durmus:bach:2017}. We prove the statement for SGLD; the adaptation for LMC, SGLDFP and SGD is immediate.
Let $\gamma\in\ooint{0,2/L}$ and $\lambda_1,\lambda_2\in\mathcal{P}_2(\mathbb{R}^d)$. By
\cite[Theorem 4.1]{VillaniTransport}, there exists a couple of random variables $(\theta^{(1)}_0, \theta^{(2)}_0)$ such that $\wasser[2]^2(\lambda_1, \lambda_2)= \expe{\norm[2]{\theta^{(1)}_0 - \theta^{(2)}_0}}$.
Let $(\theta^{(1)}_k,\theta^{(2)}_k)_{k\in\mathbb{N}}$ be the SGLD iterates starting from $\theta^{(1)}_0$ and $\theta^{(2)}_0$ respectively and driven by the same noise, \textit{i.e.}~for all $k\in\mathbb{N}$,
\begin{equation*}
  \begin{cases}
    \theta^{(1)}_{k+1} & = \theta^{(1)}_k - \gamma \defEns{\nabla U_0(\theta^{(1)}_k) + (N/p)\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(1)}_k)} + \sqrt{2\gamma} Z_{k+1} \;,\\
    \theta^{(2)}_{k+1} & = \theta^{(2)}_k - \gamma\defEns{\nabla U_0(\theta^{(2)}_k) +(N/p)\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(2)}_k)} + \sqrt{2\gamma} Z_{k+1} \;,
  \end{cases}
\end{equation*}
where $(Z_k)_{k\geq 1}$ is an i.i.d.~sequence of standard Gaussian variables and $(S_k)_{k\geq 1}$ an i.i.d.~sequence of subsamples of $\defEns{1,\ldots,N}$ of size $p$. Denote by $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the filtration associated to $(\theta^{(1)}_k,\theta^{(2)}_k)_{k\in\mathbb{N}}$. We have for $k\in\mathbb{N}$,
\begin{align*}
  &\norm[2]{\theta^{(1)}_{k+1} - \theta^{(2)}_{k+1}} =  \\
  &\norm[2]{\theta^{(1)}_{k} - \theta^{(2)}_{k}} + \gamma^2 \norm[2]{\nabla U_0(\theta^{(1)}_{k}) + \frac{N}{p}\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(1)}_{k}) - \nabla U_0(\theta^{(2)}_{k}) - \frac{N}{p}\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(2)}_{k})} \\
  &- 2\gamma\ps{\theta^{(1)}_{k} - \theta^{(2)}_{k}}{\nabla U_0(\theta^{(1)}_{k}) + \frac{N}{p}\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(1)}_{k}) - \nabla U_0(\theta^{(2)}_{k}) - \frac{N}{p}\sum_{i\inS_{k+1}} \nabla U_i(\theta^{(2)}_{k})} \;.
\end{align*}
By \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex}, $\theta\mapsto \nabla U_0(\theta) +  (N/p)\sum_{i\inS} \nabla U_i(\theta)$ is $\mathbb{P}$-\ensuremath{\text{a.s.}}~$L$-co-coercive \cite{Zhu:Marcotte:1996}. Taking the conditional expectation w.r.t.~$\mathcal{F}_k$, we obtain
\begin{equation*}
  \expecond{\norm[2]{\theta^{(1)}_{k+1} - \theta^{(2)}_{k+1}}}{\mathcal{F}_k} \leq \norm[2]{\theta^{(1)}_{k} - \theta^{(2)}_{k}} - 2\gamma\defEns{1-(\gammaL)/2} \ps{\theta^{(1)}_{k} - \theta^{(2)}_{k}}{\nabla U(\theta^{(1)}_{k}) - \nabla U(\theta^{(2)}_{k})} \;,
\end{equation*}
and by \Cref{assum:strongly-convex}
\begin{equation*}
  \expecond{\norm[2]{\theta^{(1)}_{k+1} - \theta^{(2)}_{k+1}}}{\mathcal{F}_k} \leq \defEns{1-2m\gamma(1-(\gammaL)/2)} \norm[2]{\theta^{(1)}_{k} - \theta^{(2)}_{k}} \;.
\end{equation*}
Since for all $k\geq 0$, $(\theta^{(1)}_k, \theta^{(2)}_k)$ belongs to $\Pi(\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k)$, we get by a straightforward induction
\begin{equation}\label{eq:contraction-wasserstein}
  \wasser[2]^2(\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) \leq \expe{\norm[2]{\theta^{(1)}_{k} - \theta^{(2)}_{k}}} \leq \defEns{1-2m\gamma(1-(\gammaL)/2)}^k \wasser[2]^2(\lambda_1,\lambda_2) \;.
\end{equation}
By \Cref{assum:gradient-Lipschitz}, $\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}} \in\mathcal{P}_2(\mathbb{R}^d)$ and taking $\lambda_2 = \lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}$, we get $\sum_{k=0}^{+\infty} \wasser[2]^2(\lambda_1\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k, \lambda_1\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^{k+1}) < +\infty$.By \cite[Theorem 6.16]{VillaniTransport}, $\mathcal{P}_2(\mathbb{R}^d)$ endowed with $\wasser[2]$ is a Polish space. $(\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k)_{k\geq 0}$ is a Cauchy sequence and converges to a limit $
\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_1}\in\mathcal{P}_2(\mathbb{R}^d)$. The limit $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_1}$ does not depend on $\lambda_1$ because, given $\lambda_2\in\mathcal{P}_2(\mathbb{R}^d)$, by the triangle inequality
\begin{equation*}
  \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_1}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_2}) \leq \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_1}, \lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) + \wasser[2](\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) + \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_2}, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) \;.
\end{equation*}
Taking the limit $k\to+\infty$, we get $\wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_1}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}^{\lambda_2})=0$. The limit is thus the same for all initial distributions and is denoted $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$. $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ is invariant for $\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}$ since we have for all $k\in\mathbb{N}^*$,
\begin{equation*}
  \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}) \leq \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) + \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k) \;.
\end{equation*}
Taking the limit $k\to+\infty$, we obtain $\wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}) = 0$. Using \eqref{eq:contraction-wasserstein}, $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ is the unique invariant probability measure for $\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}$.

\subsection{Proof of \Cref{thm:wasserstein-ula-sgld-sgd-fp}}
\label{subsec:proof-thm-wasserstein-ula-sgld-sgd-fp}

Proof of \ref{item:W-bound-ULA-SGLDFP}.
Let $\gamma\in\ocint{0,1/L}$ and $\lambda_1,\lambda_2\in\mathcal{P}_2(\mathbb{R}^d)$. By
\cite[Theorem~4.1]{VillaniTransport}, there exists a couple of random variables $(\theta_0, \vartheta_0)$ such that $\wasser[2]^2(\lambda_1, \lambda_2)= \expe{\norm[2]{\theta_0 - \vartheta_0}}$.
Let $(\theta_k,\vartheta_k)_{k\in\mathbb{N}}$ be the LMC and SGLDFP iterates starting from $\theta_0$ and $\vartheta_0$ respectively and driven by the same noise, \textit{i.e.}~for all $k\in\mathbb{N}$,
\begin{equation*}
  \begin{cases}
    \theta_{k+1} & = \theta_k - \gamma \nabla U(\theta_k) + \sqrt{2\gamma} Z_{k+1} \;,\\
    \vartheta_{k+1} & = \vartheta_k - \gamma\parenthese{\nabla U_0(\vartheta_k) - \nabla U_0(\theta^{\star}) +(N/p)\sum_{i\inS_{k+1}} \defEns{\nabla U_i(\vartheta_k) - \nabla U_i(\theta^{\star})}} + \sqrt{2\gamma} Z_{k+1} \;,
  \end{cases}
\end{equation*}
where $(Z_k)_{k\geq 1}$ is an i.i.d.~sequence of standard Gaussian variables and $(S_k)_{k\geq 1}$ an i.i.d.~sequence of subsamples with replacement of $\defEns{1,\ldots,N}$ of size $p$. Denote by $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the filtration associated to $(\theta_k,\vartheta_k)_{k\in\mathbb{N}}$. We have for $k\in\mathbb{N}$,
\begin{equation}\label{eq:wasserstein-ULA-FP-1}
  \expecond{\norm[2]{\theta_{k+1} - \vartheta_{k+1}}}{\mathcal{F}_k} = \norm[2]{\theta_k - \vartheta_k} - 2\gamma\ps{\theta_k - \vartheta_k}{\nabla U(\theta_k) - \nabla U(\vartheta_k)} + \gamma^2 A
\end{equation}
where
\begin{align*}
  A &= \expecond{\norm[2]{\nabla U(\theta_k) - \parenthese{\nabla U_0(\vartheta_k) - \nabla U_0(\theta^{\star}) +(N/p)\sum_{i\inS_{k+1}} \defEns{\nabla U_i(\vartheta_k) - \nabla U_i(\theta^{\star})}}}}{\mathcal{F}_k} \\
  &= A_1 + A_2 \;, \\
  A_1 & = \norm[2]{\nabla U(\theta_k) - \nabla U(\vartheta_k)} \;, \\
  A_2 & = \expecond{\norm[2]{\nabla U(\vartheta_k) - \parenthese{\nabla U_0(\vartheta_k) - \nabla U_0(\theta^{\star}) +(N/p)\sum_{i\inS_{k+1}} \defEns{\nabla U_i(\vartheta_k) - \nabla U_i(\theta^{\star})}}}}{\mathcal{F}_k} \;.
\end{align*}
Denote by $W$ the random variable equal to $\nabla U_i(\vartheta_k) - \nabla U_i(\theta^{\star}) - (1/N)\sum_{j=1}^{N} \{ \nabla U_j(\vartheta_k) - \nabla U_j(\theta^{\star}) \}$ for $i\in\defEns{1,\ldots,N}$ with probability $1/N$. By \Cref{assum:gradient-Lipschitz} and using the fact that the subsamples $(S_k)_{k\geq 1}$ are drawn with replacement, we obtain
\begin{equation*}
  A_2 = (N^2/p) \expe{\norm[2]{W} | \mathcal{F}_k} \leq (L^2 / p) \norm[2]{\vartheta_k - \theta^{\star}} \;.
\end{equation*}
Combining it with \eqref{eq:wasserstein-ULA-FP-1}, and using the $L$-co-coercivity of $\nabla U$ under \Cref{assum:gradient-Lipschitz} and \Cref{assum:strongly-convex}, we get
\begin{equation*}
  \expecond{\norm[2]{\theta_{k+1} - \vartheta_{k+1}}}{\mathcal{F}_k} \leq (1-m\gamma) \norm[2]{\theta_k - \vartheta_k} + \{(L^2\gamma^2)/p\} \norm[2]{\vartheta_k - \theta^{\star}} \;.
\end{equation*}
Iterating and using \Cref{lemma:bound-tp-tpstar}-\ref{item:bound-tp-tpstar-sgldfp-ula}, we have for $n\in\mathbb{N}$
\begin{align*}
  & \wasser[2]^2(\lambda_1 \Rker_{\scriptscriptstyle{\operatorname{LMC}}}^n, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{FP}}}^n) \leq \expe{\norm[2]{\theta_n - \vartheta_n}} \\
  &\leq (1-m\gamma)^n \wasser[2]^2(\lambda_1,\lambda_2) + \frac{L^2\gamma^2}{p} \sum_{k=0}^{n-1} (1-m\gamma)^{n-1-k} \expe{\norm[2]{\vartheta_k - \theta^{\star}}} \\
  & \leq (1-m\gamma)^n \wasser[2]^2(\lambda_1,\lambda_2) + \frac{L^2\gamma^2}{p} n (1-m\gamma)^{n-1} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda_2(\mathrm{d} \vartheta) + \frac{2L^2\gamma d}{p m^2} \;.
\end{align*}

Proof of \ref{item:W-bound-ULA-Langevin}.
Denote by $\kappa=(2\mUL)/(m+L)$.
By \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \cite[Theorem 5]{durmus:moulines:highdimULA}, we have for all $n\in\mathbb{N}$,
\begin{multline*}
  \wasser[2]^2(\lambda_1 P_{n\gamma}, \lambda_2 \Rker_{\scriptscriptstyle{\operatorname{LMC}}}^n) \leq 2\parenthese{1 - \kappa\gamma/2}^n \wasser[2]^2(\lambda_1,\lambda_2) + \frac{2L^2\gamma}{\kappa} (\kappa^{-1}+\gamma)\parenthese{2d+\frac{dL^2\gamma^2}{6}} \\
  + L^4\gamma^3 (\kappa^{-1}+\gamma) \sum_{k=1}^{n} \delta_k \defEns{1-\kappa\gamma/2}^{n-k}
\end{multline*}
where for all $k\in\defEns{1,\ldots,n}$,
\begin{equation*}
  \delta_k \leq \mathrm{e}^{-2m(k-1)\gamma} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda_1 (\mathrm{d} \vartheta) + d/m \;.
\end{equation*}
We get the result by straightforward simplifications and using $\gamma\leq 1/L$.

Proof of \ref{item:W-bound-SGLD-SGD}.
Let $\gamma\in\ocint{0,1/L}$ and $\lambda_1,\lambda_2\in\mathcal{P}_2(\mathbb{R}^d)$. By
\cite[Thereom 4.1]{VillaniTransport}, there exists a couple of random variables $(\theta_0, \vartheta_0)$ such that $\wasser[2]^2(\lambda_1, \lambda_2)= \expe{\norm[2]{\theta_0 - \vartheta_0}}$.
Let $(\theta_k,\vartheta_k)_{k\in\mathbb{N}}$ be the SGLD and SGD iterates starting from $\theta_0$ and $\vartheta_0$ respectively and driven by the same noise, \textit{i.e.}~for all $k\in\mathbb{N}$,
\begin{equation*}
  \begin{cases}
    \theta_{k+1} & = \theta_k - \gamma\parenthese{\nabla U_0(\theta_k)+(N/p)\sum_{i\inS_{k+1}} \nabla U_i(\theta_k)} + \sqrt{2\gamma} Z_{k+1} \;,\\
    \vartheta_{k+1} & = \vartheta_k - \gamma\parenthese{\nabla U_0(\vartheta_k) +(N/p)\sum_{i\inS_{k+1}} \nabla U_i(\vartheta_k)} \;,
  \end{cases}
\end{equation*}
where $(Z_k)_{k\geq 1}$ is an i.i.d.~sequence of standard Gaussian variables and $(S_k)_{k\geq 1}$ an i.i.d.~sequence of subsamples with replacement of $\defEns{1,\ldots,N}$ of size $p$. Denote by $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the filtration associated to $(\theta_k,\vartheta_k)_{k\in\mathbb{N}}$. We have for $k\in\mathbb{N}$,
\begin{multline*}
  \expecond{\norm[2]{\theta_{k+1} - \vartheta_{k+1}}}{\mathcal{F}_k} = \norm[2]{\theta_k - \vartheta_k} - 2\gamma\ps{\theta_k - \vartheta_k}{\nabla U(\theta_k) - \nabla U(\vartheta_k)} + 2\gamma d \\
  + \gamma^2 \expecond{\norm[2]{\nabla U_0(\theta_k)+(N/p)\sum_{i\inS_{k+1}} \nabla U_i(\theta_k) - \nabla U_0(\vartheta_k)-(N/p)\sum_{i\inS_{k+1}} \nabla U_i(\vartheta_k)}}{\mathcal{F}_k} \;.
\end{multline*}
By \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex}, $\theta\mapsto \nabla U_0(\theta) +  (N/p)\sum_{i\inS} \nabla U_i(\theta)$ is $\mathbb{P}$-\ensuremath{\text{a.s.}}~$L$-co-coercive and we obtain
\begin{equation*}
  \expecond{\norm[2]{\theta_{k+1} - \vartheta_{k+1}}}{\mathcal{F}_k} \leq \{1-2m\gamma(1-\gammaL/2)\} \norm[2]{\theta_k - \vartheta_k} + 2\gamma d \;,
\end{equation*}
which concludes the proof by a straightforward induction.

\subsection{Proof of \Cref{thm:W-lower-bound-SGLD-SGLDFP}}
\label{subsec:proof-thm-W-lower-bound-SGLD-SGLDFP}

Proof of \ref{item:W-lower-bound-SGLD-FP-1}.
Let $\gamma\in\ocint{0,\Sigma^{-1}\{1+N/(p\sum_{i=1}^{N}x_i^2)\}^{-1}}$, $(\theta_k)_{k\in\mathbb{N}}$ be the iterates of SGLD \eqref{eq:SGLD} started at $\theta^{\star}$ and $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the associated filtration. For all $k\in\mathbb{N}$, $\expe{\theta_k} = \theta^{\star}$. The variance of $\theta_k$ satisfies the following recursion for $k\in\mathbb{N}$
\begin{align*}
  & \expecond{(\theta_{k+1} - \theta^{\star})^2}{\mathcal{F}_k} \\
  &= \expecond{\defEns{\theta_k - \theta^{\star} - \gamma\parenthese{\Sigma(\theta_k - \theta^{\star}) + \rho(S_{k+1})(\theta_k - \theta^{\star}) + \xi(S_{k+1})} + \sqrt{2\gamma} Z_{k+1}}^2}{\mathcal{F}_k} \\
  & = \mu (\theta_k - \theta^{\star})^2 + 2\gamma + \gamma^2 A \;,
\end{align*}
where
\begin{equation*}
  \mu = \expe{\defEns{1-\gamma\parenthese{\frac{1}{\sigma^2_{\tp}}+\frac{N}{\sigma^2_y p}\sum_{i\inS} x_i^2}}^2} \;,\quad
  A = \expe{\defEns{\frac{\theta^{\star}}{\sigma^2_{\tp}} + \frac{N}{\sigdyp}\sum_{i\inS} \parenthese{x_i\theta^{\star} - y_i}x_i}^2} \;.
\end{equation*}
We have for $\mu$,
\begin{align*}
  \mu & = 1-2\gamma\Sigma + \gamma^2 \expe{\defEns{\frac{N}{\sigma^2_y p}\sum_{i\inS} x_i^2 - \frac{1}{\sigma^2_y}\sum_{i=1}^{N}x_i^2}^2} +\gamma^2\Sigma^2\\
  & = 1 - 2\gamma\Sigma + \gamma^2\defEns{\Sigma^2 + \frac{N}{\sigma^4_y p}\sum_{i=1}^{N}\parenthese{x_i^2 -\frac{1}{N}\sum_{j=1}^{N}x_j^2}} \leq 1 - \gamma\Sigma \;,
\end{align*}
and for $A$,
\begin{equation*}
  A = \frac{N}{p}\sum_{i=1}^{N}\defEns{\frac{\parenthese{x_i\theta^{\star} - y_i}x_i}{\sigma^2_y} + \frac{\theta^{\star}}{N\sigma^2_{\tp}}}^2 \;.
\end{equation*}
By a straightforward induction, we obtain that the variance of the $n^{\text{th}}$ iterate of SGLD started at $\theta^{\star}$ is for $n\in\mathbb{N}^*$
\begin{equation*}
  \int_{\mathbb{R}} (\theta-\theta^{\star})^2 \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^n(\theta^{\star}, \mathrm{d} \theta) = \frac{1-\mu^n}{1-\mu} 2\gamma + \frac{1-\mu^n}{1-\mu} \frac{N\gamma^2}{p}\sum_{i=1}^{N}\defEns{\frac{\parenthese{x_i\theta^{\star} - y_i}x_i}{\sigma^2_y} + \frac{\theta^{\star}}{N\sigma^2_{\tp}}}^2 \;.
\end{equation*}
For SGLDFP, the additive part of the noise in the stochastic gradient disappears and we obtain similarly for $n\in\mathbb{N}^*$
\begin{equation*}
  \int_{\mathbb{R}} (\theta-\theta^{\star})^2 \Rker_{\scriptscriptstyle{\operatorname{FP}}}^n(\theta^{\star}, \mathrm{d} \theta) = \frac{1-\mu^n}{1-\mu} 2\gamma \; .
\end{equation*}
To conclude, we use that for two probability measures with given mean and covariance matrices, the Wasserstein distance between the two Gaussians with these respective parameters is a lower bound for the Wasserstein distance between the two measures \cite[Theorem 2.1]{gelbrich:1990}.

The proof of \ref{item:W-lower-bound-SGLD-FP-2} is straightforward.

\section{Proofs of \Cref{subsec:variance-piula-pifp-puslgd-pisgd}}
\label{subsec:proof-section-bias-variance}

\subsection{Proof of \Cref{lemma:mean-pi-tpstar}}
\label{subsubsec:proof-lemma-mean-pi-tpstar}

Let $\theta$ be distributed according to $\pi$. By \Cref{assum:strongly-convex}, for all $\vartheta\in\mathbb{R}^d$, $U(\vartheta) \geq U(\theta^{\star}) + (m/2)\norm[2]{\vartheta - \theta^{\star}}$ and $\expe{\nabla U(\theta)}=0$. By a Taylor expansion of $\nabla U$ around $\theta^{\star}$, we obtain
\begin{equation*}
  0 = \expe{\nabla U(\theta)} = \nabla^2 U(\theta^{\star}) \parenthese{\expe{\theta} - \theta^{\star}} + (1/2) \operatorname{D}^3 U(\theta^{\star})[\expe{(\theta - \theta^{\star})^{\otimes 2}}] + \expe{\reste{1}(\theta)} \;,
\end{equation*}
where by \Cref{assum:gradient-Lipschitz}, $\reste{1}:\mathbb{R}^d \to \mathbb{R}^d$ satisfies
\begin{equation}\label{eq:proof-mean-pipistar-reste-1}
  \sup_{\vartheta\in\mathbb{R}^d} \defEns{\norm{\reste{1}(\vartheta)} / \norm[3]{\vartheta - \theta^{\star}}} \leq L / 6 \;.
\end{equation}
Rearranging the terms, we get
\begin{equation*}
  \expe{\theta} - \theta^{\star} = -(1/2) \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star})[\expe{(\theta - \theta^{\star})^{\otimes 2}}] - \nabla^2 U(\theta^{\star})^{-1} \expe{\reste{1}(\theta)} \;.
\end{equation*}
To estimate the covariance matrix of $\pi$ around $\theta^{\star}$, we start again from the Taylor expansion of $\nabla U$ around $\theta^{\star}$ and we obtain
\begin{equation}
  \label{eq:mean-pi-tpstar_1}
  \expe{\nabla U(\theta)^{\otimes 2}} = \expe{\parenthese{\nabla^2 U(\theta^{\star})(\theta - \theta^{\star}) + \reste{2}(\theta)}^{\otimes 2}}
  = \nabla^2 U(\theta^{\star})^{\otimes 2} \expe{(\theta - \theta^{\star})^{\otimes 2}} + \expe{\reste{3}(\theta)}
\end{equation}
where by \Cref{assum:gradient-Lipschitz}, $\reste{2}:\mathbb{R}^d \to \mathbb{R}^d$ satisfies
\begin{equation}\label{eq:proof-mean-pipistar-reste-2}
  \sup_{\vartheta\in\mathbb{R}^d} \defEns{\norm{\reste{2}(\vartheta)} / \norm[2]{\vartheta - \theta^{\star}}} \leq L / 2 \;,
\end{equation}
and $\reste{3}:\mathbb{R}^d \to \mathbb{R}^{d\times d}$ is defined for all $\vartheta\in\mathbb{R}^d$ by
\begin{equation}\label{eq:proof-mean-pipistar-reste-3}
  \reste{3}(\vartheta) = \nabla^2 U(\theta^{\star})(\vartheta - \theta^{\star}) \otimes \reste{2}(\vartheta) + \reste{2}(\vartheta) \otimes \nabla^2 U(\theta^{\star})(\vartheta - \theta^{\star}) + \reste{2}(\vartheta)^{\otimes 2} \;.
\end{equation}
$\expe{\nabla U(\theta)^{\otimes 2}}$ is the Fisher information matrix and by a Taylor expansion of $\nabla^2 U$ around $\theta^{\star}$ and an integration by parts,
\begin{equation*}
  \expe{\nabla U(\theta)^{\otimes 2}} = \expe{\nabla^2 U(\theta)} = \nabla^2 U(\theta^{\star}) + \expe{\reste{4}(\theta)}
\end{equation*}
where by \Cref{assum:gradient-Lipschitz}, $\reste{4}:\mathbb{R}^d \to \mathbb{R}^{d\times d}$ satisfies
\begin{equation}\label{eq:proof-mean-pipistar-reste-4}
  \sup_{\vartheta\in\mathbb{R}^d} \defEns{\norm{\reste{4}(\vartheta)} / \norm{\vartheta - \theta^{\star}}} \leq L \;.
\end{equation}
Combining this result, \eqref{eq:proof-mean-pipistar-reste-1}, \eqref{eq:mean-pi-tpstar_1}, \eqref{eq:proof-mean-pipistar-reste-2}, \eqref{eq:proof-mean-pipistar-reste-3}, \eqref{eq:proof-mean-pipistar-reste-4} and $\expeLigne{\norm{\theta-\theta^{\star}}^4} \leq d(d+2)/m^2$ by \cite[Lemma 9]{2017arXiv171005559B} conclude the proof.

\subsection{Proofs of \Cref{thm:bias-variance-tpstar-ula-fp} and \Cref{thm:bias-variance-tpstar-sgld-sgd}}
\label{subsec:proof-section-bias-variance}

First note that under \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex},  there exists $r \in \ccintLigne{0,L/(\sqrt{p} m)}$ such that
\begin{equation}
  \label{eq:assum:r-dom-nabla2Utpstar}
 \opL \preceq r^2 (\nabla^2 U(\theta^{\star}))^{\otimes 2} \;,
\end{equation}
\textit{i.e.}~for all $A\in\mathbb{R}^{d\times d}$,
\[ \operatorname{Tr}(A^{\operatorname{T}} \opL (A)) \leq r^2 \operatorname{Tr}(A^{\operatorname{T}} (\nabla^2 U(\theta^{\star}))^{\otimes 2} A) \;, \]
and where $\opL$ is defined in \eqref{eq:def-operator-L}.
In addition, if $\liminf_{N\to+\infty} N^{-1}m>0$, $r$ can be chosen independently of $N$.

Moreover, for all $\gamma\in\oointLigne{0,2/L}$, $\operatorname{H}$ defined in \eqref{eq:def-operator-H}, is invertible and for all $\gamma\in\oointLigne{0,2/\{(1+r^2)L\}}$, $\operatorname{G}$ defined in \eqref{eq:def-operator-G}, is invertible. Indeed,
\begin{align*}
  \operatorname{H} & = \nabla^2 U(\theta^{\star}) \otimes \parenthese{\operatorname{Id} - \frac{\gamma}{2} \nabla^2 U(\theta^{\star})} + \parenthese{\operatorname{Id} - \frac{\gamma}{2} \nabla^2 U(\theta^{\star})} \otimes \nabla^2 U(\theta^{\star}) \succ 0 \;, \\
  \operatorname{G} &\succeq \nabla^2 U(\theta^{\star}) \otimes \operatorname{Id} + \operatorname{Id} \otimes \nabla^2 U(\theta^{\star}) - \gamma(1+r^2) \nabla^2 U(\theta^{\star}) \otimes \nabla^2 U(\theta^{\star}) \\
  &\succeq \nabla^2 U(\theta^{\star}) \otimes \parenthese{\operatorname{Id} - \frac{\gamma(1+r^2)}{2}\nabla^2 U(\theta^{\star})} + \parenthese{\operatorname{Id} - \frac{\gamma(1+r^2)}{2}\nabla^2 U(\theta^{\star})} \otimes \nabla^2 U(\theta^{\star}) \succ 0 \;.
\end{align*}

For simplicity of notation, in this Section, we use $\epsilon(\theta)$ to denote the difference between the stochastic and the exact gradients at $\theta\in\mathbb{R}^d$. More precisely, $\epsilon$ is the null function for LMC and is defined for $\theta\in\mathbb{R}^d$ by
\begin{align}
  \label{eq:def-gesp-sgld}
  \epsilon(\theta) & = \frac{N}{p}\sum_{i\inS} \nabla U_i(\theta) - \sum_{j=1}^{N} \nabla U_j(\theta) \quad \text{for SGLD and SGD,} \\
  \label{eq:def-gesp-sgldfp}
  \epsilon(\theta) &= \nabla U_0(\theta) - \nabla U_0(\theta^{\star}) + \frac{N}{p}\sum_{i\inS} \defEns{\nabla U_i(\theta) - \nabla U_i(\theta^{\star})} - \nabla U(\theta) \quad \text{for SGLDFP,}
\end{align}
where $S$ is a random subsample of $\{1,\ldots,N\}$ with replacement of size $p\in\mathbb{N}^*$. In this setting, the update equation for LMC, SGLD and SGLDFP is given for $k\in\mathbb{N}$ by
\begin{equation}\label{eq:update-gesp}
  \theta_{k+1} = \theta_k - \parenthese{\nabla U(\theta_k) + \epsilon_{k+1}(\theta_k)} + \sqrt{2\gamma} Z_{k+1} \;,
\end{equation}
where $(Z_k)_{k\geq 1}$ is a sequence of i.i.d.~standard $d$-dimensional Gaussian variables and the sequence of vector fields $(\epsilon_k)_{k\geq 1}$ is associated to a sequence $(S_k)_{k\geq 1}$ of i.i.d.~random subsample of $\defEns{1,\ldots,N}$ with replacement of size $p\in\mathbb{N}^*$.
We also denote by $\bar{\pi}\in\mathcal{P}_2(\mathbb{R}^d)$ the invariant probability measure of LMC, SGLDFP or SGLD.

\subsubsection{Control of the moments of order $2$ and $4$ of LMC, SGLDFP and SGLD}
\label{subsec:control-moments-ula-sgldfp-sgld}

\begin{lemma}\label{lemma:bound-tp-tpstar}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  \begin{enumerate}[label=\roman*)]
    \item \label{item:bound-tp-tpstar-sgldfp-ula}
    For all initial distribution $\lambda\in\mathcal{P}_2(\mathbb{R}^d)$, $\gamma\in\ocint{0,1/L}$ and $k\in\mathbb{N}$,
    \begin{equation*}
      \expe{\norm[2]{\theta_k - \theta^{\star}}} \leq (1-m\gamma)^k \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) + (2d)/m
    \end{equation*}
    where $(\theta_k)_{k\in\mathbb{N}}$ are the iterates of SGLDFP \eqref{eq:SGLDFP} or LMC \eqref{eq:LMC}.
    \item \label{item:bound-tp-tpstar-sgld}
    For all initial distribution $\lambda\in\mathcal{P}_2(\mathbb{R}^d)$, $\gamma\in\ocint{0,1/(2L)}$ and $k\in\mathbb{N}$,
    \begin{multline*}
      \expe{\norm[2]{\theta_k - \theta^{\star}}} \leq (1-m\gamma)^k \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) + \frac{2d}{m} \\
      + \frac{2\gamma N}{m p} \sum_{i=1}^{N} \norm[2]{\nabla U_i(\theta^{\star}) - \frac{1}{N}\sum_{j=1}^{N} \nabla U_j(\theta^{\star})}
    \end{multline*}
    where $(\theta_k)_{k\in\mathbb{N}}$ are the iterates of SGLD \eqref{eq:SGLD}.
  \end{enumerate}
\end{lemma}

\begin{proof}
\ref{item:bound-tp-tpstar-sgldfp-ula}.
We prove the result for SGLDFP, the case of LMC is identical.
Let $\gamma\in\ocint{0,1/L}$, $(\theta_k)_{k\in\mathbb{N}}$ be the iterates of SGLDFP and $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the filtration associated to $(\theta_k)_{k\in\mathbb{N}}$. By \eqref{eq:SGLDFP}, we have for all $k\in\mathbb{N}$,
\begin{multline*}
  \expecond{\norm[2]{\theta_{k+1} - \theta^{\star}}}{\mathcal{F}_k} = \norm[2]{\theta_k - \theta^{\star}} - 2\gamma\ps{\theta_k - \theta^{\star}}{\nabla U(\theta_k) - \nabla U(\theta^{\star})} + 2\gamma d \\
  + \gamma^2 \expecond{ \norm[2]{\nabla U_0(\theta_k) - \nabla U_0(\theta^{\star}) +  \frac{N}{p} \sum_{i\inS_{k+1}} \defEns{\nabla U_i(\theta_k) - \nabla U_i(\theta^{\star})} }}{\mathcal{F}_k}
\end{multline*}
By \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex}, $\theta\mapsto \nabla U_0(\theta) - \nabla U_0(\theta^{\star}) +  (N/p)\sum_{i\inS} \{\nabla U_i(\theta) - \nabla U_i(\theta^{\star})\}$ is $\mathbb{P}$-\ensuremath{\text{a.s.}}~$L$-co-coercive and we obtain
\begin{equation*}
  \expecond{\norm[2]{\theta_{k+1} - \theta^{\star}}}{\mathcal{F}_k} \leq \defEns{1 - 2m\gamma(1-\gammaL/2)}\norm[2]{\theta_k - \theta^{\star}} + 2\gamma d \;.
\end{equation*}
A straightforward induction concludes the proof.

 \ref{item:bound-tp-tpstar-sgld}.
Let $\gamma\in\ocint{0,1/(2L)}$, $(\theta_k)_{k\in\mathbb{N}}$ be the iterates of SGLD and $(\mathcal{F}_k)_{k\in\mathbb{N}}$ the filtration associated to $(\theta_k)_{k\in\mathbb{N}}$. By \eqref{eq:SGLD}, we have for all $k\in\mathbb{N}$,
\begin{align*}
  \expecond{\norm[2]{\theta_{k+1} - \theta^{\star}}}{\mathcal{F}_k} &= \norm[2]{\theta_k - \theta^{\star}} - 2\gamma\ps{\theta_k - \theta^{\star}}{\nabla U(\theta_k) - \nabla U(\theta^{\star})} + 2\gamma d \\
  &+ \gamma^2 \expecond{ \norm[2]{\nabla U_0(\theta_k) +  \frac{N}{p} \sum_{i\inS_{k+1}} \nabla U_i(\theta_k) }}{\mathcal{F}_k} \\
  &\leq \norm[2]{\theta_k - \theta^{\star}} - 2\gamma\ps{\theta_k - \theta^{\star}}{\nabla U(\theta_k) - \nabla U(\theta^{\star})} + 2\gamma d \\
  &+2\gamma^2\expecond{ \norm[2]{\nabla U_0(\theta_k) - \nabla U_0(\theta^{\star}) +  \frac{N}{p} \sum_{i\inS_{k+1}} \defEns{\nabla U_i(\theta_k) - \nabla U_i(\theta^{\star})} }}{\mathcal{F}_k} \\
  &+2\gamma^2 \expecond{ \norm[2]{\nabla U_0(\theta^{\star}) +  \frac{N}{p} \sum_{i\inS_{k+1}} \nabla U_i(\theta^{\star}) }}{\mathcal{F}_k} \;.
\end{align*}
By \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex}, $\theta\mapsto \nabla U_0(\theta) +  (N/p)\sum_{i\inS} \nabla U_i(\theta)$ is $\mathbb{P}$-\ensuremath{\text{a.s.}}~$L$-co-coercive and we obtain
\begin{multline*}
  \expecond{\norm[2]{\theta_{k+1} - \theta^{\star}}}{\mathcal{F}_k} \leq \defEns{1 - 2m\gamma(1-\gammaL)}\norm[2]{\theta_k - \theta^{\star}} + 2\gamma d \\
  + \frac{2\gamma^2 N}{p} \sum_{i=1}^{N} \norm[2]{\nabla U_i(\theta^{\star}) - \frac{1}{N}\sum_{j=1}^{N} \nabla U_j(\theta^{\star})} \;.
\end{multline*}
A straightforward induction concludes the proof.
\end{proof}

\begin{lemma}\label{lemma:bound-4-tp-tpstar}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  For all initial distribution $\lambda\in\mathcal{P}_4(\mathbb{R}^d)$, $\gamma\in\ocint{0,1/\{12(L \vee 1)\}}$ and $k\in\mathbb{N}$,
  \begin{align*}
    \expe{\norm[4]{\theta_k - \theta^{\star}}} &\leq (1-2m\gamma)^k \int_{\mathbb{R}^d} \norm[4]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) \\
    &+ \defEns{12\gamma^2\expe{\norm[2]{\epsilon(\theta^{\star})}} + 2\gamma(2d+1)} k (1-m\gamma)^{k-1} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) \\
    &+ \defEns{\frac{2d+1}{m} + \frac{6\gamma}{m}\expe{\norm[2]{\epsilon(\theta^{\star})}}}^2 \\
    &+ \frac{2\gamma d(2+d)}{m} + \frac{4 \gamma^3}{m} \expe{\norm[4]{\epsilon(\theta^{\star})}} + \frac{4\gamma^2 (d+2)}{m} \expe{\norm[2]{\epsilon(\theta^{\star})}} \;.
  \end{align*}
  where $(\theta_k)_{k\in\mathbb{N}}$ are the iterates of LMC \eqref{eq:LMC}, SGLD \eqref{eq:SGLD} or SGLDFP \eqref{eq:SGLDFP}.
\end{lemma}

\begin{proof}
Let $\gamma\in\ocint{0,1/\{12(L \vee 1)\}}$, $(\theta_k)_{k\in\mathbb{N}}$ be the iterates of LMC \eqref{eq:LMC}, SGLD \eqref{eq:SGLD} or SGLDFP \eqref{eq:SGLDFP} and $(\mathcal{F}_k)_{k\in\mathbb{N}}$ be the associated filtration. By developing the square, we have
\begin{multline*}
  \norm[4]{\theta_1 - \theta^{\star}} = \Big(\norm[2]{\theta_0 - \theta^{\star}} + 2\gamma\norm[2]{Z_1} + \gamma^2 \norm[2]{\nabla U(\theta_0) + \epsilon_1(\theta_0)}  \\
  - 2\gamma \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{\theta_0 - \theta^{\star}} + \sqrt{2\gamma}\ps{\theta_0 - \theta^{\star}}{Z_1} - (2\gamma)^{3/2} \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{Z_1}\Big)^2 \;,
\end{multline*}
and taking the conditional expectation w.r.t.~$\mathcal{F}_0$,
\begin{align*}
  & \expecond{\norm[4]{\theta_1 - \theta^{\star}}}{\mathcal{F}_0} = \mathbb{E} \Big[ \norm[4]{\theta_0 - \theta^{\star}}+ 4\gamma^2 \norm[4]{Z_1} + \gamma^4 \norm[4]{\nabla U(\theta_0) + \epsilon_1(\theta_0)} \\
  &+ 4 \gamma^2 \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{\theta_0-\theta^{\star}}^2 + 2\gamma \ps{\theta_0 - \theta^{\star}}{Z_1}^2 + (2\gamma)^3 \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{Z_1}^2 \\
  &+ 4\gamma \norm[2]{Z_1}\norm[2]{\theta_0 - \theta^{\star}} + 2 \gamma^2 \norm[2]{\theta_0 - \theta^{\star}} \norm[2]{\nabla U(\theta_0) + \epsilon_1(\theta_0)} \\
  &- 4\gamma \norm[2]{\theta_0 - \theta^{\star}} \ps{\nabla U(\theta_0)}{\theta_0 - \theta^{\star}} + 4\gamma^3 \norm[2]{Z_1}\norm[2]{\nabla U(\theta_0) + \epsilon_1(\theta_0)} \\
  &- 8\gamma^2 \norm[2]{Z_1} \ps{\nabla U(\theta_0)}{\theta_0 - \theta^{\star}} - 4\gamma^3 \norm[2]{\nabla U(\theta_0) + \epsilon_1(\theta_0)} \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{\theta_0 - \theta^{\star}} \\
  &- 8\gamma^2 \ps{\theta_0 - \theta^{\star}}{Z_1} \ps{\nabla U(\theta_0) + \epsilon_1(\theta_0)}{Z_1}
  | \mathcal{F}_0 \Big] \;.
\end{align*}
By \Cref{assum:gradient-Lipschitz} and \Cref{assum:Ui-convex}, $\theta\mapsto \nabla U(\theta) +  \epsilon_1(\theta)$ is $\mathbb{P}$-\ensuremath{\text{a.s.}}~$L$-co-coercive and we have for all $\theta\in\mathbb{R}^d$, $\mathbb{P}$-\ensuremath{\text{a.s.}}~,
\begin{align*}
  \norm[2]{\nabla U(\theta) + \epsilon_1(\theta) - \epsilon_1(\theta^{\star})} & \leq L \ps{\theta - \theta^{\star}}{\nabla U(\theta) + \epsilon_1(\theta) - \epsilon_1(\theta^{\star})} \;, \\
  \norm[4]{\nabla U(\theta) + \epsilon_1(\theta) - \epsilon_1(\theta^{\star})} & \leq L^2\norm[2]{\theta-\theta^{\star}} \ps{\theta - \theta^{\star}}{\nabla U(\theta) + \epsilon_1(\theta) - \epsilon_1(\theta^{\star})} \;.
\end{align*}
Combining it with $\expe{\norm[4]{Z_1}} = d(2+d)$, we obtain
\begin{align*}
  & \expecond{\norm[4]{\theta_1 - \theta^{\star}}}{\mathcal{F}_0,S_1} \leq \norm[4]{\theta_0 - \theta^{\star}} - 4\gamma(1 - 3\gammaL - 2\gamma^3L^2)\norm[2]{\theta_0 - \theta^{\star}} \\
  &\times \ps{\theta_0 - \theta^{\star}}{\nabla U(\theta_0) + \epsilon_1(\theta_0) - \epsilon_1(\theta^{\star})} + (12\gamma^2\norm[2]{\epsilon_1(\theta^{\star})} + 2\gamma(2d+1)) \norm[2]{\theta_0 - \theta^{\star}} \\
  &+ 4\gamma^2 d(2+d) + 8 \gamma^4 \norm[4]{\epsilon_1(\theta^{\star})} + 8\gamma^3 (d+2) \norm[2]{\epsilon_1(\theta^{\star})} \\
  &- 8(d+1)\gamma^2(1-2\gammaL) \ps{\theta_0 - \theta^{\star}}{\nabla U(\theta_0) + \epsilon_1(\theta_0) - \epsilon_1(\theta^{\star})} \;.
\end{align*}
By \Cref{assum:strongly-convex} and using $\gamma\leq 1/\{12(L \vee 1)\}$, we get
\begin{align*}
  &\expecond{\norm[4]{\theta_1 - \theta^{\star}}}{\mathcal{F}_0} \leq (1-2m\gamma)\norm[4]{\theta_0 - \theta^{\star}} + \defEns{12\gamma^2\expe{\norm[2]{\epsilon_1(\theta^{\star})}} + 2\gamma(2d+1)} \norm[2]{\theta_0 - \theta^{\star}} \\
  &+ 4\gamma^2 d(2+d) + 8 \gamma^4 \expe{\norm[4]{\epsilon_1(\theta^{\star})}} + 8\gamma^3 (d+2) \expe{\norm[2]{\epsilon_1(\theta^{\star})}} \;.
\end{align*}
By a straightforward induction, we have for all $n\in\mathbb{N}$
\begin{align*}
  &\expe{\norm[4]{\theta_n - \theta^{\star}}} \leq (1-2m\gamma)^n \expe{\norm[4]{\theta_0 - \theta^{\star}}} \\
  &+ \defEns{12\gamma^2\expe{\norm[2]{\epsilon(\theta^{\star})}} + 2\gamma(2d+1)} \sum_{k=0}^{n-1}(1-2m\gamma)^{n-1-k} \expe{\norm[2]{\theta_k - \theta^{\star}}} \\
  &+ (2m\gamma)^{-1}\defEns{4\gamma^2 d(2+d) + 8 \gamma^4 \expe{\norm[4]{\epsilon(\theta^{\star})}} + 8\gamma^3 (d+2) \expe{\norm[2]{\epsilon(\theta^{\star})}}}
\end{align*}
and by \Cref{lemma:bound-tp-tpstar},
\begin{align*}
  \expe{\norm[4]{\theta_n - \theta^{\star}}} &\leq (1-2m\gamma)^n \int_{\mathbb{R}^d} \norm[4]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) \\
  &+ \defEns{12\gamma^2\expe{\norm[2]{\epsilon_1(\theta^{\star})}} + 2\gamma(2d+1)} n (1-m\gamma)^{n-1} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \lambda(\mathrm{d} \vartheta) \\
  &+ \defEns{\frac{2d+1}{m} + \frac{6\gamma}{m}\expe{\norm[2]{\epsilon(\theta^{\star})}}}^2 \\
  &+ \frac{2\gamma d(2+d)}{m} + \frac{4 \gamma^3}{m} \expe{\norm[4]{\epsilon(\theta^{\star})}} + \frac{4\gamma^2 (d+2)}{m} \expe{\norm[2]{\epsilon(\theta^{\star})}} \;.
\end{align*}
\end{proof}

Thanks to this lemma, we obtain the following corollary. The upper bound for SGD is given by \cite[Lemma 13]{dieuleveut:durmus:bach:2017}.

\begin{corollary}\label{corollary:moments-4-sgld-sgldfp-ula}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  \begin{enumerate}[label=\roman*)]
    \item \label{item:moments-4-1}
      Let $\gamma = \eta/N$ with $\eta\in\ocintLigne{0,1/\{24(\tilde{L} \vee 1)\}}$ and assume that $\liminf_{N\to+\infty} N^{-1} m >0$. Then,
      \begin{align*}
       & \int_{\mathbb{R}^d} \norm[4]{\theta - \theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{LMC}}} (\mathrm{d} \theta) = d^2 O_{N\to+\infty}(N^{-2}) \;, \\
       &\int_{\mathbb{R}^d} \norm[4]{\theta - \theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{FP}}} (\mathrm{d} \theta) = d^2 O_{N\to+\infty}(N^{-2}) \;.
       \end{align*}
    \item \label{item:moments-4-2}
      Let $\gamma = \eta/N$ with $\eta\in\ocintLigne{0,1/\{24(\tilde{L} \vee 1)\}}$ and assume that $\liminf_{N\to+\infty} N^{-1} m >0$ and that $N \geq 1/\eta$. Then,
      \begin{align*}
      &\int_{\mathbb{R}^d} \norm[4]{\theta - \theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{SGLD}}} (\mathrm{d} \theta) = d^2 O_{\eta\to 0}(\eta^2) \;, \;
      \int_{\mathbb{R}^d} \norm[4]{\theta - \theta^{\star}} \pi_{\scriptscriptstyle{\operatorname{SGD}}} (\mathrm{d} \theta) = d^2 O_{\eta \to 0}(\eta^2) \;.
      \end{align*}
  \end{enumerate}
\end{corollary}

\subsubsection{Proofs of \Cref{thm:bias-variance-tpstar-ula-fp} and \Cref{thm:bias-variance-tpstar-sgld-sgd}}
\label{subsec:proof-thm-bias-variance-tpstar}

Denote by
\begin{equation}\label{eq:def-eta-0}
  \eta_0 = \inf_{N \geq 1} \defEns{\frac{N}{12(L \vee 1)} \wedge \frac{2N}{(1+r^2) L}} > 0 \;,
\end{equation}
and set $\gamma=\eta/N$ with $\eta\in\ooint{0,\eta_0}$.
Let $\updelta\in\defEns{0,1}$ be equal to $1$ for LMC, SGLDFP and SGLD and $0$ for SGD.
Let $\theta_0$ be distributed according to $\bar{\pi}$. By \eqref{eq:update-gesp} and using a Taylor expansion around $\theta^{\star}$ for $\nabla U$, we obtain
\begin{equation*}
  \theta_1 - \theta^{\star} = \theta_0 - \theta^{\star} - \gamma\parenthese{\nabla^2 U(\theta^{\star}) (\theta_0 - \theta^{\star}) + \reste{1}(\theta_0) + \epsilon_1(\theta_0)} + \updelta \sqrt{2\gamma} Z_1 \;,
\end{equation*}
where by \Cref{assum:gradient-Lipschitz}, $\reste{1}:\mathbb{R}^d \to \mathbb{R}^d$ satisfies
\begin{equation}\label{eq:proof-bias-variance-reste-1}
  \sup_{\theta\in\mathbb{R}^d} \defEns{\norm{\reste{1}(\theta)} / \norm[2]{\theta - \theta^{\star}}} \leq L/2 \;.
\end{equation}
Taking the tensor product and the expectation, and using that $\theta_0,\epsilon_1,Z_1$ are mutually independent, we obtain
\begin{multline}
  \label{eq:bias-variance-tpstar-ula-fp_1}
  \operatorname{H} \expe{(\theta_0 - \theta^{\star})^{\otimes 2}} = 2\updelta \operatorname{Id} + \gamma \expe{\epsilon_1(\theta_0)^{\otimes 2}} + \expe{\reste{1}(\theta_0) \otimes \{ \theta_0-\theta^{\star}\} + \{ \theta_0-\theta^{\star}\} \otimes \reste{1}(\theta_0)} \\
  + \gamma \expe{ \reste{1}(\theta_0)^{\otimes 2} + \{\nabla^2 U(\theta^{\star}) (\theta_0 - \theta^{\star})\} \otimes \reste{1}(\theta_0)+ \reste{1}(\theta_0)\otimes \nabla^2 U(\theta^{\star}) (\theta_0 - \theta^{\star})} \;.
\end{multline}
For LMC, $\epsilon_1$ is the null function and by \Cref{corollary:moments-4-sgld-sgldfp-ula}-\ref{item:moments-4-1}, \eqref{eq:proof-bias-variance-reste-1} and \eqref{eq:bias-variance-tpstar-ula-fp_1}, we obtain \eqref{eq:covariance-lmc}.
Regarding SGLDFP, SGLD and SGD, by a Taylor expansion of $\epsilon_1$ around $\theta^{\star}$, we get for all $\theta\in\mathbb{R}^d$, $\mathbb{P}$-\ensuremath{\text{a.s.}}~,
\begin{equation*}
  \epsilon_1(\theta) = \epsilon_1(\theta^{\star}) + \nabla \epsilon_1(\theta^{\star})(\theta - \theta^{\star}) + \reste{2}(\theta)
\end{equation*}
where by \Cref{assum:gradient-Lipschitz},$\reste{2}:\mathbb{R}^d \to \mathbb{R}^d$ satisfies
\begin{equation}\label{eq:proof-bias-variance-reste-2}
  \sup_{\theta\in\mathbb{R}^d} \defEns{\norm{\reste{2}(\theta)} / \norm[2]{\theta - \theta^{\star}}} \leq L/2 \;.
\end{equation}
Therefore, taking the tensor product and the expectation, we obtain
\begin{equation}\label{eq:proof-development-gesp}
  \expe{\epsilon_1(\theta_0)^{\otimes 2}} = \expe{\epsilon_1(\theta^{\star})^{\otimes 2}} + \parenthese{\nabla \epsilon_1(\theta^{\star})}^{\otimes 2}\expe{(\theta_0 - \theta^{\star})^{\otimes 2}} + \expe{\reste{3}(\theta_0)}
\end{equation}
where $\reste{3}:\mathbb{R}^d \to \mathbb{R}^{d \times d}$ is defined for all $\theta\in\mathbb{R}^d$, $\mathbb{P}$-\ensuremath{\text{a.s.}}~,
\begin{align}
  \nonumber
  &\reste{3}(\theta) = \epsilon_1(\theta^{\star}) \otimes \{ \nabla \epsilon_1(\theta^{\star})(\theta - \theta^{\star})\} + \{ \nabla \epsilon_1(\theta^{\star})(\theta - \theta^{\star}) \} \otimes \epsilon_1(\theta^{\star}) \\
  \label{eq:proof-bias-variance-reste-3}
  &+ \{ \epsilon_1(\theta^{\star}) + \nabla \epsilon_1(\theta^{\star})(\theta - \theta^{\star})\} \otimes \reste{2}(\theta) + \reste{2}(\theta) \otimes \{ \epsilon_1(\theta^{\star}) + \nabla \epsilon_1(\theta^{\star})(\theta - \theta^{\star})\} + \reste{2}^{\otimes 2}(\theta) \;.
\end{align}
Note that $\opL = \expe{\parenthese{\nabla \epsilon_1(\theta^{\star})}^{\otimes 2}}$.
For SGLDFP, $\epsilon_1(\theta^{\star})=0$ \ensuremath{\text{a.s.}}~By \Cref{corollary:moments-4-sgld-sgldfp-ula}-\ref{item:moments-4-1}, \eqref{eq:proof-bias-variance-reste-1}, \eqref{eq:bias-variance-tpstar-ula-fp_1}, \eqref{eq:proof-bias-variance-reste-2}, \eqref{eq:proof-development-gesp} and \eqref{eq:proof-bias-variance-reste-3}, we obtain \eqref{eq:covariance-fp}.

Regarding SGLD and SGD, we have $\expe{\epsilon_1(\theta^{\star})^{\otimes 2}} = (N/p) \operatorname{M}$ where $\operatorname{M}$ is defined in \eqref{eq:def-opM}. By \Cref{corollary:moments-4-sgld-sgldfp-ula}-\ref{item:moments-4-2}, \eqref{eq:proof-bias-variance-reste-1}, \eqref{eq:bias-variance-tpstar-ula-fp_1}, \eqref{eq:proof-bias-variance-reste-2}, \eqref{eq:proof-development-gesp} and \eqref{eq:proof-bias-variance-reste-3}, we obtain \eqref{eq:covariance-sgld} and \eqref{eq:covariance-sgd}.

For the mean of $\pi_{\scriptscriptstyle{\operatorname{LMC}}},\pi_{\scriptscriptstyle{\operatorname{FP}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$, by a Taylor expansion around $\theta^{\star}$ for $\nabla U$ of order $3$, we obtain
\begin{multline*}
  \theta_1 - \theta^{\star} = \theta_0 - \theta^{\star} - \gamma\parenthese{\nabla^2 U(\theta^{\star}) (\theta_0 - \theta^{\star})+ (1/2)\operatorname{D}^3 U(\theta^{\star}) (\theta_0 - \theta^{\star})^{\otimes 2} + \reste{4}(\theta_0) + \epsilon_1(\theta_0)} \\
  + \updelta \sqrt{2\gamma} Z_1 \;,
\end{multline*}
where  by \Cref{assum:gradient-Lipschitz}, $\reste{4}:\mathbb{R}^d \to \mathbb{R}^d$ satisfies
\begin{equation}\label{eq:proof-bias-variance-reste-4}
  \sup_{\theta\in\mathbb{R}^d} \defEns{\norm{\reste{4}(\theta)} / \norm[3]{\theta - \theta^{\star}}} \leq L/6 \;.
\end{equation}
Taking the expectation and using that $\theta_1$ is distributed according to $\bar{\pi}$, we get
\begin{equation*}
  \expe{\theta_0} - \theta^{\star} =-(1/2) \nabla^2 U(\theta^{\star}) \operatorname{D}^3 U(\theta^{\star})[\expe{(\theta_0 - \theta^{\star})^{\otimes 2}}] - \nabla^2 U(\theta^{\star})^{-1} \expe{\reste{4}(\theta_0)} \;.
\end{equation*}
\eqref{eq:covariance-lmc}, \eqref{eq:covariance-fp}, \eqref{eq:covariance-sgld},\eqref{eq:covariance-sgd}, \eqref{eq:proof-bias-variance-reste-4} and \Cref{corollary:moments-4-sgld-sgldfp-ula} conclude the proof.

\section{Results}
\label{sec:results}

\subsection{Analysis in Wasserstein distance}
\label{subsec:wasserstein-distance}

Before presenting the results, some notations and elements of Markov chain theory have to be introduced. Denote by $\mathcal{P}_2(\mathbb{R}^d)$ the set of probability measures  with finite second moment and by $\mathcal{B}(\mathbb{R}^d)$ the Borel $\sigma$-algebra of $\mathbb{R}^d$. For  $\lambda,\nu\in\mathcal{P}_2(\mathbb{R}^d)$, define the Wasserstein distance of order $2$ by
\begin{equation*}
  \wasser[2](\lambda,\nu) = \inf_{\xi\in\Pi(\lambda,\nu)} \parenthese{\int_{\mathbb{R}^d \times \mathbb{R}^d} \norm[2]{\theta - \vartheta} \xi(\mathrm{d} \theta, \mathrm{d} \vartheta)}^{1/2} \;,
\end{equation*}
where $\Pi(\lambda,\nu)$ is the set of probability measures $\xi$ on $\mathcal{B}(\mathbb{R}^d) \otimes \mathcal{B}(\mathbb{R}^d)$ satisfying for all $\mathsf{A}\in\mathcal{B}(\mathbb{R}^d)$, $\xi(\mathsf{A} \times \mathbb{R}^d)) = \lambda(\mathsf{A})$ and $\xi(\mathbb{R}^d \times \mathsf{A}) = \nu(\mathsf{A})$.

A Markov kernel $R$ on $\mathbb{R}^d \times \mathcal{B}(\mathbb{R}^d)$ is a mapping $R: \mathbb{R}^d \times \mathcal{B}(\mathbb{R}^d) \to \ccint{0,1}$ satisfying the following conditions: (i) for every $\theta \in \mathbb{R}^d$, $R(\theta,\cdot): \mathsf{A} \mapsto R(\theta,\mathsf{A})$ is a probability measure on $\mathcal{B}(\mathbb{R}^d)$ (ii) for every $\mathsf{A} \in \mathcal{B}(\mathbb{R}^d)$, $R(\cdot,A): \theta\mapsto  R(\theta,A)$ is a measurable function.
For any probability measure $\lambda$ on $\mathcal{B}(\mathbb{R}^d)$, we define $\lambda R$ for all $\mathsf{A}\in\mathcal{B}(\mathbb{R}^d)$ by $\lambdaR(\mathsf{A}) = \int_{\mathbb{R}^d} \lambda(\mathrm{d} \theta) R(\theta,\mathsf{A})$.
For all $k\in\mathbb{N}^*$, we define the Markov kernel $R^k$ recursively by $R^1 = R$ and for all $\theta\in\mathbb{R}^d$ and $\mathsf{A}\in\mathcal{B}(\mathbb{R}^d)$,
$  R^{k+1}(\theta,\mathsf{A}) = \int_{\mathbb{R}^d} R^k(\theta, \mathrm{d} \vartheta) R(\vartheta, \mathsf{A}) $.
A probability measure $\bar{\pi}$  is invariant for $R$ if $\bar{\pi} R = \bar{\pi}$.

The  LMC, SGLD, SGD and SGLDFP algorithms defined respectively by \eqref{eq:LMC}, \eqref{eq:SGLD}, \eqref{eq:SGD} and \eqref{eq:SGLDFP} are homogeneous Markov chains with Markov kernels denoted $\Rker_{\scriptscriptstyle{\operatorname{LMC}}},\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}, \Rker_{\scriptscriptstyle{\operatorname{SGD}}}$, and $\Rker_{\scriptscriptstyle{\operatorname{FP}}}$. To avoid overloading the notations, the dependence on $\gamma$ and  $N$ is implicit.
\begin{lemma}\label{lemma:existence-pi-gamma}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  For any step size $\gamma\in\ooint{0,2/L}$, $\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}$ (respectively $\Rker_{\scriptscriptstyle{\operatorname{LMC}}},\Rker_{\scriptscriptstyle{\operatorname{SGD}}},\Rker_{\scriptscriptstyle{\operatorname{FP}}}$) has a unique invariant measure $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}\in\mathcal{P}_2(\mathbb{R}^d)$ (respectively $\pi_{\scriptscriptstyle{\operatorname{LMC}}},\pi_{\scriptscriptstyle{\operatorname{SGD}}},\pi_{\scriptscriptstyle{\operatorname{FP}}}$). In addition, for all $\gamma\in\ocint{0,1/L}$, $\theta\in\mathbb{R}^d$ and $k\in\mathbb{N}$,
  \begin{equation*}
    \wasser[2]^2(\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^k(\theta,\cdot), \pi_{\scriptscriptstyle{\operatorname{SGLD}}}) \leq (1-m\gamma)^k \int_{\mathbb{R}^d} \norm[2]{\theta - \vartheta} \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \vartheta)
  \end{equation*}
  and the same inequality holds for LMC, SGD and SGLDFP.
\end{lemma}

\begin{proof}
  The proof is postponed to \Cref{subsec:proof-lemma-existence-pi-gamma}.
\end{proof}

Under \Cref{assum:gradient-Lipschitz}, \eqref{eq:langevin-SDE} has a unique strong solution $(\theta_t)_{t\geq 0}$ for every initial condition $\theta_0\in\mathbb{R}^d$ \cite[Chapter 5, Theorems 2.5 and 2.9]{karatzas:shreve:1991}. Denote by $(P_t)_{t\geq 0}$ the semigroup of the Langevin diffusion defined for all $\theta_0 \in \mathbb{R}^d$ and $\mathsf{A} \in \mcb{\mathbb{R}^d}$ by $P_t(\theta_0,\mathsf{A}) = \mathbb{P}(\theta_t \in \mathsf{A})$.

\begin{theorem}\label{thm:wasserstein-ula-sgld-sgd-fp}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  For all $\gamma\in\ocint{0,1/L}$, $\lambda,\mu\in\mathcal{P}_2(\mathbb{R}^d)$ and $n\in\mathbb{N}$, we have the following upper-bounds in Wasserstein distance between
  \begin{enumerate}[label=\roman*)]
    \item \label{item:W-bound-ULA-SGLDFP}
    LMC and SGLDFP,
    \begin{multline*}
      \wasser[2]^2(\lambda \Rker_{\scriptscriptstyle{\operatorname{LMC}}}^n, \mu \Rker_{\scriptscriptstyle{\operatorname{FP}}}^n) \leq (1-m\gamma)^n \wasser[2]^2(\lambda,\mu) + \frac{2L^2\gamma d}{p m^2} \\
      + \frac{L^2\gamma^2}{p} n (1-m\gamma)^{n-1} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \mu(\mathrm{d} \vartheta) \;,
    \end{multline*}
    \item \label{item:W-bound-ULA-Langevin}
    the Langevin diffusion and LMC,
    \begin{align*}
      \wasser[2]^2(\lambda \Rker_{\scriptscriptstyle{\operatorname{LMC}}}^n, \mu P_{n\gamma}) &\leq 2\parenthese{1 - \frac{\mUL\gamma}{m + L}}^n \wasser[2]^2(\lambda,\mu) + d \gamma \frac{m + L}{2m}\parenthese{3+\frac{L}{m}} \parenthese{\frac{13}{6}+\frac{L}{m}} \\
      &+  n \mathrm{e}^{-(m/2)\gamma(n-1)} L^3 \gamma^3 \parenthese{1+\frac{m+L}{2m}} \int_{\mathbb{R}^d} \norm[2]{\vartheta - \theta^{\star}} \mu(\mathrm{d} \vartheta) \;,
    \end{align*}
    \item \label{item:W-bound-SGLD-SGD}
    SGLD and SGD,
    \begin{equation*}
      \wasser[2]^2(\lambda \Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^n, \mu \Rker_{\scriptscriptstyle{\operatorname{SGD}}}^n) \leq (1-m\gamma)^n \wasser[2]^2(\lambda,\mu) + (2d)/m \;.
    \end{equation*}
  \end{enumerate}
\end{theorem}

\begin{proof}
  The proof is postponed to \Cref{subsec:proof-thm-wasserstein-ula-sgld-sgd-fp}.
\end{proof}

\begin{corollary}\label{corollary-wasserstein-piula-pisgld-pisgd-pifp}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}. Set $\gamma = \eta / N$ with $\eta \in \ocintLigne{0,1/(2\tilde{L})}$ and assume that $\liminf_{N \to \infty}  m N^{-1}>0$. Then,

\begin{enumerate}[label=\roman*)]
\item   for all $n \in \mathbb{N}$,  we get $ \wasser[2](\Rker_{\scriptscriptstyle{\operatorname{LMC}}}^n(\theta^{\star},\cdot), \Rker_{\scriptscriptstyle{\operatorname{FP}}}^n(\theta^{\star},\cdot)) = \sqrt{d \eta } \, O(N^{-1/2})$ and $
    \wasser[2](\pi_{\scriptscriptstyle{\operatorname{LMC}}}, \pi_{\scriptscriptstyle{\operatorname{FP}}}) =  \sqrt{d \eta }\, O(N^{-1/2})$, \, $ \wasser[2](\pi_{\scriptscriptstyle{\operatorname{LMC}}}, \pi) = \sqrt{d \eta }\, O(N^{-1/2})$.
    \item for all $n \in \mathbb{N}$,  we get $ \wasser[2](\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^n(\theta^{\star},\cdot), \Rker_{\scriptscriptstyle{\operatorname{SGD}}}^n(\theta^{\star},\cdot)) = \sqrt{d} \, O(N^{-1/2})$ and $ \wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}}, \pi_{\scriptscriptstyle{\operatorname{SGD}}}) = \sqrt{d}\, O(N^{-1/2}) $.
\end{enumerate}
  \end{corollary}

\Cref{thm:wasserstein-ula-sgld-sgd-fp} implies that the number of iterations necessary to obtain a sample $\varepsilon$-close from $\pi$ in Wasserstein distance is the same for LMC and SGLDFP.
However for LMC, the cost of one iteration is $N d$ which is larger than $p d$ the cost of one iteration for SGLDFP. In other words, to obtain an approximate sample from the target distribution at an accuracy $O(1/\sqrt{N})$ in $2$-Wasserstein distance, LMC requires $\Theta(N)$ operations, in contrast with SGLDFP that needs only $\Theta(1)$ operations.

 We show in the sequel that $\wasser(\pi_{\scriptscriptstyle{\operatorname{FP}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}) = \Omega(1)$ when $N\to+\infty$ in the case of a Bayesian linear regression, where for two sequences $(u_N)_{N \geq 1}$, $(v_N)_{N \geq 1}$, $u_N = \Omega(v_N)$ if $\liminf_{N \to +\infty} u_N/v_N >0$.
The dataset is $\mathbf{z} = \{(y_i,x_i)\}_{i=1}^{N}$ where $y_i\in\mathbb{R}$ is the response variable and  $x_i\in\mathbb{R}^d$ are the covariates. Set $\mathbf{y} = (y_1,\ldots,y_N) \in\mathbb{R}^N$ and $\mathbf{X} \in\mathbb{R}^{N\times d}$ the matrix of covariates such that the $i^{\text{th}}$ row of $\mathbf{X}$ is $x_i$. Let $\sigma^2_y, \sigma^2_{\tp}>0$. For $i\in\defEns{1,\ldots,N}$, the conditional distribution of $y_i$ given $x_i$ is Gaussian with mean $x_i^{\operatorname{T}} \theta$ and variance $\sigma^2_y$.
The prior $\pi_0(\theta)$ is a normal distribution of mean $0$ and variance $\sigma^2_{\tp} \operatorname{Id}$. The posterior distribution $\pi$ is then proportional to $\pi(\theta) \propto \exp\parenthese{-(1/2)(\theta-\theta^{\star})^{\operatorname{T}} \Sigma (\theta-\theta^{\star})}$ where
\begin{equation*}
  \Sigma = \operatorname{Id}/\sigma^2_{\tp} +\mathbf{X}^{\operatorname{T}}\mathbf{X}/\sigma^2_y \quad \text{and} \quad \theta^{\star} = \Sigma^{-1} (\mathbf{X}^{\operatorname{T}} \mathbf{y})/\sigma^2_y \;.
\end{equation*}
We assume that $\mathbf{X}^{\operatorname{T}}\mathbf{X} \succeq m \operatorname{Id}$, with $\liminf_{N\to +\infty} m/N >0$.  Let $S$ be a minibatch of $\defEns{1,\ldots,N}$ with replacement of size $p$. Define
\begin{equation*}
  \nabla U_0(\theta) + (N/p)\sum_{i\inS} \nabla U_i(\theta) = \Sigma(\theta-\theta^{\star}) + \rho(S)(\theta-\theta^{\star}) + \xi(S)
\end{equation*}
where
\begin{equation}\label{eq:def-partm-parta}
  \rho(S) = \frac{\operatorname{Id}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y}\sum_{i\inS} x_i x_i^{\operatorname{T}} - \Sigma \;,\;
  \xi(S) = \frac{\theta^{\star}}{\sigma^2_{\tp}} + \frac{N}{p \sigma^2_y} \sum_{i\inS} \parenthese{x_i^{\operatorname{T}} \theta^{\star} - y_i} x_i \;.
\end{equation}
$\rho(S)(\theta-\theta^{\star})$ is the multiplicative part of the noise in the stochastic gradient, and $\xi(S)$ the additive part that does not depend on $\theta$. The additive part of the stochastic gradient for SGLDFP disappears since
\begin{equation*}
  \nabla U_0(\theta) - \nabla U_0(\theta^{\star}) + (N/p)\sum_{i\inS} \defEns{\nabla U_i(\theta) - \nabla U_i(\theta^{\star})}= \Sigma(\theta-\theta^{\star}) + \rho(S)(\theta-\theta^{\star}) \;.
\end{equation*}

In this setting, the following theorem shows that the Wasserstein distances between the marginal distribution of the iterates of SGLD and SGLDFP, and $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi$, is of order $\Omega(1)$ when $N\to+\infty$. This is in sharp contrast with the results of \Cref{corollary-wasserstein-piula-pisgld-pisgd-pifp} where the Wasserstein distances tend to $0$ as $N\to+\infty$ at a rate $N^{-1/2}$. For simplicity, we state the result for $d=1$.

\begin{theorem}\label{thm:W-lower-bound-SGLD-SGLDFP}
  Consider the case of the Bayesian linear regression in dimension $1$.
  \begin{enumerate}[label=\roman*)]
    \item \label{item:W-lower-bound-SGLD-FP-1}
     For all $\gamma\in\ocintLigne{0,\Sigma^{-1}\{1+N/(p\sum_{i=1}^{N}x_i^2)\}^{-1}}$ and $n\in\mathbb{N}^*$,
    \begin{multline*}
      \parenthese{\frac{1-\mu}{1-\mu^n}}^{1/2} \wasser[2](\Rker_{\scriptscriptstyle{\operatorname{SGLD}}}^n(\theta^{\star},\cdot), \Rker_{\scriptscriptstyle{\operatorname{FP}}}^n(\theta^{\star},\cdot)) \\
      \geq \defEns{2\gamma + \frac{\gamma^2 N}{p} \sum_{i=1}^{N}\parenthese{\frac{(x_i \theta^{\star} - y_i)x_i}{\sigma^2_y} + \frac{\theta^{\star}}{N\sigma^2_{\tp}}}^2}^{1/2} - \sqrt{2\gamma} \;,
    \end{multline*}
    where $\mu \in\ocint{0, 1 - \gamma\Sigma}$.
    \item \label{item:W-lower-bound-SGLD-FP-2}
    Set $\gamma = \eta/N$ with $\eta\in\ocintLigne{0, \liminf_{N\to +\infty} N\Sigma^{-1}\{1+N/(p\sum_{i=1}^{N}x_i^2)\}^{-1}}$ and assume that $\liminf_{N\to+\infty} N^{-1}\sum_{i=1}^{N} x_i^2 >0$. We have $\wasser[2](\pi_{\scriptscriptstyle{\operatorname{SGLD}}},\pi) = \Omega(1)$.

  \end{enumerate}
\end{theorem}

\begin{proof}
  The proof is postponed to \Cref{subsec:proof-thm-W-lower-bound-SGLD-SGLDFP}.
\end{proof}

The study in Wasserstein distance emphasizes the different behaviors of the LMC, SGLDFP, SGLD and SGD algorithms. When $N \to \infty$ and $\lim_{N \to +\infty} m/N>0$, the marginal distributions of the $k^{\text{th}}$ iterates of the  LMC and SGLDFP algorithm are very close to the Langevin diffusion and their invariant probability measures $\pi_{\scriptscriptstyle{\operatorname{LMC}}}$ and $\pi_{\scriptscriptstyle{\operatorname{FP}}}$ are similar to the posterior distribution of interest $\pi$. In contrast, the marginal distributions of the $k^{\text{th}}$ iterates  of SGLD and SGD are analogous and their invariant probability measures $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ are very different from $\pi$ when $N\to+\infty$.

Note that to fix the asymptotic bias of SGLD, other strategies can be considered: choosing a step size $\gamma \propto N^{-\beta}$ where $\beta>1$ and/or increasing the batch size $p \propto N^{\alpha}$ where $\alpha\in\ccint{0,1}$. Using the Wasserstein (of order 2) bounds of SGLD w.r.t.~the target distribution $\pi$, see e.g.~\cite[Theorem 3]{dalalyan:karagulyan:2017}, $\alpha + \beta$ should be equal to $2$ to guarantee the $\varepsilon$-accuracy in Wasserstein distance of SGLD for a cost proportional to $N$ (up to logarithmic terms), independently of the choice of $\alpha$ and $\beta$.

\subsection{Mean and covariance matrix of $\pi_{\scriptscriptstyle{\operatorname{LMC}}}, \pi_{\scriptscriptstyle{\operatorname{FP}}}, \pi_{\scriptscriptstyle{\operatorname{SGLD}}}$}
\label{subsec:variance-piula-pifp-puslgd-pisgd}

 We now establish an expansion of the mean and second moments of  $\pi_{\scriptscriptstyle{\operatorname{LMC}}},\pi_{\scriptscriptstyle{\operatorname{FP}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ as $N \to +\infty$, and  compare them.  We first give an expansion of the mean and second moments of $\pi$ as $N \to +\infty$.

\begin{proposition}\label{lemma:mean-pi-tpstar}
  Assume \Cref{assum:gradient-Lipschitz} and \Cref{assum:strongly-convex} and that $\liminf_{N\to+\infty} N^{-1} m  >0$. Then,
  \begin{align*}
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi(\mathrm{d} \theta) &= \nabla^2 U(\theta^{\star})^{-1} + O_{N\to+\infty}(N^{-3/2}) \;, \\
    \int_{\mathbb{R}^d} \theta \,  \pi(\mathrm{d} \theta) - \theta^{\star} &= - (1/2) \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star})[\nabla^2 U(\theta^{\star})^{-1}] + O_{N\to+\infty}(N^{-3/2}) \;.
  \end{align*}

\end{proposition}

\begin{proof}
  The proof is postponed to \Cref{subsubsec:proof-lemma-mean-pi-tpstar}.
\end{proof}

Contrary to the Bayesian linear regression  where the covariance matrices can be explicitly computed, see \Cref{sec:mean_var_gaussian}, only approximate expressions are available in the general case. For that purpose, we consider two types of asymptotic. For LMC and SGLDFP, we assume that $\lim_{N \to +\infty} m/N >0$, $\gamma = \eta /N$, for $\eta >0$, and we develop an asymptotic when $N\to+\infty$. Combining \Cref{lemma:mean-pi-tpstar} and \Cref{thm:bias-variance-tpstar-ula-fp} , we show that the biases and covariance matrices of $\pi_{\scriptscriptstyle{\operatorname{LMC}}}$ and $\pi_{\scriptscriptstyle{\operatorname{FP}}}$ are of order $\Theta(1/N)$ with remainder terms of the form $O(N^{-3/2})$, where for two sequences $(u_N)_{N \geq 1}$, $(v_N)_{N \geq 1}$, $u=\Theta(v)$ if $0< \liminf_{N \to +\infty} u_N / v_N  \leq \limsup_{N \to +\infty} u_N / v_N  < +\infty$.

Regarding SGD and SGLD, we do not have such concentration properties when $N\to+\infty$ because of the high variance of the stochastic gradients. The biases and covariance matrices of SGLD and SGD are of order $\Theta(1)$ when $N\to+\infty$. To obtain approximate expressions of these quantities, we set $\gamma = \eta/N$ where $\eta>0$ is the step size for the gradient descent over the normalized potential $U/N$. Assuming that $m$ is proportional to $N$ and $N \geq 1/\eta$, we show by combining \Cref{lemma:mean-pi-tpstar} and \Cref{thm:bias-variance-tpstar-sgld-sgd} that the biases and covariance matrices of SGLD and SGD are of order $\Theta(\eta)$ with remainder terms of the form $O(\eta^{3/2})$ when $\eta\to 0$.

Before giving the results associated to $\pi_{\scriptscriptstyle{\operatorname{LMC}}},\pi_{\scriptscriptstyle{\operatorname{FP}}},\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$, we need to introduce some notations. For any matrices $A_1,A_2\in\mathbb{R}^{d\times d}$, we denote by $A_1 \otimes A_2$ the Kronecker product defined on $\mathbb{R}^{d\times d}$ by $A_1 \otimes A_2 : Q \mapsto A_1 Q A_2$ and  $A^{\otimes 2} = A \otimes A$. Besides, for all $\theta_1\in\mathbb{R}^d$ and $\theta_2\in\mathbb{R}^d$, we denote by $\theta_1 \otimes \theta_2 \in \mathbb{R}^{d\times d}$ the tensor product of $\theta_1$ and $\theta_2$. For any matrix $A\in\mathbb{R}^{d\times d}$, $\operatorname{Tr}(A)$ is the trace of $A$.

Define $\opL:\mathbb{R}^{d\times d}\to \mathbb{R}^{d\times d}$ for all $A\in\mathbb{R}^{d\times d}$ by
\begin{equation}\label{eq:def-operator-L}
  \opL (A) = \frac{N}{p}\sum_{i=1}^{N} \parenthese{\nabla^2 U_i(\theta^{\star}) - \frac{1}{N}\sum_{j=1}^{N} \nabla^2 U_j(\theta^{\star})}^{\otimes 2} A \;.

\end{equation}
and $\operatorname{H}$ and $\operatorname{G}:\mathbb{R}^{d \times d} \to \mathbb{R}^{d \times d}$ by
\begin{align}
  \label{eq:def-operator-H}
  \operatorname{H} &= \nabla^2 U(\theta^{\star}) \otimes \operatorname{Id} + \operatorname{Id} \otimes \nabla^2 U(\theta^{\star}) - \gamma \nabla^2 U(\theta^{\star}) \otimes \nabla^2 U(\theta^{\star}) \;, \\
  \label{eq:def-operator-G}
  \operatorname{G} &= \nabla^2 U(\theta^{\star}) \otimes \operatorname{Id} + \operatorname{Id} \otimes \nabla^2 U(\theta^{\star}) - \gamma (\nabla^2 U(\theta^{\star}) \otimes \nabla^2 U(\theta^{\star}) + \opL) \;.
\end{align}
$\opL$, $\operatorname{H}$ and $\operatorname{G}$ can be interpreted as perturbations of $\nabla^2 U(\theta^{\star})^{\otimes 2}$ and $\nabla^2 U(\theta^{\star})$, respectively, due to the noise of the stochastic gradients.
It can be shown, see \Cref{subsec:proof-section-bias-variance}, that for $\gamma$ small enough, $\operatorname{H}$ and $\operatorname{G}$ are invertible.

\begin{theorem}\label{thm:bias-variance-tpstar-ula-fp}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}.
  Set $\gamma = \eta/N$ and assume that $\liminf_{N\to+\infty} N^{-1} m  >0$.
  There exists an (explicit) $\eta_0$ independent of $N$ such that for all $\eta\in\ooint{0,\eta_0}$,
  \begin{align}
    \label{eq:covariance-lmc}
    &\int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{LMC}}}(\mathrm{d} \theta) = \operatorname{H}^{-1} (2\operatorname{Id}) + O_{N\to+\infty}(N^{-3/2}) \;, \\
    \label{eq:covariance-fp}
    &\int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{FP}}}(\mathrm{d} \theta)  = \operatorname{G}^{-1} (2\operatorname{Id}) + O_{N\to+\infty}(N^{-3/2}) \;,
  \end{align}
  and
  \begin{align*}

    & \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{LMC}}}(\mathrm{d} \theta) - \theta^{\star} = - \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star}) [\operatorname{H}^{-1} \operatorname{Id}] + O_{N\to+\infty}(N^{-3/2}) \;, \\

    & \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{FP}}}(\mathrm{d} \theta) - \theta^{\star} = - \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star}) [\operatorname{G}^{-1} \operatorname{Id}] + O_{N\to+\infty}(N^{-3/2}) \;.
  \end{align*}
\end{theorem}

\begin{proof}
  The proof is postponed to \Cref{subsec:proof-thm-bias-variance-tpstar}.
\end{proof}

\begin{theorem}\label{thm:bias-variance-tpstar-sgld-sgd}
  Assume \Cref{assum:gradient-Lipschitz}, \Cref{assum:strongly-convex} and \Cref{assum:Ui-convex}. Set $\gamma = \eta/N$ and assume that $\liminf_{N\to+\infty} N^{-1} m  >0$. There exists an (explicit) $\eta_0$ independent of $N$ such that for all $\eta\in\ooint{0,\eta_0}$ and $N\geq 1/\eta$,
  \begin{align}
    \label{eq:covariance-sgld}
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \theta)& = \operatorname{G}^{-1} \defEns{2\operatorname{Id} + (\eta/p)\operatorname{M}} + O_{\eta \to 0}(\eta^{3/2}) \;, \\
    \label{eq:covariance-sgd}
    \int_{\mathbb{R}^d} (\theta - \theta^{\star})^{\otimes 2} \pi_{\scriptscriptstyle{\operatorname{SGD}}}(\mathrm{d} \theta) &= \fraca{\eta}{p} \operatorname{G}^{-1} \operatorname{M}  + O_{\eta \to 0}(\eta^{3/2}) \;,
  \end{align}
  and
  \begin{align*}
    \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{SGLD}}}(\mathrm{d} \theta) - \theta^{\star} &= - (1/2) \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star})[ \operatorname{G}^{-1}\defEns{2\operatorname{Id} + (\eta/p)\operatorname{M}}] + O_{\eta \to 0}(\eta^{3/2}) \;, \\
    \int_{\mathbb{R}^d} \theta \pi_{\scriptscriptstyle{\operatorname{SGD}}}(\mathrm{d} \theta) - \theta^{\star} &= -\fraca{\eta}{2p} \nabla^2 U(\theta^{\star})^{-1} \operatorname{D}^3 U(\theta^{\star})[\operatorname{G}^{-1}\operatorname{M}]  + O_{\eta \to 0}(\eta^{3/2}) \;,
  \end{align*}
  where
  \begin{equation}
    \label{eq:def-opM}
    \operatorname{M} = \sum_{i=1}^{N} \parenthese{\nabla U_i(\theta^{\star}) - \frac{1}{N} \sum_{j=1}^{N} \nabla U_j(\theta^{\star})}^{\otimes 2}\;,
  \end{equation}
  and $\operatorname{G}$ is defined in \eqref{eq:def-operator-G}.
\end{theorem}

\begin{proof}
  The proof is postponed to \Cref{subsec:proof-thm-bias-variance-tpstar}.
\end{proof}

Note that this result implies that the mean and the covariance matrix of $\pi_{\scriptscriptstyle{\operatorname{SGLD}}}$ and $\pi_{\scriptscriptstyle{\operatorname{SGD}}}$ stay lower bounded by a positive constant for any $\eta >0$ as $N \to +\infty$. In \Cref{sec:illustration-m1-m2}, a figure illustrates the results of \Cref{thm:bias-variance-tpstar-ula-fp} and \Cref{thm:bias-variance-tpstar-sgld-sgd} in the asymptotic $N\to+\infty$.