Theoretical guarantees for approximate sampling from smooth and log-concave densities

Introduction

These estimators are rarely available in closed-form. Therefore, optimisation techniques are used for computing the maximum-likelihood estimator while the computation of the Bayes estimator often requires sampling from a density proportional to $e^{-f(\boldsymbol{\theta})}$. In most situations, the exact computation of these two estimators is impossible and one has to resort to approximations provided by iterative algorithms. There is a vast variety of such algorithms for solving both tasks, see for example (Boyd and Vandenberghe, 2004) for optimisation and (Atchadé et al., 2011) for approximate sampling. However, a striking fact is that the convergence properties of optimisation algorithms are much better understood than those of the approximate sampling algorithms. The goal of the present work is to partially fill this gap by establishing easy-to-apply theoretical guarantees for some approximate sampling algorithms.

To be more precise, let us consider the case of a strongly convex function $f$ having a Lipschitz continuous gradient. That is, there exist two positive constants $m$ and $M$ such that

where $\nabla f$ stands for the gradient of $f$ and $\|\cdot\|_{2}$ is the Euclidean norm. There is a simple result characterising the convergence of the well-known gradient descent algorithm under the assumption (1).

This theorem implies that the convergence of the gradient descent is exponential in $k$. More precisely, it results from Eq. (3) that in order to achieve an approximation error upper bounded by $\epsilon>0$ in the Euclidean norm it suffices to perform

evaluations of the gradient of $f$. An important feature of this result is the logarithmic dependence of $k_{\epsilon}$ on $\epsilon$ but also its independence of the dimension $p$. Note also that even though the right-hand side of (4) is a somewhat conservative bound on the number of iterations, all the quantities involved in that expression are easily computable and lead to a simple stopping rule for the iterative algorithm.

The situation for approximate computation of $\boldsymbol{\theta}^{{\rm B}}$ or for approximate sampling from the density proportional to $e^{-f(\boldsymbol{\theta})}$ is much more contrasted. While there exist almost as many algorithms for performing these tasks as for the optimisation, the convergence properties of most of them are studied only empirically and, therefore, provide little theoretically grounded guidance for the choice of different tuning parameters or of the stopping rule. Furthermore, it is not clear how the rate of convergence of these algorithms scales with the growing dimension. While it is intuitively understandable that the problem of sampling from a distribution is more difficult than that of maximising its density, this does not necessarily justifies the huge gap that exists between the precision of theoretical guarantees available for the solutions of these two problems. This gap is even more surprising in light of the numerous similarities between the optimisation and approximate sampling algorithms.

where $h>0$ is a tuning parameter, often referred to as the step-size, and $\boldsymbol{\xi}^{(1)},\ldots,\boldsymbol{\xi}^{(k)},\ldots$ is a sequence of independent centered Gaussian vectors with covariance matrix equal to identity and independent of $\boldsymbol{\vartheta}^{(0)}$. It is well known that under some assumptions on $f$, when $h$ is small and $k$ is large (so that the product $kh$ is large), the distribution of $\boldsymbol{\vartheta}^{(k,h)}$ is close in total variation to the distribution with density proportional to $e^{-f(\boldsymbol{\theta})}$, hereafter referred to as the target distribution. The goal of the present work is to establish a nonasymptotic upper bound, involving only explicit and computable quantities, on the total variation distance between the target distribution and its approximation by the distribution of $\boldsymbol{\vartheta}^{(k,h)}$. We will also analyse a variant of the LMC, termed LMCO, which makes use of the Hessian of $f$.

In order to give the reader a foretaste of the main contributions of the present work, we summarised in Table 1 some guarantees established and described in detail in the next sections. To keep things simple, we translated all the nonasymptotic results into asymptotic ones for large dimension $p$ and small precision level $\epsilon$ (the $O^{*}$ notation ignores the dependence on constant and logarithmic factors). The complexity of one iteration of the LMC indicated in the table corresponds to the computation of the gradient and generation of a Gaussian $p$-vector, whereas the complexity of one iteration of the LMCO is the cost of performing a singular values decomposition on the Hessian matrix of $f$, which is of size $p\times p$.

Background on the Langevin Monte Carlo algorithm

Under assumption (1), for any probability density $\nu$,

The proof of this lemma, postponed to Section 8, is based on the bounds on the spectral gap established in (Chen and Wang, 1997, Remark 4.14), see also (Bakry et al., 2014, Corollary 4.8.2). In simple words, inequality (7) shows that for large values of $t$, the distribution of $\boldsymbol{L}_{t}$ approaches exponentially fast to the target distribution, and the idea behind the LMC is to approximate $\boldsymbol{L}_{t}$ by $\boldsymbol{\vartheta}^{(k,h)}$ for $t=kh$. Note that inequalities of type (7) can be obtained under conditions (such as the curvature-dimension condition, see Bakry et al. (2014, Definition 1.16.1 and Theorem 4.8.4)) weaker than the strong log-concavity required in the present work. However, we decided to restrict ourselves to the strong log-concavity condition since it is easy to check and is commonly used in machine learning and optimisation.

The first and probably the most influential work providing probabilistic analysis of asymptotic properties of the LMC algorithm is (Roberts and Tweedie, 1996). However, one of the recommendations made by the authors of that paper is to avoid using Langevin algorithm as it is defined in (5), or to use it very cautiously, since the ergodicity of the corresponding Markov chain is very sensitive to the choice of the parameter $h$. Even in the cases where the Langevin diffusion is geometrically ergodic, the inappropriate choice of $h$ may result in the transience of the Markov chain $\{\boldsymbol{\vartheta}^{(k,h)}\}$. These findings have very strongly influenced the subsequent studies since all the ensuing research focused essentially on the Metropolis adjusted version of the LMC, known as Metropolis adjusted Langevin algorithm (MALA), and its modifications (Roberts and Rosenthal, 1998; Stramer and Tweedie, 1999a, b; Jarner and Hansen, 2000; Roberts and Stramer, 2002; Pillai et al., 2012; Xifara et al., 2014).

In contrast to this, we show here that under the strong convexity assumption imposed on $f$ (or, equivalently, on $-\log\pi$) coupled with the Lipschitz continuity of the gradient of $f$, one can ensure the non-transience of the Markov chain $\boldsymbol{\vartheta}^{(k,h)}$ by simply choosing $h\leq 1/M$. In fact, the non-explosion of this chain follows from the following proposition the proof of which is very strongly inspired by the one of Theorem 1.

where $\boldsymbol{\theta}^{*}$ stands for the point of (global) minimum of $f$. As a consequence, the sequence $\boldsymbol{\vartheta}^{(k,h)}$ produced by the LMC algorithm is bounded in $L^{2}$ provided that $h\leq 1/M$.

A crucial step in analyzing the long-time behaviour of the LMC algorithm is the assessment of the distance between the distribution of the random variable $\boldsymbol{L}_{Kh}$ and that of $\boldsymbol{\vartheta}^{(K,h)}$. It is intuitively clear that for a fixed $K$ this distance should tend to zero when $h$ tends to zero. However, in order to get informative bounds we need to quantify the rate of this convergence. To this end, we follow the ideas presented in (Dalalyan and Tsybakov, 2009, 2012) which consist in performing the following two steps. First, a continuous-time Markov process $\{\boldsymbol{D}_{t}:t\geq 0\}$ is introduced such that the distribution of the random vectors \big{(}\boldsymbol{\vartheta}^{(0)},\boldsymbol{\vartheta}^{(1,h)},\ldots,\boldsymbol{\vartheta}^{(K,h)}\big{)} and \big{(}\boldsymbol{D}_{0},\boldsymbol{D}_{h},\ldots,\boldsymbol{D}_{Kh}\big{)} coincide. Second, the distance between the distributions of the variables $\boldsymbol{D}_{Kh}$ and $\boldsymbol{L}_{Kh}$ is bounded from above by the distance between the distributions of the continuous-time processes $\{\boldsymbol{D}_{t}:t\in[0,Kh]\}$ and $\{\boldsymbol{L}_{t}:t\in[0,Kh]\}$.

To be more precise, we introduce a diffusion-type continuous-time process $\boldsymbol{D}$ obeying the following stochastic differential equation:

with the (nonanticipative) drift $\boldsymbol{b}_{t}(\boldsymbol{D})=-\sum_{k=0}^{\infty}\nabla f(\boldsymbol{D}_{kh})\mathds{1}_{[kh,(k+1)h[}(t)$. By integrating the last equation on the interval $[kh,(k+1)h]$, we check that the increments of this process satisfy $\boldsymbol{D}_{(k+1)h}-\boldsymbol{D}_{kh}=-h\nabla f(\boldsymbol{D}_{kh})+\sqrt{2h}\boldsymbol{\zeta}^{(k+1)}$, where $\boldsymbol{\zeta}^{(k+1)}=(\boldsymbol{W}_{(k+1)h}-\boldsymbol{W}_{kh})/\sqrt{h}$. Since the Brownian motion is a Gaussian process with independent increments, we conclude that $\{\boldsymbol{\zeta}^{(k)}:k=1,\ldots,K\}$ is a sequence of iid standard Gaussian random vectors. This readily implies the equality of the distributions of the random vectors \big{(}\boldsymbol{\vartheta}^{(0)},\boldsymbol{\vartheta}^{(1,h)},\ldots,\boldsymbol{\vartheta}^{(K,h)}\big{)} and \big{(}\boldsymbol{D}_{0},\boldsymbol{D}_{h},\ldots,\boldsymbol{D}_{Kh}\big{)}.

Note that the specific form of the drift $\boldsymbol{b}$ used in the LMC algorithm has the advantage of meeting the following two conditions. First, $\boldsymbol{b}_{t}(\boldsymbol{L})$ is close to $-\nabla f(\boldsymbol{L}_{t})$, the drift of the Langevin diffusion. Second, it is possible to sample from the distribution $\mathbf{P}^{h}_{\!\!\boldsymbol{D}}(\boldsymbol{x},\>\cdot\>)$, where $h$ is the step of discretisation used in the LMC algorithm. Any nonanticipative drift function satisfying these two conditions may be used for defining a version of the LMC algorithm. Such an example, the LMC algorithm with Ozaki discretisation, is considered in Section 5.

It is worth emphasising that the last inequality remains valid when the initial values of the processes $\boldsymbol{D}$ and $\boldsymbol{L}$ are random but have the same distribution.

Note that the idea of discretising the diffusion process in order to approximately sample from its invariant density is not new. It can be traced back at least to (Lamberton and Pagès, 2002), see also the thesis (Lemaire, 2005) for an overview. The results therein are stated for more general discretisation with variable step-sizes but are of asymptotic nature. This point of view has been adopted and extended to the nonasymptotic case in the recent work (Durmus and Moulines, 2015).

Nonasymptotic bounds on the error of the LMC algorithm

We are now in a position to establish a nonasymptotic bound with explicit constants on the distance between the target distribution $\mathbf{P}_{\pi}$ and the one produced by the LMC algorithm. As explained earlier, the bound is obtained by controlling two types of errors: the error of approximating $\mathbf{P}_{\pi}$ by the distribution of the Langevin diffusion $\boldsymbol{L}_{Kh}$ (6) and the error of approximating the Langevin diffusion by its discretised version $\boldsymbol{D}$ given by (10). The first error is a decreasing function of $T=Kh$: in order to make this error small it is necessary to choose a large $T$. A rather precise quantitative assessment of this error is given by Lemma 1 in the previous section. The second error vanishes when the step-size $h$ goes to zero, provided that $T=Kh$ is fixed. Thus, it is in our interest to choose a small $h$. However, our goal is not only to minimise the error, but also to reduce, as much as possible, the computational cost of the algorithm. For a fixed $T$, if we choose a small value of $h$ then a large number of steps $K$ is necessary for getting close to the target distribution. Therefore, the computational complexity is a decreasing function of $h$. In order to find a value of $h$ leading to a reasonable trade-off between the computational complexity and the approximation error, we need to complement Lemma 1 with a precise bound on the second approximation error. This is done in the following lemma.

Let us set $T=Kh$. Since it simplifies the mathematical formulae and is possible to achieve in practice in view of Theorem 1, we assume in the sequel that the initial value of the LMC algorithm is drawn at random from the Gaussian distribution with mean $\boldsymbol{\theta}^{*}$, a stationary point of $f$, and covariance matrix $M^{-1}\mathbf{I}_{p}$. Then, in view of (12) and the convexity of the Kullback-Leibler divergence, we get (for $\nu=\mathcal{N}_{p}(\boldsymbol{\theta}^{*},M^{-1}\mathbf{I}_{p})$)

for every $K\geq\alpha$ and $h\leq 1/(\alpha M)$. We can now state the main result of this section, the proof of which is postponed to Section 8.

The second term in the right-hand side of (14) tends to infinity when the time horizon $T$ goes to infinity while the step-size $h$ remains fixed. Since the total variation is always bounded by one, the obtained bound is not sharp for large values of $T$. The main reason for this is the fact that we upper bound the total variation distance by the Kullback-Leibler divergence. Improving this argument in order to get a tighter upper bound is a challenging open problem.

We provide here a simple consequence of the last theorem that furnishes easy-to-apply rules for choosing the time horizon $T$ and the step-size $h$.

Let $p\geq 2$, $f$ satisfy (1) and $\epsilon\in(0,1/2)$ be a target precision level. Let the time horizon $T$ and the step-size $h$ be defined by

where $\alpha=(1+MpT\epsilon^{-2})/2$. Then the output of the $K$-step LMC algorithm, with $K=\lceil T/h\rceil$, satisfies \big{\|}\nu\mathbf{P}_{\boldsymbol{\vartheta}}^{K}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\epsilon.

The choice of $T$ and $h$ implies that the two summands in the right-hand side of (14) are bounded by $\epsilon/2$. Furthermore, one easily checks that $\alpha=(1+MpT\epsilon^{-2})/2$ is larger than one and satisfies $h\leq 1/(\alpha M)$. In addition, $K\geq T/h\geq\alpha MT\geq 2\alpha(M/m)\log(1/\epsilon)\geq\alpha\log 4$, which ensures the applicability of Theorem 2. ∎

Let us first remark that the claim of Corollary 1 can be simplified by taking $\alpha=1$. However, for this value of $\alpha$ the factor $(2\alpha-1)/\alpha$ equals one, whereas for the slightly more complicated choice recommended by Corollary 1, this factor is close to two. In practice, increasing $h$ by a factor $2$ results in halving the running time, which represents a non-negligible gain.

Besides providing concrete and easily applicable guidance for choosing the step of discretisation and the stopping rule for the LMC algorithm to achieve a prescribed error rate, the last corollary tells us that in order to get an error smaller than $\epsilon$, it is enough to perform K=O(T^{2}p/\epsilon^{2})=O\big{(}\epsilon^{-2}(p^{3}+p\log^{2}(1/\epsilon))\big{)} evaluations of the gradient of $f$. To the best of our knowledge, this is the first result that establishes polynomial in $p$ guarantees for sampling from a log-concave density using the LMC algorithm. We discuss the relation of this and subsequent results to earlier work in Section 7.

Possible extensions

In this section, we state some extensions of the previous results that do not require any major change in the proofs, but might lead to improved computational complexity or be valid under relaxed assumptions in some particular cases.

The choice of the distribution $\nu$ of the initial value $\boldsymbol{\theta}^{(0)}$ has a significant impact on the convergence of the LMC algorithm. If $\nu$ is close to $\pi$, smaller number of iterations might be enough for making the TV-error smaller than $\epsilon$. The goal of this section is to present quantitative bounds characterising the influence of $\nu$ on the convergence and, as a consequence, on the computational complexity of the LMC algorithm.

The first observation that can be readily deduced from (12) is that for any $h\leq 1/(2M)$,

Combining this bound with (38), Lemma 1 and (40) we get

Elaborating on this inequality, we get the following result.

satisfies, for $K=[T/h]\geq 2$, the inequality \big{\|}\nu\mathbf{P}_{\boldsymbol{\vartheta}}^{K}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\epsilon.

The proof of this proposition is immediate and, therefore, is left to the reader. What we infer from this result is that the choice of the initial distribution $\nu$ has a strong impact on the convergence of the LMC algorithm. For instance, if for some specific $\pi$ we are able to sample from a density $\nu$ satisfying, for some $\varrho>0$, the relation $\chi^{2}(\nu\|\pi)=O(p^{\varrho})$ as $p\to\infty$, then the time horizon $T$ for approximating the target density $\pi$ within $\epsilon$ is $O(\log(p\vee\epsilon^{-1}))$ and the step-size satisfies $h^{-1}=O(\epsilon^{-2}p\log(p\vee\epsilon^{-1}))$. Thus, in such a situation, one needs to perform $[T/h]=O(\epsilon^{-2}p\log^{2}(p\vee\epsilon^{-1}))$ evaluations of the gradient of $f$ to get a sampling density within a distance of $\epsilon$ of the target, which is substantially smaller than O\big{(}\epsilon^{-2}(p^{3}+p\log^{2}(1/\epsilon))\big{)} obtained in the previous section in the general case.

2 Preconditioning

i.e., the approximation error of the LMC algorithm with a preconditioner $\mathbf{A}$ is characterised by Corollary 1. This means that if the function $g$ satisfies condition (1) with constants $(m_{\mathbf{A}},M_{\mathbf{A}})$, then the number of steps $K$ after which the preconditioned LMC algorithm has an error bounded by $\epsilon$ is given by K=(M_{\mathbf{A}}/m_{\mathbf{A}})^{2}p\epsilon^{-2}\big{(}2\log(1/\epsilon)+(p/2)\log(M_{\mathbf{A}}/m_{\mathbf{A}})\big{)}^{2}. Hence, the preconditioner $\mathbf{A}$ yielding the best guaranteed computational complexity for the LMC algorithm is the matrix $\mathbf{A}$ minimising the ratio $M_{\mathbf{A}}/m_{\mathbf{A}}$.

The impact of preconditioning can be measured, for instance, in the case of multidimensional logistic regression considered in Section 6 below. In this case, the ratio $M_{\mathbf{A}}/m_{\mathbf{A}}$ is up to some constant factor equal to the condition number of the matrix $\mathbf{A}\Sigma_{\mathbf{X}}\mathbf{A}$, where $\Sigma_{\mathbf{X}}$ is the Gram matrix of the covariates.

3 Nonstrongly log-concave densities

Theoretical guarantees developed in previous sections assume that the logarithm of the target density is strongly concave, cf. assumption (1). However, they can also be used for approximate sampling from a density which is log-concave but not necessarily strongly log-concave; we call these densities nonstrongly log-concave. The idea is then to approximate the target density by a strongly log-concave one and to apply the LMC algorithm to the latter instead of the former one.

As a consequence, \big{\|}\mathbf{P}_{\bar{\pi}}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\frac{1}{2}\|\bar{f}-f\|_{L^{2}(\pi)}.

Using the formula for the Kullback-Leibler divergence, we get

Applying successively the inequalities $\log u\leq u-1$ and $e^{-u}\leq 1-u+\frac{1}{2}u^{2}$ for every $u\geq 0$, we upper bound the second term in the right-hand side of (20) as follows:

Combining this inequality with (20), we get the first claim. The last claim of the lemma follows from the Pinsker inequality. ∎

For $\bar{f}$ given by (18), we get \big{\|}\mathbf{P}_{\bar{\pi}}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\frac{\gamma}{4}\big{(}\int_{B^{c}}(\|\boldsymbol{x}-\boldsymbol{x}_{0}\|_{2}-R)^{4}\,\pi(\boldsymbol{x})\,d\boldsymbol{x}\big{)}^{1/2}. Choosing the parameter $\gamma$ sufficiently small and the parameter $R$ sufficiently large to ensure that $\|\mathbf{P}_{\bar{\pi}}-\mathbf{P}_{\pi}\|_{\rm TV}\leq\epsilon/2$ and assuming that $\pi$ has bounded fourth-order moment, we derive from this inequality and Corollary 1 the following convergence result for the approximate LMC algorithm.

Let $f$ be a twice differentiable function satisfying $m_{R}\mathbf{I}_{p}\preceq\nabla^{2}f(\boldsymbol{x})\preceq M\mathbf{I}_{p}$ for every $\boldsymbol{x}\in B_{R}(\boldsymbol{x}_{0})$ and for every $R\in[0,+\infty]$. Let $\epsilon\in(0,1/2)$ be a target precision level. Assume that for some known value $\mu_{R}$ we have $\int_{B_{R}(\boldsymbol{x}_{0})^{c}}(\|\boldsymbol{x}-\boldsymbol{x}_{0}\|_{2}-R)^{4}\pi(\boldsymbol{x})\,d\boldsymbol{x}\leq p^{2}\mu_{R}^{2}$ and define $\bar{m}=m_{2R}\wedge(m_{\infty}+0.5\gamma)$, $\bar{M}=M+\gamma$ for some $\gamma\leq 2\epsilon/(p\mu_{R})$. Set the time horizon $T$ and the step-size $h$ as follows:

Then the output of the $K$-step LMC algorithm (5) applied to the approximation $\bar{f}$ provided by (18), with $K=\lceil T/h\rceil$, satisfies \big{\|}\nu\mathbf{P}_{\boldsymbol{\vartheta}}^{K}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\epsilon.

Let us comment this result in the case $R=0$ which concerns nonstrongly log-concave densities. Then the previous result implies that $K=O(p^{5}\epsilon^{-4}\log^{2}(p\vee\epsilon^{-1}))$. Clearly, the dependence of $K$ both on the dimension $p$ and on the acceptable error level $\epsilon$ gets substantially deteriorated as compared to the strongly log-concave case. Some improvements are possible in specific cases. First, we can improve the dependence of $K$ on $p$ if we are able to simulate from a distribution $\nu$ that is not too far from $\bar{\pi}$ in the sense of $\chi^{2}$ divergence. More precisely, repeating the arguments of Section 4.1 we get the following result: if the initial distribution of the LMC algorithm satisfies \chi^{2}(\nu\|\bar{\pi})=O\big{(}(p/\gamma)^{\varrho}\big{)} for some $\varrho>0$ then one needs at most K=O\big{(}p^{3}\epsilon^{-4}\log^{2}(p/\epsilon)\big{)} steps of the LMC algorithm for getting an error bounded by $\epsilon$. Second, in some cases the dependence of $K$ on $p$ can be further improved by using a preconditioner and/or by replacing the penalty $\|\boldsymbol{x}\|_{2}^{2}$ in (18) by $\|\mathbf{M}\boldsymbol{x}\|_{2}^{2}$, where $\mathbf{M}$ is a properly chosen $p\times p$ matrix.

This being said, our intuition is that Corollary 2 is more helpful in the case of convex functions $f$ that are strongly convex in a neighbourhood of their minimum point $\boldsymbol{\theta}^{*}$. In such a situation, our recommendation is to set $\boldsymbol{x}_{0}=\boldsymbol{\theta}^{*}$ and to choose $R$ by maximising the quantity $\bar{m}=m_{2R}\wedge(m_{\infty}+\epsilon/(p\mu_{R}))$. We showcase this approach in Section 6 on the example of logistic regression.

Note that the convergence of the MCMC methods for sampling from log-concave densities was also studied in (Brooks, 1998), where a strategy for defining the stopping rule is proposed. However, as the computational complexity of that strategy increases exponentially fast in the dimension $p$, its scope of applicability is limited.

Ozaki discretisation and guarantees for smooth Hessian matrices

For convex log-densities $f$ which are not only continuously differentiable but also have a smooth Hessian matrix $\nabla^{2}f$, it is possible to take advantage of the Ozaki discretisation (Ozaki, 1992) of the Langevin diffusion which is more accurate than the Euler discretisation analysed in the foregoing sections. It consists in considering the diffusion process $\boldsymbol{D}^{O}$ defined by (10) with the drift function

where, as previously, $h$ is the step-size and $K$ is the number of iterations to attain the desired time horizon $T=Kh$. This expression leads to a diffusion process having linear drift function on each interval $[kh,(k+1)h[$. Such a diffusion admits a closed-form formula. The resulting MCMC algorithm (Stramer and Tweedie, 1999b), hereafter referred to as LMCO algorithm (for Langevin Monte Carlo with Ozaki discretisation), is defined by an initial value $\bar{\boldsymbol{\vartheta}}^{(0)}$ and the following update rule. For every $k\geq 0$, we set $\mathbf{H}_{k}=\nabla^{2}f(\bar{\boldsymbol{\vartheta}}^{(k,h)})$, which is an invertible $p\times p$ matrix since $f$ is strongly convex, and define

The proof of this theorem is deferred to Section 8. Let us state now a direct consequence of the last theorem, which provides sufficient conditions on the number of steps for the LMCO algorithm to achieve a prescribed precision level $\epsilon$. The proof of the corollary is trivial and, therefore, is omitted.

Let $f$ satisfy (1) with a Hessian that is Lipschitz-continuous with constant $L_{f}$. For every $\epsilon\in(0,1/2)$, if the time horizon $T$ and the step-size $h$ are chosen so that

then the distribution of the outcome of the LMCO algorithm with $K=[T/h]$ steps fulfils \big{\|}\nu\mathbf{P}_{\bar{\boldsymbol{\vartheta}}}^{K}-\mathbf{P}_{\pi}\big{\|}_{\rm TV}\leq\epsilon.

This corollary provides simple recommendation for the choice of the parameters $h$ and $T$ in the LMCO algorithm. It also ensures that for the recommended choice of the parameters, it is sufficient to perform K=O\big{(}(p+\log(1/\epsilon))^{3/2}p\epsilon^{-1}\big{)} number of steps of the LMCO algorithm in order to reach the desired precision level $\epsilon$. This number is much smaller than that provided earlier by Corollary 1, which was of order O\big{(}(p+\log(1/\epsilon))^{2}p\epsilon^{-2}\big{)}. However, one should pay attention to the fact that each iteration of the LMCO requires computing the exponential of the Hessian of $f$ at the current state and, therefore, the computational complexity of each iteration is usually much larger for the LMCO as compared to the LMC ($O(p^{3})$ versus $O(p)$). This implies that the LMCO would most likely be preferable to the LMC only in situations where $p$ is not too large, and the required precision level $\epsilon$ is very small. For instance, the arguments of this paragraph advocate for using the LMCO instead of the LMC when $\epsilon=o(p^{-3/2})$.

This being said, it is worth noting that for some functions $f$ the cost of performing a singular values decomposition on the Hessian of $f$, which is the typical way of computing the matrix exponential, might be much smaller than the aforementioned worst-case complexity $O(p^{3})$. This is, in particular, the case for the first example considered in the next section. One can also approximate the matrix exponentials by matrix polynomials. For second-order polynomials, this amounts to replacing the updates (24) by

Establishing guarantees for such a modified LMCO is out of scope of the present work. We will limit ourselves to an empirical assessment of the quality of this approximation on the example of logistic regression considered in Section 6.

To close this section, let us remark that in the case a warm start is available, the number of iterations for the LMCO algorithm to reach the precision $\epsilon$ may be reduced to $O^{*}(p\epsilon^{-1})$. Indeed, if the $\chi^{2}$ divergence between the initial distribution and the target is bounded by a quantity independent of $p$, or increasing not faster than a polynomial in $p$, then the time horizon can be chosen as $O^{*}(1)$ and the choice of $h$ provided by Corollary 3 leads to a number of iterations $K$ satisfying $K=O^{*}(p\epsilon^{-1})$.

Numerical experiments

To illustrate the results established in the previous sections, we carried out some experiments on synthetic data. The experiments were conducted on a HP Elitebook PC with the following configuration: Intel (R) Core (TM) i7-3687U with 2.6 GHz CPU and 16 GB of RAM. The code, written in Matlab, does not use parallelisation. We considered two examples; both satisfy all the assumptions required in previous sections. This implies that Corollaries 1 and 3 apply and guarantee that the choices of $h$ and $T$ suggested by these corollaries allow us to generate random vectors having a distribution which is within a prescribed distance $\epsilon$, in total variation, of the target distribution.

The goal of this first experiment is merely to show on a simple example the validity of our theoretical findings. That is, we check below that the LMC and the LMCO algorithms with the values of time horizon $T$ and step-size $h$ recommended by Corollaries 1 and 3 produce samples distributed approximately as the target distribution within a reasonable running time. To this end, we consider the simple task of sampling from the density $\pi$ defined by

Using the fact that 0\leq 4e^{2\boldsymbol{x}^{\top}\mathbf{a}}\big{(}1+e^{2\boldsymbol{x}^{\top}\mathbf{a}}\big{)}^{-2}\leq 1, we infer that for $\|\mathbf{a}\|_{2}<1$, the function $f$ is strongly convex and satisfies (1) with $m=1-\|\mathbf{a}\|_{2}^{2}$ and $M=1$. Furthermore, the Hessian matrix is Lipschitz continuous with the constant $L_{f}=\frac{1}{2}\|\mathbf{a}\|_{2}^{3}$. Hence, both algorithms explored in the previous sections, LMC and LMCO, can be used for sampling from the density $\pi$ defined by (26). Note also that one can sample directly from $\pi$ by drawing independently at random a Bernoulli$(1/2)$ random variable $Y$ and a standard Gaussian vector $\boldsymbol{Z}\sim\mathcal{N}(0,\mathbf{I}_{p})$ and by computing $\boldsymbol{X}=Y\cdot(\boldsymbol{Z}-\mathbf{a})+(1-Y)\cdot(\boldsymbol{Z}+\mathbf{a})$. The density of the random vector $\boldsymbol{X}$ defined in such a way coincides with $\pi$. One can check that the unique minimum of $f$ is achieved at $\boldsymbol{\theta}^{*}=c^{*}\cdot\mathbf{a}$, where $c^{*}$ is the unique solution of the equation $c=1-2(1+e^{2c\|\mathbf{a}\|_{2}^{2}})^{-1}$. Choosing $\mathbf{a}$ so that $\|\mathbf{a}\|_{2}^{2}=1/2$, we get $\boldsymbol{\theta}^{*}=0$.

To illustrate the dependence on the dimension $p$ of the computational complexity of the proposed sampling strategies, we report in Table 2 the number of iterations and the overall running times for generating $N=10^{3}$ independent samples by the LMC and the LMCO for the target specified by (26), when the dimension $p$ varies in $\{4,8,12,16,20,30,40,60\}$. One may observe that the computational time is much smaller for the LMCO than for the LMC algorithm, which is mainly explained by the fact that the singular vectors of the Hessian of the function $f$, in the example under consideration, do not depend on the value $\boldsymbol{x}$ at which the Hessian is computed.

This example confirms our theoretical findings in that it shows that (a) the samples drawn from the LMC and the LMCO algorithms with the parameters $T$ and $h$ suggested by theoretical considerations have distributions that are very close to the target distribution and that (b) the running-times for these algorithms remain reasonable even for moderately large values of dimension $p$.

Example 2: Binary logistic regression

where $\boldsymbol{Y}=(Y_{1},\ldots,Y_{n})^{\top}\in\{0,1\}^{n}$ and $\mathbf{X}$ is the $n\times p$ matrix having the feature $\boldsymbol{X}_{i}$ as $i^{\rm th}$ row. The first two terms in the exponential correspond to the log-likelihood of the logistic model, whereas the last term comes from the log-density of the prior and can be seen as a penalty term. The parameter $\lambda>0$ is usually specified by the practitioner. Many authors have studied this model from a Bayesian perspective, see for instance (Holmes and Held, 2006; Roy, 2012), and it seems that there is no compelling alternative to the MCMC algorithms for computing the Bayesian estimators in this model. Furthermore, even for the MCMC approach, although geometric ergodicity under some strong assumptions is established, there is no theoretically justified rule for assessing the convergence and, especially, ensuring that the convergence is achieved in polynomial time. Such guarantees are provided by our results, when either the LMC or the LMCO is used.

we get the setting described in the Introduction. It is useful here to apply the preconditioning technique of Section 4.2 with the preconditioner $\mathbf{A}=\boldsymbol{\Sigma}_{\mathbf{X}}^{-1/2}$. Thus, the LMC and the LMCO can be used with the function $f$ replaced by $g(\boldsymbol{\theta})=f(\mathbf{A}\boldsymbol{\theta})$. One checks that $g$ and $f$ are infinitely differentiable and

For the function $g$, since $\nabla^{2}g(\boldsymbol{\theta})=\mathbf{A}\nabla^{2}f(\mathbf{A}\boldsymbol{\theta})\mathbf{A}$, we can infer from these relations that (1) holds with $m_{\mathbf{A}}=\lambda$ and $M_{\mathbf{A}}=\lambda+0.25n$. Note here that if we do not use any preconditioner, the constants $m$ and $M$ would be given by $m=\lambda\cdot\nu_{\min}(\boldsymbol{\Sigma}_{\mathbf{X}})$ and $M=(\lambda+0.25n)\cdot\nu_{\max}(\boldsymbol{\Sigma}_{\mathbf{X}})$, where $\nu_{\min}(\boldsymbol{\Sigma})$ and $\nu_{\max}(\boldsymbol{\Sigma})$ are respectively the smallest and the largest eigenvalues of $\boldsymbol{\Sigma}$. This implies that the ratio $\nu_{\max}(\boldsymbol{\Sigma}_{\mathbf{X}})/\nu_{\min}(\boldsymbol{\Sigma}_{\mathbf{X}})$ quantifies the gain of efficiency obtained by preconditioning. This ratio might be large especially when $p$ is large and the covariates are strongly correlated.

Furthermore, $\nabla^{2}g$ is Lipschitz with a constant $L_{g}$ provided by the following formula (the proof of which is postponed to Section 8):

In our second experiment, for a set of values of $p$ and $n$, we randomly drew $n$ iid samples $(\boldsymbol{X}_{i},Y_{i})$ according to the following data generating device. The features $\boldsymbol{X}_{i}$ were drawn from a Rademacher distribution (i.e., each coordinate takes the values $\pm 1$ with probability $1/2$), and then renormalised to have an Euclidean norm equal to one. Each label $Y_{i}$, given $\boldsymbol{X}_{i}=\boldsymbol{x}$, was drawn from a Bernoulli distribution with parameter $r({\boldsymbol{\theta}^{\rm true}},\boldsymbol{x})$. The true vector $\boldsymbol{\theta}^{\rm true}$ was set to ${\bf 1}_{p}=(1,1,\ldots,1)^{\top}$. For each value of $p$ and $n$, we generated $100$ samples $(\mathbf{X},\boldsymbol{Y})$. For each sample, we computed the MLE using the gradient descent as described in Theorem 1 with a precision level $\epsilon=10^{-6}$. Following the recommendation of (Hanson et al., 2014), the parameter $\lambda$ was set to $3p/\pi^{2}$. We carried out two sub-experiments with well specified distinct purposes: to empirically assess the gain obtained by applying the trick of strong-convexification described in Subsection 4.3 and to evaluate the loss of accuracy caused by applying to the LMCO algorithm the second-order approximation (25).

In the first sub-experiment, we applied the strategy outlined in Subsection 4.3 for various values of $n,p$ and $\epsilon$. To this end, we exploited the following formulae

where $\Gamma(p;x)=\int_{x}^{\infty}t^{p-1}e^{-t}\,dt$ is the upper incomplete gamma function and $C_{4}^{j}$ stands for the binomial coefficient. The proof of the fact that the quantities $m_{R}$ and $\mu_{R}$ defined by these formulae satisfy all the assumptions of Subsection 4.3 is provided in the supplementary material. In this experiment, we used two values of $\epsilon$ ($0.1$ and $0.01$), three values of dimension $p$ ($2$, $5$ and $20$), and five values for the sample size $n$ (500, 1000, 2000, 4000 and 8000). We reported in Table 3 the number of iterates using the LMC algorithm ($K$) and the average number of iterates of the modified LMC algorithm as described in Subsection 4.3 ($K^{\prime}$). Note that in the case of modified LMC algorithm, the number of iterates depends on the original data $(\mathbf{X},\boldsymbol{Y})$. Therefore, the numbers $K^{\prime}$ reported in Table 3 are those obtained by averaging over 100 independent trials.

The results of Table 3 show clearly the advantage of using the strong-convexification trick. For instance, when $\epsilon=0.1$, $p=5$ and $n=1000$, the gain is very impressive since the number of iterations is reduced from nearly $7.5\times 10^{5}$ to $2.2\times 10^{3}$. This represents a reduction by a factor close to 340. The gain is less significant in the case when the ratio $p/n$ is larger. Our explanation of this phenomenon is that for a small ratio $p/n$, the posterior density has a very strong peak at its mode. Therefore, even for a relatively large radius $R$ the condition number $M/m_{R}$ is not too large. Thus, small $p/n$ is the typical situation in which the strong-convexification trick is likely to lead to considerable savings in running-time.

In the second sub-experiment, we aimed at verifying the validity of the second-order approximation of the LMCO algorithm, hereafter referred to as LMCO’, obtained by applying the update rule (25). To this end, for $\epsilon=0.1$, $p\in\{2,5,10\}$ and for $n\in\{200,300,400,500\}$, we generated $N_{\text{MC}}=100$ Monte-Carlo samples using the LMC algorithm and the LMCO’ algorithm. To check the closeness of the distributions of these two $p$-dimensional samples, we compared several aspects of them. More precisely, we compared their marginal means, marginal medians and marginal quartiles. Mathematically speaking, for each data-set $\mathcal{D}_{\text{data}}=(\mathbf{X},\boldsymbol{Y})$, we generated $N_{\text{MC}}$ samples $\mathcal{D}_{\text{MC}}=\{\boldsymbol{\theta}^{1},\ldots,\boldsymbol{\theta}^{N_{\text{MC}}}\}$ and $\bar{\mathcal{D}}_{\text{MC}}=\{\bar{\boldsymbol{\theta}}^{1},\ldots,\bar{\boldsymbol{\theta}}^{N_{\text{MC}}}\}$ using the LMC and the LMCO’, respectively. We then computed the normalised distance between their marginal means: $d_{\text{mean}}=\frac{1}{p}\|\mathop{\rm mean}(\mathcal{D}_{\text{MC}})-\mathop{\rm mean}(\bar{\mathcal{D}}_{\text{MC}})\|_{1}$. We also computed the quantities $d_{\rm median}$, $d_{Q_{1}}$ and $d_{Q_{3}}$, which are defined analogously by replacing the mean by the coordinate-wise median, first quartile and third quartile, respectively. The idea for considering these quantities is that, for large $N_{\rm MC}$ and small $\epsilon$, all the aforementioned distances should be close to zero.

We opted for the boxplot representation of 100 values of each of these distances obtained over 100 independent replications of the data-set $\mathcal{D}_{\text{data}}$. These boxplots are drawn in Fig. 2. They show that the distances are small—at most of the order of $10^{-1}$—which may be considered as an argument in favor of the modification proposed in (25). Indeed, with $\epsilon=0.1$ and $N_{\rm MC}=100$, we could not expect to have an error of smaller order. This is very promising since this modified LMCO algorithm has a significantly smaller computational complexity than the original LMCO: each iteration has a worst-case accuracy $O(p^{2})$ instead of $O(p^{3})$, thanks to the fact that matrix exponentials as well as the inversion of the Hessian are replaced by the computation of the Hessian and its product with vectors.

Summary and conclusion

We have established easy-to-implement, nonasymptotic theoretical guarantees for approximate sampling from log-concave and strongly log-concave probability densities. To this end, we have analysed the Langevin Monte Carlo (LMC) algorithm and its Ozaki discretised version LMCO. These algorithms can be regarded as the natural counterparts—when the task of optimisation is replaced by the task of sampling—of the gradient descent algorithm, widely studied in convex optimisation. Despite its broad applicability in the framework of Bayesian statistics and beyond, to the best of our knowledge, there were no theoretical result in the literature proving that the computational complexity of the aforementioned algorithms scales at most polynomially in dimension and in $\epsilon^{-1}$, the inverse of the desired precision level. The results proved in the present work fill this gap by showing that in order to achieve a precision (in total variation) bounded from above by $\epsilon$, the LMC needs no more than O\big{(}\epsilon^{-2}(p^{3}+p\log(\epsilon^{-1}))\big{)} evaluations of the gradient when the target density is strongly log-concave and O\big{(}\epsilon^{-4}p^{5}\log^{2}(p\vee\epsilon^{-1}))\big{)} evaluations of the gradient when the target density is nonstrongly log-concave. Further improvement of the rates can be achieved if a “warm start” is available. More precisely, if there is an efficiently samplable distribution $\nu$ such that the chi-squared divergence between $\nu$ and the target scales polynomially in $p$, then the LMC with an initial value drawn from $\nu$ needs no more than O\big{(}\epsilon^{-2}p\log^{2}(p\vee\epsilon^{-1}))\big{)} evaluations of the gradient when the target density is strongly log-concave and O\big{(}\epsilon^{-4}p^{3}\log^{2}(p\vee\epsilon^{-1}))\big{)} gradient evaluations when the target density is nonstrongly log-concave. An important advantage of our results is that all the bounds come with explicit numerical constants of reasonable magnitude.

The search for tractable theoretical guarantees for MCMC algorithms is an active topic of research not only in probability and statistics but also in theoretical computer science and in machine learning. To the best of our knowledge, first computable bounds on the constants involved in the geometric convergence of Markov chains were derived in (Meyn and Tweedie, 1994), see also subsequent work (Rosenthal, 2002; Douc et al., 2004) and the survey paper (Roberts and Rosenthal, 2004). However, because of the broad generality of the considered Markov processesThe authors do not confine their study to the log-concave densities., their results are difficult to implement for getting tight bounds on the constants in the context of high dimensionality. In particular, we did not succeed in deriving from their results convergence rates for the LMC algorithm (neither for its Metropolis-Hastings-adjusted version, MALA) that are polynomial in the dimension $p$ and hold for every strongly log-concave target density. Note also that some nonasymptotic convergence results for the MALA were obtained by Bou-Rabee and Hairer (2013), where strongly log-concave four times continuously differentiable functions $f$ were considered. Unfortunately, the constants involved in their bounds are not explicit and cannot be used for our purposes.

The problem of sampling from log-concave distributions is not new. It has been considered in early references (Frieze et al., 1994) and (Frieze and Kannan, 1999). An important progress in this topic was made in a series of papers by Lovázs and Vempala (see, in particular, Lovász and Vempala (2006b, a) for the sharpest results), which are perhaps the closest to our work. They investigated the problem of sampling from a log-concave density with a compact support and derived nonasymptotic bounds on the number of steps that are sufficient for approximating the target density; the best bounds are obtained for the hit-and-run algorithm. The analysis they carried out is very different from the one presented in the present work and the constants in their results are prohibitively large (for instance, $10^{31}$ in (Lovász and Vempala, 2006b, Corollary 1.2)), which makes the established guarantees of little interest for practice. On the positive side, one of the most remarkable features of the results proved in (Lovász and Vempala, 2006b, a) is that the number of steps required to achieve the level $\epsilon$ scales polylogarithmically in $1/\epsilon$. This is of course much better than the dependence on $\epsilon$ in our bounds. However, the logarithm of $1/\epsilon$ in their result is raised to power $5$, which for most interesting values of $\epsilon$ behaves itself as a linear function of $1/\epsilon$. On the down side, the dependence on the dimension in the results of Lovász and Vempala (2006b, a), when no warm start is available, scales as $p^{4}$, which is worse than $p^{3}$ inferred from our analysis. A difference worth being stressed between our framework and that of Lovász and Vempala (2006b, a) is that the LMC algorithm we have analysed here is based on the evaluations of the gradient of $f$, whereas the algorithms studied in (Lovász and Vempala, 2006b, a) need to sample from the restriction of $\pi_{f}$ on the lines. On a related note, building on the results by Lovàzs and Vempala, Belloni and Chernozhukov (2009) provided polynomial guarantees for sampling from a distribution which converges asymptotically to a Gaussian one.

After the submission of the present paper, the manuscript (Durmus and Moulines, 2015) has been posted on arXiv, which refines our results in various directions. In particular, the authors of that manuscript manage to assess more accurately the impact of the initial distribution on the final precision of the LMC algorithm and investigate an Euler scheme with nonconstant step-size. Roughly speaking, they prove that the rate we obtained in the case of a warm start is valid for any starting point which is not too far away from the mode of the density. On a related note, we focus in the present work only on the total variation distance between some MCMC algorithms and the target distribution, whereas in many applications one may be only interested in approximating integrals with respect to the target distribution. Clearly, guarantees on the total variation distance imply guarantees on the approximations of integrals, at least when the integrands are bounded functions. However, since the problem of approximating integrals is, in some sense, easier than sampling from a distribution, one could hope to get tighter bounds for the former problem. This and related questions are thoroughly investigated in (Durmus and Moulines, 2015).

Although the main contribution of the present work is of theoretical nature, we can also draw some conclusions which might be of interest for practitioners. First of all, our results show that the heuristic choice of the stopping rule for the MCMC algorithms is not the only possible option: it is also possible to have theoretically grounded guidelines for choosing the stopping time. The resulting algorithm will be of polynomial complexity both in dimension and in the precision level. Second, the results reported in this work show that there is no need to apply Metropolis-Hastings correction to the Langevin algorithm and its various variants in order to ensure their convergence. Third, when the dimension is not very high and a high level of precision is required (i.e., when $p^{3/2}\epsilon$ is small), the LMCO algorithm is preferable to the LMC algorithm, and the modified LMCO using the update rule of Eq. (25) is even better. Note, however, that this last claim was checked empirically but comes without any theoretical justification.

Finally, we would like to mention that, in recent years, several studies making the connection between convex optimisation and MCMC algorithms were carried out. They mainly focused on proposing new algorithms of approximate sampling (Girolami and Calderhead, 2011; Schreck et al., 2013; Pereyra, 2014) inspired by the ideas coming from convex optimisation. We hope that the present work will stimulate a more extensive investigation of the relationship between approximate sampling and optimisation, especially in the aim of establishing user friendly theoretical guarantees for the MCMC algorithms.

Postponed proofs and some technical results

After a suitable rearrangement of the terms we get the claim of Lemma 5. ∎

2 Proofs of results concerning the LMC

Instead of proving Proposition 1, we prove below the following stronger result.

Throughout this proof, we use the shorthand notation $f^{(k)}=f(\boldsymbol{\vartheta}^{(k,h)})$ and $\nabla f^{(k)}=\nabla f(\boldsymbol{\vartheta}^{(k,h)})$. In view of the relation (5) and the Taylor expansion, we have

Taking the expectations of both sides, we get

Applying this inequality to $\boldsymbol{x}=\boldsymbol{\vartheta}^{(k,h)}$ and combining it with (8.2), whenever $h<2/M$ we get

Let us set $\gamma=mh(2-Mh)\in(0,1)$ for any $h\in(0,2/M)$. Subtracting $f^{*}$ from the both sides of (33) we arrive at

Inequality (30) follows by replacing $\gamma$ by $mh(2-Mh)$. To prove (31), it suffices to combine (30) with the first inequality in (1), Lemma 4 and the inequality $(1-mh)^{k}\leq e^{-mhk}$. ∎

Let $h\leq 1/\alpha M$ with $\alpha\geq 1$ and $K\geq 1$ be an integer. Under the conditions of Proposition 1, it holds

Using inequality (8.2) and the fact that $2-Mh\geq(2\alpha-1)/\alpha$, we get

Summing up these inequalities for $k=0,\ldots,K-1$ and using the obvious bound $f^{(K)}\geq f^{*}$, we get

To complete the proof, it suffices to remark that in view of Lemma 4, it holds 2\mathbf{E}\big{[}f^{(0)}-f^{*}\big{]}\leq M\mathbf{E}\big{[}\|\boldsymbol{\vartheta}^{(0)}-\boldsymbol{\theta}^{*}\|_{2}^{2}\big{]}. ∎

Using the Cauchy-Schwarz inequality, we get

For every fixed Borel set $A$, if we set $\varphi(\boldsymbol{x})=\mathds{1}_{A}(\boldsymbol{x})-\pi(A)$ and use (36), we obtain that

Since $\nabla f$ is Lipschitz continuous with Lipschitz constant $M$, we have

Applying Corollary 4, the desired inequality follows. ∎

In view of the triangle inequality, we have

The first term in the right-hand side is what we call first type error. The source of this error is the finiteness of time, since it would be equal to zero if we could choose $T=Kh=+\infty$. The second term in the right-hand side of (38) is the second type error, which is caused by the practical impossibility to take the step-size $h$ equal to zero. These two errors can be evaluated as follows.

For the first type error, apply Lemma 1 to get $\|\nu\mathbf{P}_{\boldsymbol{L}}^{T}-\mathbf{P}_{\pi}\|_{\rm TV}\leq\frac{1}{2}\chi^{2}(\nu\|\pi)^{1/2}e^{-Tm/2}$. Since $\nu$ is a Gaussian distribution, the expectation in the above formula is not difficult to evaluate. The corresponding result, provided by Lemma 5, yields

To evaluate the second type error, we use the Pinsker inequality:

Combining this inequality with (3), we get the desired result. ∎

3 Proofs of results concerning the LMCO

Using the same arguments as those of the proof of Theorem 2. This leads to the inequality

Since on each interval $[kh,(k+1)h[$ the function $t\mapsto b_{t}$ is linear, for every $t\in[kh,(k+1)h[$, we get \big{\|}\nabla f(\boldsymbol{D}_{t}^{O})+b_{t}(\boldsymbol{D}^{O})\big{\|}_{2}^{2}=\big{\|}\nabla f(\boldsymbol{D}_{t}^{O})-\nabla f(\boldsymbol{D}_{kh}^{O})-\nabla^{2}f(\boldsymbol{D}_{kh}^{O})(\boldsymbol{D}^{O}_{t}-\boldsymbol{D}_{kh}^{O})\big{\|}_{2}^{2}. Using the mean-value theorem and the Lipschitz continuity of the Hessian of $f$, we derive from the above relation that

for every $t\in[kh,(k+1)h[$. Note now that equation (24) provides the conditional distribution of $\boldsymbol{D}_{(k+1)h}^{O}$ given $\boldsymbol{D}_{kh}^{O}$. An analogous formula holds for the conditional distribution of $\boldsymbol{D}_{t}^{O}-\boldsymbol{D}_{kh}^{O}$ given $\boldsymbol{D}_{kh}^{O}$, which is multivariate Gaussian with mean \big{(}\mathbf{I}_{p}-e^{-(t-kh)\mathbf{H}_{k}}\big{)}\mathbf{H}_{k}^{-1}\nabla f\big{(}\boldsymbol{D}_{kh}^{O}\big{)} and covariance matrix \boldsymbol{\Sigma}_{k}=\big{(}\mathbf{I}_{p}-e^{-2(t-hk)\mathbf{H}_{k}}\big{)}\mathbf{H}_{k}^{-1}, where $\mathbf{H}_{k}=\nabla^{2}f(\boldsymbol{D}_{kh}^{O})$. Under convexity condition on $f$, we have \|\big{(}\mathbf{I}_{p}-e^{-s\mathbf{H}_{k}}\big{)}\mathbf{H}_{k}^{-1}\|\leq s for every $s>0$. Therefore, conditioning with respect to $\boldsymbol{D}_{kh}^{O}$ and using the inequality $(a+b)^{4}\leq 8(a^{4}+b^{4})$, for every $t\in[kh,(k+1)h[$ we get

This inequality, in conjunction with (42) and (43) yields

To bound the last expectation, we use the fact that $\boldsymbol{D}_{kh}^{O}$ equals $\bar{\boldsymbol{\vartheta}}^{(k,h)}$ in distribution, and the next lemma (the proof of which is provided in the supplementary material).

If $p\geq 2$, $T\geq 4/(3M)$ and $h\leq 1/(8M)$, then the iterates of the LMCO algorithm satisfy \mathbf{E}\big{[}\big{(}\sum_{k=0}^{K-1}\|\nabla f(\bar{\boldsymbol{\vartheta}}^{(k,h)})\|_{2}^{2}\big{)}^{2}\big{]}\leq\frac{32}{3}\big{(}{TMp}/h\big{)}^{2}.

Acknowledgments

The work of the author was partially supported by the grant Investissements d’Avenir (ANR-11-IDEX-0003/Labex Ecodec/ANR-11-LABX-0047).

References