Lattice Gaussian Sampling by Markov Chain Monte Carlo: Bounded Distance Decoding and Trapdoor Sampling

TL;DR

This paper advances MCMC-based lattice Gaussian sampling, deriving spectral gaps, analyzing decoding performance, and proposing new algorithms that improve convergence and enable parallel implementation for practical cryptographic applications.

Contribution

It introduces new spectral gap bounds, analyzes bounded distance decoding trade-offs, and proposes the MTMK algorithm for faster convergence in lattice Gaussian sampling.

Findings

Independent MHK converges faster due to spectral gap bounds.

Decoding performance shows a trade-off between radius and complexity.

MTMK enhances convergence rate and supports parallel implementation.

Abstract

Sampling from the lattice Gaussian distribution plays an important role in various research fields. In this paper, the Markov chain Monte Carlo (MCMC)-based sampling technique is advanced in several fronts. Firstly, the spectral gap for the independent Metropolis-Hastings-Klein (MHK) algorithm is derived, which is then extended to Peikert's algorithm and rejection sampling; we show that independent MHK exhibits faster convergence. Then, the performance of bounded distance decoding using MCMC is analyzed, revealing a flexible trade-off between the decoding radius and complexity. MCMC is further applied to trapdoor sampling, again offering a trade-off between security and complexity. Finally, the independent multiple-try Metropolis-Klein (MTMK) algorithm is proposed to enhance the convergence rate. The proposed algorithms allow parallel implementation, which is beneficial for practical applications.

Tables

Table 1. TABLE I: Values of ϑ 3 subscript italic-ϑ 3 \vartheta_{3} .

${\vartheta }_3\left(1\right)$	${\sum }_{n=-\infty }^{+\infty }{e}^{-1\pi n^2}$	$\frac{\sqrt[4]{\pi }}{\Gamma \left(\frac{3}{4}\right)}$	1.087
${\vartheta }_3\left(2\right)$	${\sum }_{n=-\infty }^{+\infty }{e}^{-2\pi n^2}$	$\frac{\sqrt[4]{6\pi +4\sqrt{2}\pi }}{2\Gamma \left(\frac{3}{4}\right)}$	1.0039
${\vartheta }_3\left(3\right)$	${\sum }_{n=-\infty }^{+\infty }{e}^{-3\pi n^2}$	$\frac{\sqrt[4]{27\pi +18\sqrt{3}\pi }}{3\Gamma \left(\frac{3}{4}\right)}$	1.00037
${\vartheta }_3\left(4\right)$	${\sum }_{n=-\infty }^{+\infty }{e}^{-4\pi n^2}$	$\frac{\sqrt[4]{8\pi }+2\sqrt[4]{\pi }}{4\Gamma \left(\frac{3}{4}\right)}$	1.0002
${\vartheta }_3\left(5\right)$	${\sum }_{n=-\infty }^{+\infty }{e}^{-5\pi n^2}$	$\frac{\sqrt[4]{225\pi +100\sqrt{5}\pi }}{5\Gamma \left(\frac{3}{4}\right)}$	1.0001

Equations

$$\Lambda =\{\mathbf{Bx}:\mathbf{x}\in {\mathbb{Z}}^n\},$$

$${\rho }_{\sigma ,\mathbf{c}}(\mathbf{z})=e^{-\frac{\Vert \mathbf{z}-\mathbf{c}{\Vert }^2}{2{\sigma }^2}},$$

$$D_{\Lambda ,\sigma ,\mathbf{c}}(\mathbf{x})=\frac{{\rho }_{\sigma ,\mathbf{c}}(\mathbf{Bx})}{{\rho }_{\sigma ,\mathbf{c}}(\Lambda )}=\frac{e^{-\frac{1}{2{\sigma }^2}\parallel \mathbf{Bx}-\mathbf{c}{\parallel }^2}}{{\sum }_{\mathbf{x}\in {\mathbb{Z}}^n}{e}^{-\frac{1}{2{\sigma }^2}\parallel \mathbf{Bx}-\mathbf{c}{\parallel }^2}}$$

$$\mathbf{c}=\mathbf{Bx}+\mathbf{w}$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}\in {\mathbb{Z}}^n}{\arg \min }\Vert \mathbf{c}-\mathbf{Bx}{\Vert }^2$$

$$\alpha =\text{min}\left\{1,\frac{\pi (\mathbf{y})q(\mathbf{y},\mathbf{x})}{\pi (\mathbf{x})q(\mathbf{x},\mathbf{y})}\right\},$$

$$P(\mathbf{x},\mathbf{y})=\{\begin{array}{l}q(\mathbf{x},\mathbf{y})\alpha \text{if}\mathbf{y}\ne \mathbf{x},\\ 1-{\sum }_{\mathbf{z}\ne \mathbf{x}}q(\mathbf{x},\mathbf{z})\alpha \text{if}\mathbf{y}=\mathbf{x}.\end{array}$$

$$t_{\text{mix}}(ϵ)=\text{min}\{t:\text{max}\Vert P^t(\mathbf{x},\cdot )-\pi (\cdot ){\Vert }_{TV}\le ϵ\},$$

$$P(y_i|{\overline{\mathbf{y}}}_{[-i]})$$

$$=\frac{e^{-\frac{1}{2{\sigma }^2}\parallel {\overline{\mathbf{c}}}^{\prime }-\overline{\mathbf{R}}\overline{\mathbf{y}}{\parallel }^2}}{{\sum }_{y_i\in \mathbb{Z}}{e}^{-\frac{1}{2{\sigma }^2}\parallel {\overline{\mathbf{c}}}^{\prime }-\overline{\mathbf{R}}\overline{\mathbf{y}}{\parallel }^2}}$$

$$=\frac{e^{-\frac{1}{2{\sigma }^2}\parallel c_i^{\prime }-{\sum }_{j=i}^n{r}_{i,j}{y}_j{\parallel }^2}}{{\sum }_{y_i\in \mathbb{Z}}{e}^{-\frac{1}{2{\sigma }^2}\parallel c_i^{\prime }-{\sum }_{j=i}^n{r}_{i,j}{y}_j{\parallel }^2}}$$

$$=\frac{e^{-\frac{1}{2{\sigma }_i^2}\parallel y_i-{\tilde{y}}_i{\parallel }^2}}{{\sum }_{y_i\in \mathbb{Z}}{e}^{-\frac{1}{2{\sigma }_i^2}\parallel y_i-{\tilde{y}}_i{\parallel }^2}}$$

$$=D_{\mathbb{Z},{\sigma }_i,{\tilde{y}}_i}(y_i),$$

$$q(\mathbf{x},\mathbf{y})$$

$$=\frac{{\rho }_{\sigma ,\mathbf{c}}(\mathbf{By})}{{\prod }_{i=1}^n{\rho }_{{\sigma }_{n+1-i},{\tilde{y}}_{n+1-i}}(\mathbb{Z})}$$

$$=q(\mathbf{y}),$$

$$\alpha$$

$$\gamma \ge \delta \triangleq \frac{{\rho }_{\sigma ,\mathbf{c}}(\mathbf{\Lambda })}{{\prod }_{i=1}^n{\rho }_{{\sigma }_i}(\mathbb{Z})}.$$

$$P(\mathbf{x},\mathbf{y})$$

$$w(\mathbf{x})=\frac{\pi (\mathbf{x})}{q(\mathbf{x})}.$$

$$P(\mathbf{x},\mathbf{y})$$

$$w({\mathbf{x}}_1)\ge w({\mathbf{x}}_2)\ge \cdots \ge w({\mathbf{x}}_{\infty }).$$

$$\mathbf{P}=\left[\begin{array}[]{ccccc}q(\mathbf{x}_{1})+\eta_{1}&\frac{\pi(\mathbf{x}_{2})}{w(\mathbf{x}_{1})}&\frac{\pi(\mathbf{x}_{3})}{w(\mathbf{x}_{1})}&\cdots&\frac{\pi(\mathbf{x}_{\infty})}{w(\mathbf{x}_{1})}\\ q(\mathbf{x}_{1})&q(\mathbf{x}_{2})+\eta_{2}&\frac{\pi(\mathbf{x}_{3})}{w(\mathbf{x}_{2})}&\cdots&\frac{\pi(\mathbf{x}_{\infty})}{w(\mathbf{x}_{2})}\\ q(\mathbf{x}_{1})&q(\mathbf{x}_{2})&q(\mathbf{x}_{3})+\eta_{3}&\cdots&\frac{\pi(\mathbf{x}_{\infty})}{w(\mathbf{x}_{3})}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ q(\mathbf{x}_{1})&q(\mathbf{x}_{2})&q(\mathbf{x}_{3})&\cdots&q(\mathbf{x}_{\infty})\end{array}\right]$$

$${\eta }_j=\sum _{i=j}^{\infty }\left(q({\mathbf{x}}_i)-\frac{\pi ({\mathbf{x}}_i)}{w({\mathbf{x}}_j)}\right)$$

$$\mathbf{P}=\mathbf{G}+\mathbf{e}{\mathbf{q}}^T,$$

$$\mathbf{G}=\left[\begin{array}[]{cccc}\eta_{1}&\frac{\pi(\mathbf{x}_{2})}{w(\mathbf{x}_{1})}-q(\mathbf{x}_{2})&\cdots&\frac{\pi(\mathbf{x}_{\infty})}{w(\mathbf{x}_{1})}-q(\mathbf{x}_{\infty})\\ 0&\eta_{2}&\cdots&\frac{\pi(\mathbf{x}_{\infty})}{w(\mathbf{x}_{2})}-q(\mathbf{x}_{\infty})\\ \vdots&\vdots&\ddots&\vdots\\ 0&0&\cdots&0\end{array}\right].$$

$${\tau }_1={\eta }_1$$

$$1>|{\eta }_1|\ge |{\eta }_2|\ge \cdots >0,$$

$${\tau }_1$$

$$=\sum _{i=1}^{\infty }q({\mathbf{x}}_i)-\frac{1}{w({\mathbf{x}}_1)}\cdot \sum _{i=1}^{\infty }\pi ({\mathbf{x}}_i)$$

$$=1-\frac{1}{w({\mathbf{x}}_1)}$$

$$=1-\frac{q({\mathbf{x}}_1)}{\pi ({\mathbf{x}}_1)}.$$

$$\frac{q(\mathbf{x})}{\pi (\mathbf{x})}\ge \delta$$

$$\Vert P^t(\mathbf{x},\cdot )-D_{\Lambda ,\sigma ,\mathbf{c}}(\cdot ){\Vert }_{TV}\le (1-\delta {)}^t,$$

$$t_{\text{mix}}(ϵ)=\frac{\text{ln}ϵ}{\text{ln}(1-\delta )}<(-\text{ln}ϵ)\cdot \left(\frac{1}{\delta }\right),ϵ<1$$

Full text

Keywords: Lattice decoding, lattice Gaussian sampling, Markov chain Monte Carlo, bounded distance decoding, large-scale MIMO detection, trapdoor sampling.

I Introduction

Nowadays, lattice Gaussian sampling has drawn a lot of attention in various research fields. In mathematics, Banaszczyk was the first to apply it to prove the transference theorems for lattices [1]. In coding, lattice Gaussian distribution was employed to obtain the full shaping gain for lattice coding [2, 3], and to achieve the capacity of the Gaussian channel [4]. It was also used to achieve information-theoretic security in the Gaussian wiretap channel [5, 6] and in the bidirectional relay channel [7], respectively. In cryptography, the lattice Gaussian distribution has become a central tool in the construction of many primitives [8, 9, 10]. Specifically, lattice Gaussian sampling lies at the core of signature schemes in the Gentry, Peikert and Vaikuntanathan (GPV) paradigm [11]. Furthermore, lattice Gaussian sampling with a suitable variance allows to solve the closest vector problem (CVP) and the shortest vector problem (SVP) [12, 13].

However, in sharp contrast to the continuous Gaussian density, it is by no means trivial even to sample from a low-dimensional discrete Gaussian distribution. For some special lattices, there are rather efficient algorithms for Gaussian sampling [4, 14]. As the default sampling algorithm for general lattices, Klein’s algorithm [15] only works when the standard deviation $\sigma=\sqrt{\omega(\text{log}\ n)}\cdot\text{max}_{1\leq i\leq n}\|\mathbf{\widehat{b}}_{i}\|$ [11], where $\omega(\text{log}\ n)$ is a superlogarithmic function, $n$ denotes the lattice dimension and $\mathbf{\widehat{b}}_{i}$’s are the Gram-Schmidt vectors of the lattice basis $\mathbf{B}$. Peikert gave an efficient lattice Gaussian sampler in [16] for parallel implementation, which however requires larger values of $\sigma$. On the other hand, the lattice Gaussian sampling algorithm proposed by Aggarwal et al. in [12, 13] to solve CVP and SVP has a lower bound $2^{n}$ on both space and time complexity; it actually obtains samples for small $\sigma$ by combining original samples for $\sigma=\sqrt{\omega(\text{log}\ n)}\cdot\text{max}_{1\leq i\leq n}\|\mathbf{\widehat{b}}_{i}\|$. Although the algorithm in [17] provides a trade-off between (exponential) time and space complexity, its complexity is still too high to be practical.

In order to sample from a target lattice Gaussian distribution with arbitrary $\sigma>0$, Markov chain Monte Carlo (MCMC) methods were introduced in [18]. In principle, it randomly generates the next Markov state conditioned on the previous one; after the burn-in time, which is normally measured by the mixing time, the Markov chain will step into a stationary distribution, when samples from the target distribution can be obtained [19]. It has been demonstrated that Gibbs sampling, which employs univariate conditional sampling to build the chain, yields an ergodic Markov chain [20]. In [18], we proposed an independent Metropolis-Hastings (MH) algorithm incorporating Klein’s algorithm (namely, the independent MHK algorithm) to generate a proposal distribution, which is shown to be uniformly ergodic (converging exponentially fast to the stationary distribution). Meanwhile, the associated convergence rate of the Markov chain is derived, resulting in a tractable estimation of the mixing time. Differently from the algorithms of [12, 13, 17], the independent MHK sampling algorithm only requires polynomial space. In this paper, we advance the state of the art of MCMC-based lattice Gaussian sampling in several fronts.

Firstly, we refine the analysis and extend the independent MHK algorithm of [18]. We obtain the spectral gap of the transition matrix and demonstrate uniformly ergodicity. We extend the independent MH algorithm to a version where Peikert’s algorithm [16] is used to generate the proposal distribution. We then compare these MCMC algorithms with rejection sampling from statistics. By deriving their rates of convergence, we show the advantage of the independent MHK. Rejection sampling achieves the same convergence rate only if its normalizing constant is carefully chosen, which is generally rather difficult.

Secondly, we apply the independent MHK algorithm to bounded distance decoding (BDD). BDD is a variant of the CVP where the input is within a certain distance to the lattice. With a careful selection of the standard deviation $\sigma$ during the sampling process, we improve the result of Klein from $\eta=O(1/n)$ to $\eta=O(\sqrt{\log n}/n)$ in terms of $\eta-$BDD111In $\eta$-BDD ($\eta<1/2$), we are given a lattice basis $\mathbf{B}$ and a query point $\mathbf{c}$, and we are asked to find a lattice point within distance $\eta\cdot\lambda_{1}$ from the target, where $\lambda_{1}$ denotes the first minimum of the lattice.. References [21, 22] achieved a larger value $\eta=O(\sqrt{\log n/n})$, at the expense of a pre-processing stage where Gaussian samples are taken from the dual lattice with standard deviation $\sigma$ equal to its smoothing parameter. However, sampling at the smoothing parameter is in general a difficult problem with no efficient solutions nowadays. For algorithms of general SVP/CVP such as enumeration and sieving, we refer the readers to the comprehensive survey [23].

Thirdly, we examine the impact of MCMC to trapdoor sampling in the GPV paradigm. In cryptographic applications, the standard deviation $\sigma$ of the sampler is the main parameter governing the security level. Namely, the smaller $\sigma$, the higher security. This is because for a signature system to be secure, it must be hard for an adversary to find lattice points of length about $\sigma\sqrt{n}$. We show that, at moderate costs of increased complexity, MCMC is able to sample with smaller $\sigma$, thereby increasing the security level relative to Klein’s algorithm [11] and Peikert’s algorithm [16].

Finally, to improve the convergence rate of the Markov chain, the independent multiple-try Metropolis-Klein (MTMK) algorithm is proposed, which fully exploits the trial samples generated from the proposal distribution. Uniform ergodicity is demonstrated and the enhanced convergence rate is also given. Since independent MHK is only a special case of independent MTMK, the decoding performance can also be improved due to the usage of trial samples. The proposed sampling algorithm allows a parallel implementation and is easily adopted to MIMO detection to achieve near-optimal performance. With the development of 5G, the demand for large-scale MIMO systems will increase in the next decade, which has triggered research activities towards low complexity decoding algorithms for large-scale MIMO detection [24, 25, 26]. Therefore, there has been considerable interest in MCMC sampling for the efficient decoding of MIMO systems [27, 28, 29, 30, 31, 32].

The rest of this paper is organized as follows. Section II introduces the lattice Gaussian distribution and briefly reviews the basics of MCMC. In Section III, we derive the spectral gaps of the Markov chains associated with independent MHK and rejection sampling-based lattice Gaussian sampling, and show their uniform ergodicity as well as convergence rates. An extension to Peikert’s algorithm is also given. Then, the decoding complexity of BDD using independent MHK algorithm is derived in Section IV. Section V addresses trapdoor sampling using MCMC. In Section VI, the independent MTMK algorithm is proposed to further strength the convergence performance. Simulation results for MIMO detection are presented in Section VII. Finally, Section VIII concludes the paper.

Notation: Matrices and column vectors are denoted by upper and lowercase boldface letters, and the transpose, inverse, pseudoinverse of a matrix $\mathbf{B}$ by $\mathbf{B}^{T},\mathbf{B}^{-1},$ and $\mathbf{B}^{{\dagger}}$, respectively. $\mathbf{I}$ denotes the identity matrix. We use $\mathbf{b}_{i}$ for the $i$th column of the matrix $\mathbf{B}$, $\mathbf{\widehat{b}}_{i}$ for the $i$th Gram-Schmidt vector of the matrix $\mathbf{B}$, $b_{i,j}$ for the entry in the $i$th row and $j$th column of the matrix $\mathbf{B}$. A symmetric matrix $\mathbf{B}$ is written as $\mathbf{B}\succ\mathbf{0}$ if it is positive definite. Similarly, we say $\mathbf{B}_{1}\succ\mathbf{B}_{2}$ if $(\mathbf{B}_{1}-\mathbf{B}_{2})\succ\mathbf{0}$. $\lceil x\rfloor$ denotes rounding to the integer closest to $x$. If $x$ is a complex number, $\lceil x\rfloor$ rounds the real and imaginary parts separately. In addition, we use the standard small omega notation $\omega(\cdot)$, i.e., $|\omega(g(n))|>k\cdot|g(n)|$ for every fixed positive number $k>0$. Finally, in this paper, the computational complexity is measured by the number of Markov moves.

II Preliminaries

In this section, we introduce the background and mathematical tools needed to describe and analyze the proposed lattice Gaussian sampling algorithms.

II-A Lattice Gaussian Distribution

Let matrix $\mathbf{B}=[\mathbf{b}_{1},\ldots,\mathbf{b}_{n}]\subset\mathbb{R}^{n}$ consist of $n$ linearly independent column vectors. The $n$-dimensional lattice $\Lambda$ generated by $\mathbf{B}$ is defined by

[TABLE]

where $\mathbf{B}$ is called the lattice basis. We define the Gaussian function centered at $\mathbf{c}\in\mathbb{R}^{n}$ for standard deviation $\sigma>0$ as

[TABLE]

for all $\mathbf{z}\in\mathbb{R}^{n}$. When $\mathbf{c}$ or $\sigma$ are not specified, we assume that they are $\mathbf{0}$ and $1$ respectively. Then, the discrete Gaussian distribution over $\Lambda$ is defined as

[TABLE]

for all $\mathbf{x}\in\mathbb{Z}^{n}$, where $\rho_{\sigma,\mathbf{c}}(\Lambda)\triangleq\sum_{\mathbf{\mathbf{Bx}}\in\Lambda}\rho_{\sigma,\mathbf{c}}(\mathbf{Bx})$ is just a scaling to obtain a probability distribution. We remark that this definition differs slightly from the one in [8], where $\sigma$ is scaled by a constant factor $\sqrt{2\pi}$ (i.e., $s=\sqrt{2\pi}\sigma$). In fact, the discrete Gaussian resembles a continuous Gaussian distribution, but is only defined over a lattice. It has been shown that discrete and continuous Gaussian distributions share similar properties, if the flatness factor is small [5].

II-B Decoding by Sampling

Consider the decoding of an $n\times n$ real-valued system. The extension to the complex-valued system is straightforward [33]. Let $\mathbf{x}\in\mathbb{Z}^{n}$ denote the transmitted signal. The corresponding received signal $\mathbf{c}$ is given by

[TABLE]

where $\mathbf{w}$ is the noise vector with zero mean and variance $\sigma_{w}^{2}$, $\mathbf{B}$ is an $n\times n$ full column-rank matrix of channel coefficients. Typically, the conventional maximum likelihood (ML) reads

[TABLE]

where $\|\cdot\|$ denotes the Euclidean norm. Clearly, ML decoding corresponds to the CVP. If the received signal $\mathbf{c}$ is the origin, then ML decoding reduces to SVP.

Intuitively, the CVP given in (5) can be solved by the lattice Gaussian sampling. Since the distribution is centered at the query point $\mathbf{c}$, the closest lattice point $\mathbf{Bx}$ to $\mathbf{c}$ is assigned the largest sampling probability. Therefore, by multiple samplings, the solution of CVP is the most likely to be returned. It has been demonstrated that lattice Gaussian sampling is equivalent to CVP via a polynomial-time dimension-preserving reduction [34]. Meanwhile, by adjusting the sample size, the sampling decoder enjoys a flexible trade-off between performance and complexity.

In [15], Klein introduced an algorithm which performs sampling from a Gaussian-like distribution (see Algorithm 1). It is shown in [15, 33, 35] that Klein’s algorithm is able to find the closest lattice point when it is close to the input vector: this technique is known as BDD in coding literature, which corresponds to a restricted variant of CVP.

II-C Classical MH Algorithms

In [36], the original Metropolis algorithm was extended to a more general scheme known as the Metropolis-Hastings (MH) algorithm. In particular, let us consider a target invariant distribution $\pi$ together with a proposal distribution $q(\mathbf{x},\mathbf{y})$. Given the current state $\mathbf{x}$ for Markov chain $\mathbf{X}_{t}$, a state candidate $\mathbf{y}$ for the next Markov move $\mathbf{X}_{t+1}$ is generated from the proposal distribution $q(\mathbf{x},\mathbf{y})$. Then the acceptance ratio $\alpha$ is computed by

[TABLE]

and $\mathbf{y}$ will be accepted as the new state with probability $\alpha$. Otherwise, $\mathbf{x}$ will be retained. In this way, a Markov chain $\{\mathbf{X}_{0},\mathbf{X}_{1},\ldots\}$ is established with the transition probability $P(\mathbf{x},\mathbf{y})$ as follows:

[TABLE]

It is interesting that in MH algorithms, the proposal distribution $q(\mathbf{x},\mathbf{y})$ can be any fixed distribution from which we can conveniently draw samples. Therefore, there is large freedom in the choice of $q(\mathbf{x},\mathbf{y})$ but it is challenging to find a suitable one with satisfactory convergence. In fact, Gibbs sampling can be viewed as a special case of the MH algorithm, whose proposal distribution is a univariate conditional distribution.

As an important parameter to measure the time required by a Markov chain to get close to its stationary distribution, the mixing time is defined as [19]

[TABLE]

where $P^{t}(\mathbf{x},\cdot)$ denotes a row of the transition matrix $\mathbf{P}$ for $t$ Markov moves and $\|\cdot\|_{TV}$ represents the total variation distance.

II-D Independent MHK Algorithm

From the MCMC perspective, lattice Gaussian distribution can be viewed as a complex target distribution lacking direct sampling methods. In order to obtain samples from $D_{\Lambda,\sigma,\mathbf{c}}(\mathbf{x})$, the independent MHK sampling was proposed in [18]. Specifically, a state candidate $\mathbf{y}$ for the next Markov move $\mathbf{X}_{t+1}$ is generated by Klein’s algorithm, via the following backward one-dimensional conditional sampling (for $i=n,n-1,\ldots,1$):

[TABLE]

where $\widetilde{y}_{i}=\frac{c^{\prime}_{i}-\sum^{n}_{j=i+1}r_{i,j}y_{j}}{r_{i,i}}$, $\sigma_{i}=\frac{\sigma}{|r_{i,i}|}$, $\mathbf{c}^{\prime}=\mathbf{Q}^{{\dagger}}\mathbf{c}$ and $\mathbf{B}=\mathbf{QR}$ by QR decomposition with $\|\mathbf{\widehat{b}}_{i}\|=|r_{i,i}|$. Note that $\overline{\mathbf{y}}_{[-i]}=[y_{i+1},\ldots,y_{n}]$, $\mathbf{\overline{R}}$, $\mathbf{\overline{c}^{\prime}}$ and $\mathbf{\overline{y}}$ are the $(n-i+1)$ segments of $\mathbf{R}$, $\mathbf{c}^{\prime}$ and $\mathbf{y}$ respectively (i.e., $\mathbf{\overline{R}}$ is a $(n-i+1)\times(n-i+1)$ submatrix of $\mathbf{R}$ with $r_{i,i}$ to $r_{n,n}$ in the diagonal).

Given the current state $\mathbf{x}$, the proposal distribution $q(\mathbf{x},\mathbf{y})$ in the independent MHK sampling is given by

[TABLE]

where the proposal distribution $q(\mathbf{x},\mathbf{y})$ is actually independent of $\mathbf{x}$. Therefore, the connection between two consecutive Markov moves is only due to the decision stage.

With the state candidate $\mathbf{y}$, the acceptance ratio $\alpha$ is obtained by substituting (10) into (6)

[TABLE]

where $\widetilde{x}_{i}=\frac{c^{\prime}_{i}-\sum^{n}_{j=i+1}r_{i,j}x_{j}}{r_{i,i}}$ and we note that $\pi=D_{\Lambda,\sigma,\mathbf{c}}$ in (6) (these notations will be followed throughput the context). The sampling procedure is summarized in Algorithm 2. Note that the initial state $\mathbf{x}_{0}$ for $\mathbf{X}_{0}$ can be chosen from $\mathbb{Z}^{n}$ arbitrarily or from the output of a suboptimal algorithm.

Thanks to the celebrated coupling technique, the uniformly ergodicity was demonstrated in [18]. Nevertheless, the spectral gap of the transition matrix, which serves as an important metric for the mixing time of the underlying Markov chain, has not been determined yet.

III Convergence Analysis

In this section, the spectrum of the Markov chain induced by independent MHK sampling is analyzed, followed by the extensions to Peikert’s algorithm and rejection sampling. As a common way to evaluate the mixing time, the spectral gap $\gamma=1-|\tau_{1}|>0$ of the transition matrix is preferred for convergence analysis in MCMC [19]. Here, $\tau_{1}$ represents the second largest eigenvalue in magnitude of the transition matrix $\mathbf{P}$ [37].

III-A Spectral Gap of Independent MHK Algorithm

Theorem 1**.**

Given the invariant lattice Gaussian distribution $D_{\Lambda,\sigma,\mathbf{c}}$, the Markov chain induced by independent MHK sampling exhibits a spectral gap

[TABLE]

Proof.

From (10) and (11), the transition probability $P(\mathbf{x},\mathbf{y})$ of each Markov move in the independent MHK sampling is given by

[TABLE]

For notational simplicity, we define the importance weight $w(\mathbf{x})$ as

[TABLE]

Then the transition probability can be rewritten as

[TABLE]

Without loss of generality, we label the countably infinite state space $\Omega=\mathbb{Z}^{n}$ as $\Omega=\{\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{\infty}\}$ and assume that these states are sorted according to their importance weights, namely,

[TABLE]

From (15) and (16), the transition matrix $\mathbf{P}$ of the Markov chain can be exactly expressed as

[TABLE]

where

[TABLE]

stands for the probability of being rejected in the decision stage with the current state $\mathbf{x}_{j}$ for $\mathbf{X}_{t}$.

Let $\mathbf{q}=[q(\mathbf{x}_{1}),q(\mathbf{x}_{2}),\ldots]^{T}$ denote the vector of proposal probabilities. Then by decomposition, it follows that

[TABLE]

where $\mathbf{e}=[1,1,\ldots]^{T}$ and $\mathbf{G}$ is an upper triangular matrix of the form

[TABLE]

It is well-known that for a Markov chain, the largest eigenvalue of the transition matrix $\mathbf{P}$ always equals 1. Here, as $\mathbf{e}$ is a common right eigenvector for both $\mathbf{P}$ and $\mathbf{P}-\mathbf{G}$, it naturally corresponds to the largest eigenvalue 1. Meanwhile, since the rank of $\mathbf{P}-\mathbf{G}$ is 1, the other eigenvalues of $\mathbf{G}$ are exactly the same as those of $\mathbf{P}$.

Thanks to the ascending order in (16), it is easy to verify that the spectral radius $\tau_{1}$ is exactly given by

[TABLE]

and

[TABLE]

thereby raising the interest of identifying the value of $\tau_{1}$.

Therefore, according to (17) and (19), we can easily get that

[TABLE]

In other words, the spectral gap $1-\tau_{1}$ is exactly captured by the ratio $q(\mathbf{x}_{1})/\pi(\mathbf{x}_{1})$. Next, we invoke the following Lemma to lower bound the ratio $q(\mathbf{x})/\pi(\mathbf{x})$ for $\mathbf{x}\in\mathbb{Z}^{n}$.

Lemma 1** ([18]).**

In the independent MHK algorithm

[TABLE]

for all $\mathbf{x}\in\mathbb{Z}^{n}$, where $\delta$ is defined in (12).

The proof is completed by combining (21) and (22).

∎

By using the coupling technique, it is shown in [18] that the Markov chain converges exponentially fast to the stationary distribution in total variational distance:

[TABLE]

The mixing time of the Markov chain is given by

[TABLE]

which is proportional to $1/\delta$, and becomes $O(1)$ if $\delta\to 1$.

III-B Extension to Peikert’s Algorithm

Klein’s sampling algorithm is a randomized variant of Babai’s nearest-plane algorithm for lattice decoding [38]. Babai also proposed a simpler decoding scheme by direct rounding222In communications, Babai’s nearest-plane algorithm is known as successive interference cancelation (SIC) while the direct rounding algorithm is referred to as zero-forcing (ZF)., which was further randomized by Peikert in [16]. Although Peikert’s algorithm requires a higher value of $\sigma$, it is parallelizable and can be more attractive in practical implementation. In fact, Peikert’s algorithm can also be incorporated into the Metropolis-Hastings algorithm to overcome the limitation of $\sigma$.

Specifically, given the standard deviation $\sigma>0$ and a basis $\mathbf{B}$, one chooses a positive definite matrix $\Sigma_{1}=r^{2}\cdot\mathbf{B}\mathbf{B}^{T}\prec\Sigma=\sigma^{2}\cdot\mathbf{I}$ for $r>0$ (i.e., $\Sigma_{2}=\Sigma-\Sigma_{1}$ is positive definite). Then, the proposed sample $\mathbf{z}\in\Lambda$ is taken from the distribution $\mathbf{c}+\mathbf{z}^{\prime}+D_{\Lambda-\mathbf{c}-\mathbf{z}^{\prime},\sqrt{\Sigma_{1}}}$, where $\mathbf{z}^{\prime}\in\mathbb{R}^{n}$ is sampled from the continuous distribution $D_{\sqrt{\Sigma_{2}}}$. Note the lattice Gaussian distribution $D_{\Lambda-\mathbf{c}-\mathbf{z}^{\prime},\sqrt{\Sigma_{1}}}$ is expressed as

[TABLE]

with

[TABLE]

The joint probability distribution of $\mathbf{z}\in\Lambda$ and $\mathbf{z}^{\prime}\in\mathbb{R}^{n}$ is given by

[TABLE]

where (a) is due to the symmetry of $\rho_{\sqrt{\Sigma_{1}}}$, and (b) follows from [16, Fact 2.1] with positive definite matrix $\Sigma_{3}^{-1}=\Sigma_{1}^{-1}+\Sigma_{2}^{-1}$ and $\mathbf{c}^{\prime}=\Sigma_{3}\Sigma_{1}^{-1}(\mathbf{z}-\mathbf{c})$. Consequently, the marginal distribution of $\mathbf{z}$ is

[TABLE]

As $\mathbf{z}=\mathbf{Bx}$ for $\mathbf{x}\in\mathbb{Z}^{n}$, we have

[TABLE]

Clearly, $P(\mathbf{\cdot})$ can be used as a proposal distribution $q(\cdot)$ in the MH algorithm to obtain the state candidate $\mathbf{y}\in\mathbb{Z}^{n}$. In this case, the acceptance ratio $\alpha$ can be calculated by

[TABLE]

followed by a decision to accept $\mathbf{X}_{t+1}=\mathbf{y}$ or not. To summarize, its operation procedure is shown in Algorithm 3.

Lemma 2**.**

In the independent MH algorithm using Peikert’s algorithm, there exists a constant $\delta^{\prime}>0$ such that

[TABLE]

for all $\mathbf{x}\in\mathbb{Z}^{n}$, where

[TABLE]

Proof.

To start with, we have

[TABLE]

where inequality $(c)$ comes from the fact that $\rho_{\sqrt{\Sigma}}(\Lambda-\mathbf{c})\leq\rho_{\sqrt{\Sigma}}(\Lambda)$.

The Lemma is proven by showing that

[TABLE]

Here, $(d)$ and $(e)$ follow from the properties of determinant that

[TABLE]

and

[TABLE]

respectively, for square matrices $\mathbf{A}$ and $\mathbf{B}$ of equal sizes.

∎

To satisfy the condition that $\sigma^{2}\mathbf{I}\succ r^{2}\cdot\mathbf{B}\mathbf{B}^{T}$, we require

[TABLE]

where $s_{1}(\mathbf{B})$ denotes the largest singular value of the basis $\mathbf{B}$. It is readily verified that

[TABLE]

Lemma 3**.**

For independent MH samplings based on Peikert’s algorithm and on Klein’s algorithm, the following relation holds:

[TABLE]

Proof.

According to (12) and (32), in order to show $\delta^{\prime}\leq\delta$, we need to prove that

[TABLE]

Next, by recalling the Jacobi theta function $\vartheta_{3}(\tau)=\sum^{+\infty}_{n=-\infty}e^{-\pi\tau n^{2}}$ with $\tau>0$, we have

[TABLE]

and the left-hand side of (40)

[TABLE]

where $(f)$ utilizes the symmetry property of Theta series for isodual lattice $\mathbb{Z}$

[TABLE]

Moreover, as $\vartheta_{3}(\tau)$ is monotone decreasing with $\tau$, the following relation holds:

[TABLE]

due to (37) and (38).

Hence, we finally have that the left-hand side of (40)

[TABLE]

thus completing the proof.

∎

Similarly to independent MHK, it is easy to verify that the proposed algorithm is also uniformly ergodic.

Theorem 2**.**

Given $D_{\Lambda,\sigma,\mathbf{c}}$, the Markov chain induced by independent MH sampling using Peikert’s algorithm converges exponentially fast:

[TABLE]

By Lemma 3, we can see that the independent MH sampling based on Peikert’s algorithm converges slower than that based on Klein’s algorithm. This is numerically confirmed in Fig. 2 for checkerboard lattice $D_{4}$, where a comparison of the coefficients $1/\delta$ and $1/\delta^{\prime}$ is given. Clearly, in the whole range of $r$, the independent MH-Peikert sampling requires more iterations than independent MHK.

III-C Extension to Rejection Sampling

The classic rejection sampling is able to generate independent samples from the target distribution, but requires a normalizing constant for the application of a proposal distribution [39]. Given the target distribution $\pi(\mathbf{x})=D_{\Lambda,\sigma,\mathbf{c}}(\mathbf{x})$, its operation consists of the following three steps:

1)  Generate a candidate sample $\mathbf{y}$ from distribution $q(\mathbf{y})$ using Klein’s algorithm or Peikert’s algorithm.

2)  Calculate a normalizing constant $\omega_{0}$ such that

[TABLE]

for all $\mathbf{x}\in\mathbb{Z}^{n}$.

3)  Output $\mathbf{y}$ with probability

[TABLE]

and otherwise repeat.

Generally, rejection sampling is not directly comparable with MCMC sampling as it requires the normalizing constant $\omega_{0}$ for calibrating, which is not realistic in many cases of interest. Nevertheless, with a certain choice of $\omega_{0}$, it is possible to interpret it as a particular MCMC algorithm.

Definition 1**.**

Given the target distribution $\pi(\mathbf{x})=D_{\Lambda,\sigma,\mathbf{c}}(\mathbf{x})$, the Markov chain arising from the above rejection sampler with $\omega_{0}\geq\pi(\mathbf{x})/q(\mathbf{x})$ for all $\mathbf{x}\in\mathbb{Z}^{n}$ is reversible, irreducible and aperiodic, with transition probability

[TABLE]

Clearly, the algorithm based on rejection sampling converges when the first acceptance takes place. The samples after the acceptance are naturally independently and identically distributed (i.i.d.). Similarly to the setting in (16), the transition matrix $\mathbf{P}_{\text{r}}$ of this Markov chain is exactly given by

[TABLE]

which can be further decomposed into

[TABLE]

where

[TABLE]

[TABLE]

and

[TABLE]

Here, $P_{\text{accept}}$ and $P_{\text{reject}}$ denote the acceptance and rejection probabilities of a new candidate in the next move.

Similarly to the analysis of independent MHK, we have the following Lemma, whose proof is omitted due to simplicity.

Lemma 4**.**

The eigenvalues $\eta_{i}$’s of the transition matrix $\mathbf{P}_{\text{r}}$ satisfy that

[TABLE]

with

[TABLE]

for $i=1,\ldots,\infty$.

Furthermore, we arrive at the following Theorem.

Theorem 3**.**

Given the invariant lattice Gaussian distribution $\pi=D_{\Lambda,\sigma,\mathbf{c}}$, the Markov chain induced by rejection sampling converges exponentially fast as

[TABLE]

where the spectral radius $\tau_{1}=\eta_{1}=1-\frac{1}{\omega_{0}}$.

Proof.

Let $A_{t}$ denote the number of acceptances during consecutive $t$ moves. Then

[TABLE]

[TABLE]

[TABLE]

where $P(A_{t}=0)=(1-1/\omega_{0})^{t}$, $P(A_{t}>0)=1-(1-1/\omega_{0})^{t}$, and $P(\mathbf{x},\cdot|A_{t}>0)$ has converged to $D_{\Lambda,\sigma,\mathbf{c}}$ after the first acceptance. ∎

According to Theorem 3, the convergence rate of rejection sampling depends on the choice of the normalizing constant $\omega_{0}$. Because $\omega_{0}\geq\pi(\mathbf{x})/q(\mathbf{x})$ for all $\mathbf{x}\in\mathbb{Z}^{n}$, the spectral radius $\tau_{1}=\eta_{1}$ of rejection sampling achieves the minimum when $\omega_{0}=\omega_{\text{max}}(\mathbf{x})=\omega(\mathbf{x}_{1})$, namely,

[TABLE]

thus leading to

[TABLE]

From (23) and (59), it is worth noting that only when $\omega_{0}=\omega(\mathbf{x}_{1})$, rejection sampling and independent MH have the same convergence rate. However, the former requires the knowledge of $\omega_{0}$ while the latter does not.

Remark 1**.**

Another algorithm for lattice Gaussian sampling based on rejection sampling was proposed in [40]. However, it was only concerned with values of $\sigma$ required by Klein’s algorithm. Its goal is to use rejection sampling to produce exact Gaussian samples, since Klein’s algorithm only approximates the target distribution. In contrast, our goal is to sample with smaller values of $\sigma$. The algorithm of [40] computes a certain normalizing constant in polynomial time and needs just a few steps on average to produce an exact sample. It is possible to extend their algorithm to smaller values of $\sigma$, but its running time will blow up.

IV Complexity of BDD

In this section, we apply the independent MHK sampling to BDD and analyze its complexity. The analysis for independent MH-Peikert and rejection sampling is similar, by changing the value $\delta$. As mentioned before, the decoding complexity of MCMC is evaluated by the number of Markov moves.

In MCMC, samples from the stationary distribution tend to be correlated with each other. Thus one leaves a gap, which is the mixing time $t_{\text{mix}}$, to pick up the desired independent samples (alternatively, one can run multiple Markov chains in parallel to guarantee i.i.d. samples). Therefore, we define the complexity of solving BDD by MCMC as follows.

Definition 2**.**

Let $d(\Lambda,\mathbf{c})=\min_{\mathbf{x}\in\mathbb{Z}^{n}}\|\mathbf{Bx}-\mathbf{c}\|$ denote the Euclidean distance between the query point $\mathbf{c}$ and the lattice $\Lambda$ with basis $\mathbf{B}$, and let $\widehat{\mathbf{x}}$ be the lattice point achieving $d(\Lambda,\mathbf{c})$. The complexity (i.e., the number of Markov moves $t$) of solving BDD by MCMC is

[TABLE]

Then, $C_{\mathrm{BDD}}$ can be upper bounded by

[TABLE]

where

[TABLE]

Theorem 4**.**

The complexity of solving BDD by the independent MHK algorithm is bounded above as

[TABLE]

Proof.

To start with, let us examine the numerator in (62)

[TABLE]

where we apply the Jacobi theta function $\vartheta_{3}$ [41].

By substituting (64) to (62), the complexity $C$ is upper bounded as

[TABLE]

Now, let us recall some facts about Jacobi theta function $\vartheta_{3}(\tau)$. $\vartheta_{3}(\tau)$ is monotonically decreasing with $\tau$, and particularly

[TABLE]

By simple calculation, we can get that

[TABLE]

where $\Gamma(\cdot)$ stands for the Gamma function. Clearly, if

[TABLE]

it turns out that the following term

[TABLE]

is rather small even for values of $n$ up to hundreds (e.g., $1.0039^{100}=1.4467$). The key point here is that the pre-exponential factor is close to 1. For better accuracy, $\vartheta_{3}(3)=1.00037$ (or $\vartheta_{3}(4)$ etc.) can be applied so that $1.00037^{1000}=1.4476$. More options about $\vartheta_{3}$ can be found in Table I.

Therefore, if $\sigma$ satisfies the condition (68), namely

[TABLE]

then we have

[TABLE]

Setting $\sigma=\min_{i}\|\widehat{\mathbf{b}}_{i}\|/(2\sqrt{\pi})$, we finally arrive at the following result

[TABLE]

completing the proof. ∎

Let us highlight the significance of lattice reduction. Lattice reduction is able to significantly improve $\min_{i}\|\widehat{\mathbf{b}}_{i}\|$ while reducing $\max_{i}\|\widehat{\mathbf{b}}_{i}\|$ [42]. Therefore, increasing $\min_{i}\|\widehat{\mathbf{b}}_{i}\|$ will significantly decrease the complexity shown above.

Remark 2**.**

In fact, such an analysis also holds for Klein’s algorithm, where the probability of sampling $\mathbf{x}$ follows a Gaussian-like distribution[15]

[TABLE]

Klein chose $\sigma=\min_{i}\|\mathbf{\widehat{b}}_{i}\|/\sqrt{2\log n}$, which corresponds to $O(n^{d^{2}(\Lambda,\mathbf{c})/\min_{i}\|\mathbf{\widehat{b}}_{i}\|^{2}})$ complexity. Here, we have shown that the decoding complexity can be further reduced to $O(e^{d^{2}(\Lambda,\mathbf{c})/\min_{i}^{2}\|\widehat{\mathbf{b}}_{i}\|})$, by setting $\sigma=\min_{i}\|\widehat{\mathbf{b}}_{i}\|/(2\sqrt{\pi})$. With the help of HKZ reduction, $\min_{i}\|\widehat{\mathbf{b}}_{i}\|\geq\frac{1}{n}\lambda_{1}(\Lambda)$ [43]. Thus, Klein’s algorithm allows to solve the $\eta$-BDD with $\eta=O(1/n)$ in polynomial time, while our result shown in (72) improves it to $\eta=O(\sqrt{\log n}/n)$.

According to (63), we have

[TABLE]

where $a=\log\left(\frac{1}{\epsilon}\right)\cdot 1.0039^{n}\approx\log\left(\frac{1}{\epsilon}\right)$. Clearly, the decoding radius increases with $C_{\mathrm{BDD}}$, implying a flexible trade-off between the decoding performance and complexity. In addition, the significance of lattice reduction can be seen due to an increased value of $\min_{i}\|\widehat{\mathbf{b}}_{i}\|$.

V Trapdoor Sampling

The core technique underlying GPV’s signature scheme is discrete Gaussian sampling over a trapdoor lattice [11]. Its security crucially relies on the property that the output distribution of discrete Gaussian sampling is oblivious to any particular basis used in the sampling process, therefore preventing leakage of the private key. The original GPV signature scheme was based on Klein’s algorithm, which was subsequently extended to Peikert’s algorithm [16] (see also [44, Chap. 6] for sampling over structured lattices). In fact, any good Gaussian sampling algorithms can be applied to GPV signatures. In this Section, we demonstrate the security advantage of MCMC in GPV signatures, thanks to smaller parameters it can reach.

Firstly, we provide a high-level introduction to the GPV signature (see [11, 16] for details). In key generation, one generates a hard public basis for a random lattice $\Lambda$, together with a short private basis of $\Lambda$. The public basis serves as the public key, while the private basis serves as the private key. Given a message $\mathbf{m}$ (or rather a digest of $\mathbf{m}$), one uses the private basis to sample a point $\mathbf{x}$ from $D_{\Lambda+\mathbf{m},\sigma}$ with parameter $\sigma$. The signature of $\mathbf{m}$ is $\mathbf{x}$. The verifier checks that $\mathbf{x}$ is short and that $\mathbf{x}-\mathbf{m}\in\Lambda$ using the public basis.

It is shown in [11] that the security of GPV signing can be reduced to the hardness of the inhomogeneous short integer solution (ISIS) problem333In the language of coding theory, this is to find a short vector in a coset of a linear code. with approximation factor $\sqrt{n}\sigma$. Therefore, the width $\sigma$ is the most important property of a discrete Gaussian sampler in this context.

Obviously, there is a tradeoff between security and running time in trapdoor sampling with MCMC. A small parameter $\sigma$ gives higher security, but require longer running time. Next, we examine the impact of decreasing $\sigma$ on the mixing time. Again, we focus on the independent MHK algorithm. Recall its the mixing time is proportional to

[TABLE]

Our intuition here is that if a good basis is available (as in the case of trapdoor sampling), then $\frac{1}{\delta}$ will not blow up as $n$ grows. To give an impression, Fig. 3 shows $\frac{1}{\delta}$ as a function of $n$ for checkerboard lattice $D_{n}$ with $\mathbf{c}=\mathbf{0}$ and $\sigma^{2}=-8$ dB, using its well-known basis [41, p.117, (86)]. It is seen that $\frac{1}{\delta}$ merely grows to 12 for $n$ up to 1000.

What if $\mathbf{c}\neq\mathbf{0}$? Then the denominator of (75) can be unpredictable in general. Fortunately, it can be bounded if $\sigma$ is above the smoothing parameter. Recall that for a lattice $\Lambda$ and for $\varepsilon>0$, the smoothing parameter444Note again the difference from the definition in [8], where $\sigma$ is scaled by a constant factor $\sqrt{2\pi}$. $\eta_{\varepsilon}(\Lambda)$ is defined as the smallest $\sigma>0$ such that $\sum_{{\bf{x}^{*}}\in\Lambda^{*}\setminus\{\mathbf{0}\}}e^{-2\pi^{2}\sigma^{2}\|{\bf{x}^{*}}\|^{2}}\leq\varepsilon$. If $\varepsilon<1$, we have $\frac{\rho_{\sigma,\mathbf{c}}(\mathbf{\Lambda})}{(\sqrt{2\pi}\sigma)^{n}}\in\frac{1}{\mathrm{Vol}(\Lambda)}[1-\varepsilon,1+\varepsilon]$, $\forall\mathbf{c}$.

Here, we are concerned with the parameter region $\sigma\in[\eta_{\varepsilon}(\Lambda),\sqrt{\omega({\log}\ n)}\cdot{\max}\|\mathbf{\widehat{b}}_{i}\|]$, below GPV’s parameter [11] but above the smoothing parameter. This is because we anticipate only moderate growth in mixing time but significant increase of security for values of $\sigma$ just below GPV’s parameter.

Let $I$ denote the subset of indexes $i$ with $\sqrt{2\pi}\sigma_{i}>1$ (i.e., $\sqrt{2\pi}\sigma>\|\widehat{\mathbf{b}}_{i}\|$), $i\in\{1,2,\ldots,n\}$, $|I|=m$. It is not difficult to derive the following bound, similarly to [18, Proposition 4]:

[TABLE]

where we use the identity $\vartheta_{3}\left(\frac{1}{\tau^{2}}\right)=\tau\vartheta_{3}(\tau^{2})$ and assume $\varepsilon<1/2$ in the second step, and $\vartheta_{3}(\tau)\leq 1+\sqrt{\frac{1}{\tau}}$ in the last step.

Particularly, if $\sqrt{2\pi}\sigma\geq\sqrt{\alpha}\max_{1\leq i\leq n}\|\widehat{\mathbf{b}}_{i}\|$ for some $\alpha\geq 1$, we derive

[TABLE]

Again, our key observation is that the mixing time $\vartheta_{3}(\alpha)^{n}$ grows rather slowly for values of $\alpha$ that are not too small. For example, when $\alpha=2$, we have $\vartheta_{3}(2)^{n}=1.0039^{1000}=49$ for $n=1000$. This means that with roughly $49$ iterations, our MCMC sampler is able to reduce the parameter from $\sqrt{\omega(\log n)}\max_{1\leq i\leq n}\|\widehat{\mathbf{b}}_{i}\|$ to $\frac{1}{\sqrt{\pi}}\max_{1\leq i\leq n}\|\widehat{\mathbf{b}}_{i}\|$. Therefore, if one is willing to use a slower signature scheme in return for higher security, MCMC offers such an option.

Example 1** (FALCON).**

FALCON [45] is a GPV signature scheme instantiated by NTRU lattices. Let $m$ be a power of two, $n=\varphi(m)$ where $\varphi(\cdot)$ is Euler’s totient function, $q\in\mathbb{N}$. The secret key consists of two polynomials $f$ and $g$ in ring $R=\mathbb{Z}[x]/(x^{n}+1)$ where $f$ is invertible. Find $G$ and $F$ such that

[TABLE]

The NTRU lattice of dimension $2n$ is generated by the private basis

[TABLE]

where $\mathcal{C}(\cdot)$ denotes an $n\times n$ nega-cyclic matrix whose first row consists of the coefficients of a polynomial. The public basis is given by

[TABLE]

where $h=g/f\mod q$. Both bases $\mathbf{B}$ and $\mathbf{A}$ generate the same lattice

[TABLE]

We consider the parameters $n=512$ and $q=12289$ in FALCON. The coefficients of polynomials $f$ and $g$ are randomly sampled from $D_{\mathbb{Z},4.05}$. For a particular instance randomly generated, we find $\max_{1\leq i\leq n}\|\widehat{\mathbf{b}}_{i}\|=127$. In Fig. 4, we show as a function of $\sigma$ the term $\prod_{i=1}^{2n}\theta_{3}(2\pi\sigma_{i}^{2})$ in (V), which characterizes the complexity $1/\delta$ above the smoothing parameter. It is seen that MCMC is able to significantly reduce the parameter $\sigma$, with quite moderate increase in complexity. Specifically, the term $\prod_{i=1}^{2n}\theta_{3}(2\pi\sigma_{i}^{2})$ merely grows to about 20, even if $\sigma$ is halved relative to $\max_{1\leq i\leq n}\|\widehat{\mathbf{b}}_{i}\|$. Recall that GPV sampling requires $\sigma=\sqrt{\omega(\log\ n)}\cdot\max_{1\leq i\leq n}\|\mathbf{\widehat{b}}_{i}\|$.

Note that it is possible for MCMC to incorporate the fast Fourier sampler [45], which would speed up the sampling process for structured lattices. The security levels of various samplers have been evaluated in [44, Chap. 6]. We leave evaluation of the concrete security of MCMC samplers to future work.

VI Multiple-Try Metropolis-Klein Algorithm

In this section, the independent multiple-try Metropolis-Klein (MTMK) algorithm is proposed to enhance the mixing. We firstly prove its validity and then show its uniform ergodicity with an improved convergence rate.

VI-A Multiple-Try Metropolis Method

Rather than directly generating the state candidate $\mathbf{y}$ from the proposal distribution $q(\mathbf{x},\mathbf{y})$, the multiple-try Metropolis (MTM) method selects $\mathbf{y}$ among a set of i.i.d. trial samples from $q(\mathbf{x},\mathbf{y})$, which significantly expands the searching region of proposals [46]. In particular, the MTM method consists of the following steps:

1)  Given the current state $\mathbf{X}_{t}=\mathbf{x}$, draw $k$ i.i.d. state candidates $\mathbf{y}_{1},\ldots,\mathbf{y}_{k}$ from the proposal distribution $q(\mathbf{x},\mathbf{y})$.

2)  Select $\mathbf{y}=\mathbf{y}_{c}$ among {$\mathbf{y}_{1},\ldots,\mathbf{y}_{k}$} with probability proportional to the weight

[TABLE]

where $\lambda(\mathbf{y},\mathbf{x})$ is a nonnegative symmetric function of $\mathbf{y}$ and $\mathbf{x}$ defined initially.

3)  Draw $k-1$ i.i.d. reference candidates $\mathbf{x}_{1},\ldots,\mathbf{x}_{k-1}$ from the proposal distribution $q(\mathbf{y},\mathbf{x})$ and let $\mathbf{x}_{k}=\mathbf{x}$.

4)  Accept $\mathbf{y}=\mathbf{y}_{c}$ as the state of $\mathbf{X}_{t+1}$, i.e., $\mathbf{X}_{t+1}=\mathbf{y}$ with probability

[TABLE]

otherwise, with probability $1-\alpha_{\text{MTM}}$, let $\mathbf{X}_{t+1}=\mathbf{X}_{t}=\mathbf{x}$.

By exploring the search region more thoroughly, an improvement of convergence can be achieved by MTM. Based on a number of trial samples generated from the proposal distribution, the Markov chain enjoys a large step-size jump within every single move without lowering the acceptance rate. It should be noticed that the $k-1$ reference samples $\mathbf{x}_{i}$’s are involved only for the validity of MTM by satisfying the detailed balance condition [46]

[TABLE]

Clearly, the efficiency of MTM relies on the number of trial samples $k$ while the traditional MH sampling is a special case with $k=1$. Similar to MH sampling, there is considerable flexibility in the choice of the proposal distribution $q(\mathbf{x},\mathbf{y})$ in MTM [47]. Actually, it is even possible to use different proposal distributions to generate trial samples without altering the ergodicity of the Markov chain [48]. Meanwhile, the nonnegative symmetric function $\lambda(\mathbf{x},\mathbf{y})$ in (78) is also flexible, where the only requirement is that $\lambda(\mathbf{x},\mathbf{y})>0$ whenever $q(\mathbf{x},\mathbf{y})>0$.

VI-B The Proposed Algorithm

With the great flexibility offered by $q(\mathbf{x},\mathbf{y})$ and $\lambda(\mathbf{x},\mathbf{y})$, we now propose the independent multiple-try Metropolis-Klein (MTMK) algorithm, which is described by the following steps:

1)  Given the current state $\mathbf{X}_{t}=\mathbf{x}$, use Klein’s algorithm to draw $k$ i.i.d. state candidates $\mathbf{y}_{1},\ldots,\mathbf{y}_{k}$ from the independent proposal distribution in (10)

[TABLE]

2)  Let $\lambda(\mathbf{x},\mathbf{y})=[q(\mathbf{x},\mathbf{y})q(\mathbf{y},\mathbf{x})]^{-1}=[q(\mathbf{y})q(\mathbf{x})]^{-1}$. Then select $\mathbf{y}=\mathbf{y}_{c}$ among {$\mathbf{y}_{1},\ldots,\mathbf{y}_{k}$} with probability proportional to the importance weight

[TABLE]

3)  Accept $\mathbf{y}=\mathbf{y}_{c}$ as the state of $\mathbf{X}_{t+1}$ with acceptance rate

[TABLE]

otherwise, with probability $1-\alpha_{\text{MTM}}$, let $\mathbf{X}_{t+1}=\mathbf{X}_{t}=\mathbf{x}$.

In the proposed algorithm, the basic formulation of MTM is modified in three aspects. First, Klein’s algorithm is applied to generate trial state candidates from the independent proposal distribution $q(\mathbf{x},\mathbf{y})=q(\mathbf{y})$. Then, by setting $\lambda(\mathbf{x},\mathbf{y})=[q(\mathbf{x},\mathbf{y})q(\mathbf{y},\mathbf{x})]^{-1}$, $\omega(\mathbf{x},\mathbf{y})$ becomes the importance weight of $\mathbf{x}$ that we have defined in (14). Finally and interestingly, thanks to the independent proposals, the generation of reference samples $\mathbf{x}_{i}$’s can be removed without changing the ergodicity of the chain.

In the case of independent proposals, because both the trial samples $\mathbf{y}_{i}$’s and the reference samples $\mathbf{x}_{i}$’s are generated independently from the identical distribution $q(\cdot)$, the generation of reference samples can be greatly simplified by just setting $\mathbf{x}_{i}=\mathbf{y}_{i}$ for $i=1,\ldots,c-1,c+1,\ldots,k$ and $\mathbf{x}_{c}=\mathbf{x}$. Actually, with the same arguments, the trial samples generated in the previous Markov moves can also be used by $\mathbf{x}_{i}$ [49].

It is well known that a Markov chain which is irreducible and aperiodic will be ergodic if the detailed balance condition is satisfied [19]. Since irreducible and aperiodic are easy to verify, we show the validity of the proposed algorithm by demonstrating the detailed balance condition.

Theorem 5**.**

Given the target lattice Gaussian distribution $D_{\Lambda,\sigma,\mathbf{c}}$, the Markov chain induced by the independent MTMK algorithm is ergodic.

Proof.

To start with, let us specify the transition probability $P(\mathbf{x},\mathbf{y})$ of the underlying Markov chain. For ease of presentation, we only consider the case of $\mathbf{x}\neq\mathbf{y}$, since the case $\mathbf{x}=\mathbf{y}$ is trivial. The transition probability $P(\mathbf{x},\mathbf{y})$ can be expressed as

[TABLE]

Here, $p(\mathbf{y}_{c}|\mathbf{x},c=i)$ represents the probability of accepting $\mathbf{y}=\mathbf{y}_{c}$ as the new state of $\mathbf{X}_{t+1}$ given the previous one $\mathbf{X}_{t}=\mathbf{x}$ when the $c$th candidate among $\mathbf{y}_{i}$’s is selected. Moreover, as $\mathbf{y}_{i}$ is exchangeable and independent, it follows that $p(\mathbf{y}_{i}|\mathbf{x},i)=p(\mathbf{y}_{j}|\mathbf{x},j)$ by symmetry, namely,

[TABLE]

In contrast to MH algorithms, the generation of the state candidate $\mathbf{y}=\mathbf{y}_{c}$ for Markov move $\mathbf{X}_{t+1}$ in MTM actually follows a distribution formed by $q(\mathbf{x},\mathbf{y})$ and $\omega(\mathbf{y},\mathbf{x})$ together [46]. More precisely, $p(\mathbf{y}_{c}|\mathbf{x},c)$ can be further expressed as (86), where the terms inside the sum correspond to $q(\mathbf{x},\mathbf{y})$, $\omega(\mathbf{y},\mathbf{x})$ and $\alpha$ respectively.

From (86), it is straightforward to verify the term $\pi(\mathbf{x})p(\mathbf{y}_{c}|\mathbf{x},c)$ is symmetric in $\mathbf{x}$ and $\mathbf{y}_{c}$, namely

[TABLE]

Then, by simple substitution, the detailed balance condition is satisfied as

[TABLE]

completing the proof. ∎

VI-C Convergence Analysis

Theorem 6**.**

Given the invariant lattice Gaussian distribution $D_{\Lambda,\sigma,\mathbf{c}}$, the Markov chain induced by the independent MTMK sampling algorithm converges exponentially fast to the stationary distribution:

[TABLE]

*with *

[TABLE]

The proof of Theorem 6 is provided in Appendix A.

From (90), it can be observed that with the increase of the trial sample size $k$, the exponential decay coefficient $\delta_{\text{MTM}}=\frac{k}{k-1+\frac{1}{\delta}}$ will approach 1. In other words, with a sufficiently large $k$, sampling from the target distribution can be realized efficiently. More importantly, the generation of $k$ trial samples at each Markov move not only allows a fully parallel implementation, but also can be carried out in a preprocessing stage, which is beneficial in practice.

Now, given $\delta_{\text{MTM}}=\frac{k}{k-1+\frac{1}{\delta}}$, the mixing time of the underlying Markov chain can be estimated. Specifically, according to (8) and (89), we obtain

[TABLE]

where we again use the bound $\text{ln}(1-\alpha)<-\alpha$ for $0<\alpha<1$ in $(g)$. Clearly, the mixing time is proportional to $\frac{1}{k\delta}$, and becomes $O(1)$ if $k\delta\to 1$. Overall, compared with the mixing time given in (24), the mixing time of the independent MTMK is significantly reduced by a factor of $k$. Since the independent MTMK inherits all the formulations of the independent MHK, we have

[TABLE]

for $\sigma=\min_{i}\|\widehat{\mathbf{b}}_{i}\|/(2\sqrt{\pi})$.

Following the afore-mentioned derivation, the decoding radius of the independent MTMK algorithm can be easily obtained as

[TABLE]

Remark 3**.**

Although the independent MTMK algorithm is able to reduce the mixing time, its complexity in each move increases due to multiple calls of trial samples. Therefore, parallel implementation or preprocessing is highly desired to ease the complexity burden.

Moreover, it is possible to have a varying $k$ at each Markov move, thereby resulting in an adaptive independent MTMK algorithm as

[TABLE]

where $\delta^{i}_{\text{MTM}}=\frac{k_{i}}{k_{i}-1+\frac{1}{\delta}}$ and $k_{i}$ denotes the size of trial samples at each Markov move [50].

VII Experiments of MIMO Detection

In this section, performance of the MCMC decoding algorithms is evaluated in MIMO detection. Specifically, we present simulation results for an $n\times n$ MIMO system with a square channel matrix. Here, the $i$th entry of the transmitted signal $\mathbf{x}$, denoted as $x_{i}$, is a modulation symbol taken independently from an $M$-QAM constellation $\mathcal{X}$ with Gray mapping. Meanwhile, we assume a flat fading environment, where the channel matrix $\mathbf{H}$ contains uncorrelated complex Gaussian fading gains with unit variance and remains constant over each frame duration. Let $E_{b}$ represents the average power per bit at the receiver, then the signal-to-noise ratio (SNR) $E_{b}/N_{0}=n/(\text{log}_{2}(M)\sigma_{w}^{2})$ where $M$ is the modulation level and $\sigma_{w}^{2}$ is the noise power. Then, we can express the system model as

[TABLE]

Typically, in the case of Gaussian noise $\mathbf{w}$ with zero mean and variance $\sigma_{w}^{2}$, it follows from (72) that

[TABLE]

as $\|\mathbf{Bx}-\mathbf{c}\|^{2}\approx n\sigma_{w}^{2}$ by the law of large numbers. Therefore, the decoding complexity $C$ decrease with the SNR. Note that the noise variance $\sigma_{w}^{2}$ is different from the standard deviation $\sigma$ of the lattice Gaussian distribution555In [27], the noise variance $\sigma_{w}^{2}$ is used as the sampling variance, but this would lead to a stalling problem at high SNRs [51]..

On the other hand, soft-output decoding for MIMO bit interleaver coded modulation (BICM) system is also possible using the samples generated by MCMC. Specifically, the sample candidates can be used to approximate the log-likelihood ratio (LLR), as in [52]. For bit $b_{i}\in\{0,1\}$, the approximated LLR is computed as

[TABLE]

where $b_{i}$ is the $i$th information bit associated with sample $\mathbf{x}$. The notation $\mathbf{x}:b_{i}=\mu$ means the set of all vectors $\mathbf{x}$ for which $\mathbf{x}:b_{i}=\mu$.

Fig. 5 shows the bit error rate (BER) of MCMC decoding in a $8\times 8$ uncoded MIMO system with 16-QAM, where all the samples generated by MCMC algorithms are taken into account for decoding. This corresponds to a lattice decoding scenario with dimension $n=16$. The performance of zero-forcing (ZF) and maximum-likelihood (ML) decoding are shown as benchmarks. For a fair comparison, sequential Gibbs sampling is applied here, which performs 1-dimensional conditional sampling of $x_{i}$ in a backward order666A forward update of $x_{i}$ in sequential Gibbs sampling is also possible., completing a full iteration [27]. This corresponds to one Markov move in the independent MHK and MTMK algorithms, which also update $n$ components of $\mathbf{x}$ in one iteration.

As expected, with $t=50$ Markov moves (i.e., iterations), independent MHK outperforms Gibbs sampling. As $\sigma$ has a vital impact on the sampling algorithms, Gibbs sampling is illustrated by tuning $\sigma$ with different values. Note that the detection performance may be affected due to the finite constellation. Furthermore, as shown in (74), under the help of LLL reduction, the decoding radius of the independent MHK sampling is significantly strengthened by a larger size of $\min_{i}\|\widehat{\mathbf{b}}_{i}\|$, thereby leading to a much better decoding performance. As a comparison, LLL reduction is applied in Gibbs sampling as a preprocessing stage to yield the high quality initial starting point. Additionally, compared to independent MHK, further decoding gain can be obtained by the independent MTMK algorithm, where cases with $k=5$ and $k=10$ trial samples are illustrated respectively.

On the other hand, in Fig. 6, the BERs of MCMC sampling detectors are evaluated against the number of Markov moves (i.e., iterations) in a $8\times 8$ uncoded MIMO system with 16-QAM. The SNR is fixed as $E_{b}/N_{0}=15$ dB. Clearly, the performances of all the MCMC detectors improve with the number of Markov moves. Meanwhile, with the same number of Markov moves, a substantial performance gain is obtained by LLL reduction. By increasing number of trial samples, better decoding performance can be obtained by the independent MTMK algorithm due to a larger decoding radius shown in (93).

VIII Conclusions

In this paper, the MCMC-based lattice Gaussian sampling was studied in full details. The spectral gap of the transition matrix of the independent MHK algorithm was derived and analyzed, which leads to a tractable exponential convergence rate of the Markov chain. A comparison with the extensions to Peikert’s algorithm and rejection sampling illustrated the advantages of independent MHK. With the tractable mixing time, the decoding complexity of BDD using MCMC was derived and a trade-off between the decoding radius and complexity was established. The potential of MCMC was further demonstrated in trapdoor sampling. After that, by exploiting the potential of trial samples, the independent MTMK algorithm was proposed to enhance the convergence. It supports parallel implementation due to the independent proposal distribution, thus making independent MTMK algorithm promising in practice.

Acknowledgment

The authors would like to thank Dr. Thomas Prest for helpful discussions.

Appendix A Proof of Theorem 6

Proof.

To begin with, let us take a careful look at the term $\text{min}\{\cdot,\cdot\}$ from (86). Here, for ease of presentation, we define

[TABLE]

and

[TABLE]

Meanwhile, because the $k$ trial samples from the proposal distribution $q(\cdot)$ are independent of each other, a set $\Xi$ is defined which contains the $k-1$ trial samples $\mathbf{y}_{j}$, $j\neq c$.

Then we can express $A$ and $B$ as

[TABLE]

and

[TABLE]

Here, $Q(\Xi)=\prod_{j=1,j\neq c}^{k}q(\mathbf{y}_{j})$ represents a probability distribution that takes all $q(\mathbf{y}_{j})$, $j\neq c$ into account as a whole. On the other hand, $F_{A}(\Xi)$ and $F_{B}(\Xi)$ stand for the functions about $\Xi$, namely,

[TABLE]

and

[TABLE]

where

[TABLE]

Now, let us focus on the term $A$, and we arrive at

[TABLE]

Here, $\mathbb{E}_{u(x)}[v(x)]$ represents the expectation of function $v(x)$ while $x$ is sampled from the distribution $u(x)$, $(h)$ comes from the Jensen’s inequality in the multi-variable case. Moreover, thanks to the $k-1$ independent samples from $q(\cdot)$, $(i)$ follows the derivations shown below,

[TABLE]

Similar to $A$, we can rewrite $B$ as

[TABLE]

Therefore, from (105) and (107), we get

[TABLE]

where $\delta_{\text{MTM}}=k/(k-1+\frac{1}{\delta})$ and

[TABLE]

for $\mathbf{x}\in\mathbb{Z}^{n}$ from (22) in Lemma 1. From (108), it is straightforward to see that all the Markov transitions have a component of size $\delta_{\text{MTM}}$ in common. Then, uniform ergodicity can be easily demonstrated through spectral gap or coupling technique, which is omitted here for simplicity.

∎