Sampling from a lattice Gaussian distribution is emerging as an important problem in coding and cryptography. This paper gives a further analysis of the independent Metropolis-Hastings-Klein (MHK) algorithm we presented at ISIT 2015. We derive the exact spectral gap of the induced Markov chain, which dictates the convergence rate of the independent MHK algorithm. Then, we apply the independent MHK algorithm to lattice decoding and obtained the decoding complexity for solving the CVP as Õ(e∥Bx-c∥2 / mini ∥b̂i∥2). Finally, the tradeoff between decoding radius and complexity is also established

Ling, C

Wang, Z

English

Zheng Wang

Cong Ling

Crossref

Further results on independent Metropolis-Hastings-Klein sampling

Spiral - Imperial College Digital Repository

Further Results on IndependentMetropolis-Hastings-Klein SamplingZheng Wang and Cong LingDepartment of EEEImperial College LondonLondon, SW7 2AZ, United KingdomEmail: z.wang10, c.ling@imperial.ac.ukAbstract—Sampling from a lattice Gaussian distribution isemerging as an important problem in coding and cryptography.This paper gives a further analysis of the independent Metropolis-Hastings-Klein (MHK) algorithm we presented at ISIT 2015. Wederive the exact spectral gap of the induced Markov chain, whichdictates the convergence rate of the independent MHK algorithm.Then, we apply the independent MHK algorithm to the closestvector problem (CVP) in lattice decoding, and establish the trade-off between decoding radius and complexity.Index Terms—Lattice Gaussian sampling, MCMC methods,Metropolis-Hastings algorithm, closest vector problem.I. INTRODUCTIONAs a core problem of lattice theory, the closest vector prob-lem (CVP) has drawn a lot of attention. The CVP is normallysolved by sphere decoding (SD) or its variants. However, SDis inapplicable in high dimensions due to the exponentiallyincreased complexity. To this end, Klein introduced a samplingalgorithm to solve the CVP [1], which performs the decodingby sampling over a Gaussian-like distribution. Klein showedits decoding complexity is O(n‖Bx−c‖2/mini ‖b̂i‖2), where‖Bx−c‖ denotes the Euclidean distance from the query pointc to the lattice and b̂i’s are the Gram-Schmidt vectors of thelattice basis B. In [2], [3], improved sampling algorithms wereproposed for near-optimal decoding performance. Neverthe-less, these algorithms only yield an approximately Gaussiandistribution.As a basic scheme in Markov chain Monte Carlo (MCMC),Gibbs sampling has been introduced to lattice decoding, whichemploys univariate conditional sampling to build a Markovchain [4]–[9]. As the Markov chain converges to the stationarydistribution, the optimal solution of the CVP can be encoun-tered by sampling. However, to the best of our knowledge, theanalysis of the convergence rate for Gibbs sampling in theseapplications is still lacking.In [10], we proposed independent Metropolis-Hastings-Klein (MHK) sampling and derived its rate of convergenceusing the conventional coupling technique. In this paper, wefurther study the independent MHK algorithm and examine itscomplexity in solving the CVP. We derive the exact spectralgap of the transition matrix, thus precisely determine theconvergence rate of the underlying Markov chain. Moreover,we show its decoding complexity for solving the CVP isÕ(e‖Bx−c‖2/mini ‖b̂i‖2). We demonstrate that the independentMHK sampling decoder enjoys a flexible trade-off betweendecoding radius and complexly.The rest of this paper is organized as follows. Section IIgives a briefly review of the independent MHK sampling algo-rithm. In Section III, we derive the spectral gap of the Markovchain associated with the independent MHK algorithm. InSection IV, we examine the decoding complexity and establishthe flexible decoding trade-off. Finally, Section V concludesthe paper.II. INDEPENDENT MHK ALGORITHMLet B = [b1, . . . ,bn] ⊂ Rn consist of n linearly indepen-dent vectors. The n-dimensional lattice Λ generated by B isdefined byΛ = {Bx : x ∈ Zn}. (1)We define the Gaussian function centered at c for standarddeviation σ > 0 asρσ,c(z) = e− ‖z−c‖22σ2 , (2)for all z ∈ Rn. When c or σ are not specified, we assume thatthey are 0 and 1 respectively. Then, the discrete Gaussiandistribution over Λ is defined asDΛ,σ,c(x) =ρσ,c(Bx)ρσ,c(Λ)=e−12σ2‖Bx−c‖2∑x∈Zn e− 12σ2‖Bx−c‖2(3)for all Bx ∈ Λ, where ρσ,c(Λ) ,∑Bx∈Λ ρσ,c(Bx) is just ascaling to yield a probability distribution. Intuitively, the latticeGaussian distribution can be used to sample the closest latticepoint. The closer Bx to c, the larger probability it will besampled, indicating that the optimal solution of x is the mostlikely to be sampled.From the MCMC perspective, the discrete Gaussian distri-bution can be viewed as a complex target distribution lackingof direct sampling methods. In order to obtain samples fromDΛ,σ,c(x), the independent MHK sampling was proposed in[10]. Specifically, in the independent MHK sampling, a statecandidate y for the next Markov move Xt+1 is generatedby Klein’s algorithm [1], via the following backward one-dimensional conditional sampling:P (yi|y[−i]) = P (yi|yi+1, . . . , yn)=e−12σ2‖c′−Ry‖2∑yi∈Z e− 12σ2‖c′−Ry‖2=e−12σ2‖c′i−∑nj=i ri,jyj‖2∑yi∈Z e− 12σ2‖c′i−∑nj=i ri,jyj‖2=e− 12σ2i‖yi−ỹi‖2∑yi∈Z e− 12σ2i‖yi−ỹi‖2= DZ,σi,ỹi(yi) (4)where ỹi =c′i−∑nj=i+1 ri,jyjri,i, σi = σ|ri,i| , c′ = Q†c andB = QR. Note that y[−i] = [yi+1, . . . , yn], R, c′ and y arethe (n− i+ 1) segments of R, c′ and y respectively (i.e., Ris a (n− i+ 1)× (n− i+ 1) submatrix of R with ri,i to rn,nin the diagonal).Given the current state x of the Markov chain Xt, the pro-posal distribution q(x,y) in the independent MHK samplingis given byq(x,y) =n∏i=1P (yn+1−i|y[−(n+1−i)])=ρσ,c(By)∏ni=1 ρσn+1−i,ỹn+1−i(Z). (5)Actually, the proposal distribution is independent of x fromXt, that isq(x,y) = q(y), (6)implying the connection between two consecutive Markovmoves is only due to following the decision stage.More precisely, with the state candidate y, the acceptanceratio α is computed byα(x,y) = min{1,π(y)q(y,x)π(x)q(x,y)}= min{1,∏ni=1 ρσn+1−i,ỹn+1−i(Z)∏ni=1 ρσn+1−i,x̃n+1−i(Z)}, (7)where x̃i =c′i−∑nj=i+1 ri,jxjri,iand π = DΛ,σ,c (these nota-tions will be followed throughout the context). Then, y willbe accepted as the new state by Xt+1 with probability α.Otherwise, x will be retained by Xt+1. In this way, a Markovchain {X0,X1, . . .} is established, which is summarized inAlgorithm 1.Note that in MH algorithms, the proposal distributionq(x,y) can be any fixed distribution from which we canconveniently draw samples. Meanwhile, the initial state x0for X0 can be chosen from Zn arbitrarily or from the outputof a suboptimal algorithm.III. SPECTRAL GAPThe proof of the following theorem was given in [10]using the coupling technique. Here, we reformulate it in termsof the spectral gap of the Markov chain and give a morestraightforward proof following [11].Algorithm 1 Independent Metropolis-Hastings-Klein Sam-pling AlgorithmInput: B, σ, c,x0;Output: x ∼ DΛ,σ,c;1: let X0 = x02: for t =1,2, . . . , do3: let x denote the state of Xt−14: sample y from the proposal distribution q(x,y) in (5)5: calculate the acceptance ratio α(x,y) in (7)6: generate a sample u from the uniform density U [0, 1]7: if u ≤ α(x,y) then8: let Xt = y9: else10: Xt = x11: end if12: output x = Xt13: end forTheorem 1 (Uniform ergodicity of independent MHK). Giventhe invariant lattice Gaussian distribution DΛ,σ,c, the Markovchain induced by independent MHK sampling converges expo-nentially fast to the stationary distribution in total variationaldistance:‖P t(x, ·)−DΛ,σ,c(·)‖TV ≤ (1− δ)t, (8)where δ is a lower bound on the spectral gap to be given inLemma 1.Proof. From (5) and (7), the transition probability P (x,y) ofeach Markov move in the independent MHK sampling is givenbyP (x,y) =min{q(y), π(y)q(x)π(x)}if y 6= x,q(x)+∑z6=xmax{0,q(z)− π(z)q(x)π(x)}if y = x.(9)For notational simplicity, here we define the importanceweight w(x) asw(x) =π(x)q(x). (10)Then the transition probability can be rewritten asP (x,y)=min{q(y), π(y)w(x)}if y 6= x,q(x)+∑z6=xmax{0,q(z)− π(z)w(x)}if y = x.(11)Without loss of generality, we label the countably infinitestate space Ω = Zn as Ω = {x1,x2, . . . ,x∞} (i.e., y = xi ∈Ω ) and assume that these states are sorted according to themagnitudes of their importance weights, namely,w(x1) ≥ w(x2) ≥ · · · ≥ w(x∞). (12)From (11) and (12), the transition matrix P of the Markovchain can be exactly expressed asP =q(x1) + λ1π(x2)w(x1)π(x3)w(x1)· · · π(x∞)w(x1)q(x1) q(x2) + λ2π(x3)w(x2)· · · π(x∞)w(x2)q(x1) q(x2) q(x3) + λ3 · · · π(x∞)w(x3).......... . ....q(x1) q(x2) q(x3) · · · q(x∞)whereλj =∞∑i=j(q(xi)−π(xi)w(xj))(13)stands for the probability of being rejected in the decisionstage with the current state xj for Xt.Let q = [q(x1), q(x2), . . .]T denote the vector of proposalprobabilities. Then by decomposition, it follows thatP = G + eqT , (14)where e = [1, 1, . . .]T and G is an upper triangular matrix ofthe formG =λ1π(x2)w(x1)− q(x2) · · · π(x∞)w(x1) − q(x∞)0 λ2 · · · π(x∞)w(x2) − q(x∞)....... . ....0 0 · · · 0 .It is well-known that for a Markov chain, the largesteigenvalue of the transition matrix P always equals 1. Here,as e is a common right eigenvector for both P and P−G, itnaturally corresponds to the largest eigenvalue 1. Meanwhile,since the rank of P−G is 1, the other eigenvalues of G areexactly the same as those of P.On the other hand, the mixing time, which measures thetime required by a Markov chain to get close to its stationarydistributiontmix(ε) = min{t : max‖P t(x, ·)− π(·)‖TV ≤ ε}, (15)is asymptotically dominated by the second largest eigenvalueλ∗ [12], that istmix(ε) ≤11− |λ∗|log(1ε). (16)Thanks to the ascending order in (12), it is easy to verify thatλ∗ = λ1 (17)and1 > |λ1| ≥ |λ2| ≥ · · · > 0, (18)thereby raising the interest of identifying the value of λ1.Therefore, according to (13), we can easily get thatλ1 =∞∑i=1(q(xi)−π(xi)w(x1))=∞∑i=1q(xi)−1w(x1)·∞∑i=1π(xi)= 1− 1w(x1)= 1− q(x1)π(x1). (19)In other words, the spectral gap 1 − λ1 is exactly capturedby the ratio q(x1)/π(x1).Next, we invoke the following Lemma to lower bound theratio q(x)/π(x).Lemma 1 ([10]). In the independent MHK algorithm, thereexists δ > 0 such thatq(x)π(x)≥ δ (20)for all x ∈ Ω, whereδ ,ρσ,c(Λ)∏ni=1 ρσi(Z). (21)Therefore, from Lemma 1, we can immediately getλ1 ≤ 1− δ, (22)and by substituting (22) to (16), it follows thattmix(ε) ≤1δlog(1ε). (23)Moreover, given the mixing time shown in (15), we can easilyderive that‖P t(x, ·)− π(·)‖TV ≤ (1− δ)t, (24)completing the proof.The rate of convergence derived here is the same as thatgiven in [10]. Since the convergence rate derived from thespectral gap is precise, this means that the analysis in [10]based on the coupling technique is tight. More details on theexponential decay coefficient δ can be found in [13].IV. DECODING COMPLEXITYIn this section, we apply the independent MHK samplingto lattice decoding and examine its complexity. We evaluatethe complexity with the number of Markov moves, since thecomplexity of each Markov move is often insignificant inpractice.Consider the decoding of an n×n real-valued system. Theextension to the complex-valued system is straightforward [2].Let x denote the transmitted signal ∈ Zn. The correspondingreceived signal c is given byc = Bx + w (25)where w is the noise vector with zero mean and variance σ2w.The ML decoding is given byx̂ = arg minx∈Zn‖c−Bx‖2 (26)where ‖ · ‖ denotes the Euclidean norm. Intuitively, this CVPcan be solved with the independent MHK algorithm, sincethe closest lattice point will be sampled with the highestprobability. Note that another lattice problem — the shortestvector problem (SVP), is essentially a special case of CVP byc = 0.Because samples tend to be correlated with each other,we leave a gap tmix to pick up samples from the stationarydistribution. Therefore, we define the complexity of solvingCVP as follows:Definition 1. With independent MHK sampling, the complexityof solving CVP isCMHK =tmixDΛ,σ,c(x). (27)Note that the mixing time tmix serves here as the gapfor independent samples, which is equivalent to operate theMarkov chains in parallel. In other words, it corresponds tothe decoding in the worst case for MCMC. Therefore, it maybe possible to reduce the gap to obtain lower complexity.With the mixing time given in (23), CMHK can be upperbounded byCMHK < log(1ε)·∏ni=1 ρσn+1−i,x̃n+1−i(Z)ρσ,c(Λ)· ρσ,c(Λ)ρσ,c(Bx)= a ·∏ni=1 ρσn+1−i,x̃n+1−i(Z)ρσ,c(Bx)= a · C. (28)As a = log(1ε)is a constant depending on the value of ε, wepay attention to the complexityC =∏ni=1 ρσn+1−i,x̃n+1−i(Z)ρσ,c(Bx). (29)For the term∏ni=1 ρσn+1−i,x̃n+1−i(Z), we haven∏i=1ρσn+1−i,x̃n+1−i(Z)(a)≤n∏i=1ρσi(Z)=n∏i=1∑xi∈Ze− 12σ2i‖xi‖2=n∏i=1∑xi∈Ze−α‖b̂i‖2π2πασ2‖xi‖2(b)≤n∏i=1(2πασ2)n2∑xi∈Ze−α‖b̂i‖2π‖xi‖2(c)= (√2πασ)n2·n∏i=1ϑ3(α‖b̂i‖2). (30)Here, (a) holds due to the fact that [14]ρσ,c(Λ) ≤ ρσ(Λ), (31)(c) applies the Jacobi theta function ϑ3 [15]ϑ3(τ) =+∞∑n=−∞e−πτn2, (32)and (b) follows the inequality from [16, Lemma 1.4]∑x∈Λe−πs−1x2 ≤ sn2 ·∑x∈Λe−πx2, (33)where s ≥ 1.Note the condition for (b) in (30), namely,(2πασ2) ≥ 1, (34)where the equality holds ifσ =1√2πα. (35)By substituting (30) to (29), the complexity C is upperbounded asC ≤(√2πασ)n2 ·∏ni=1 ϑ3(α‖b̂i‖2)e−12σ2‖Bx−c‖2. (36)It is natural that σ is chosen to minimize the above lowerbound. However, as can be seen clearly, once 2πασ2 > 1, thenthe power of 2πασ2 will be the unaffordable n2. Therefore,in order to remove the impact of such an item, we choose σto make2πασ2 = 1 (37)andα =12πσ2. (38)Then, we can rewrite the upper bound (36) asC ≤∏ni=1 ϑ3(‖b̂i‖2/2πσ2)e−12σ2‖Bx−c‖2. (39)Now we recall some facts about the Jacobi theta functionϑ3(τ) from (32). It is obvious that ϑ3(τ) is monotonicallydecreasing with τ , and particularly the minimum of ϑ3(τ) is1, namely,limτ→∞inf ϑ3(τ) = 1. (40)By simple calculation, we can get thatϑ3(2) =+∞∑n=−∞e−2πn2=4√6π + 4√2π2Γ( 34 )= 1.0039, (41)where Γ(·) stands for the Gamma function. Clearly, ifmin1≤i≤n ‖b̂i‖22πσ2≥ 2 (42)then it turns out thatn∏i=1ϑ3(‖b̂i‖2/2πσ2) ≤ ϑn3 (2) = 1.0039n, (43)is really small even for values of n up to hundreds (e.g.,1.0039100 = 1.4759).Therefore, let σ satisfy the condition given in (42), that isσ ≤ min1≤i≤n‖b̂i‖/(2√π), (44)then we haveC ≤ 1.0039n · e12σ2‖Bx−c‖2 . (45)Here, we can also apply ϑ3(3) = 1.0002 (or ϑ3(4), etc.)so that 1.00021000 = 1.2214. The key point is that thepre-exponential factor is rather small. Therefore, set σ =mini ‖b̂i‖/(2√π), then the decoding complexity is given byC = Õ(e‖Bx−c‖2/mini ‖b̂i‖2). (46)By law of large numbers, ‖Bx − c‖2 ≈ nσ2w. Thus C ≈O(enσ2w/mini ‖b̂i‖2). We highlight the significance of lattice re-duction (i.e., LLL reduction) here, since increasing mini ‖b̂i‖will significantly decrease the complexity.Remark 1. In fact, such an analysis also holds for Klein’salgorithm, where the probability of sampling x is [1]P (x) ≥ e− 12σ2‖Bx−c‖2∏ni=1 ϑ3(‖b̂i‖2/2πσ2), (47)which is exactly the inverse of (39). Klein chose σ =mini ‖b̂i‖/√log n and showed the decoding complexityO(n‖Bx−c‖2/mini ‖b̂i‖2).Here, we have shown that the decoding complexity canbe reduced to Õ(e‖Bx−c‖2/mini ‖b̂i‖2), by setting σ =mini ‖b̂i‖/(2√π).In general, the complexity of solving the CVP is exponentialor higher. On the other hand, one may be interested in theperformance with a fixed number of Markov moves C  1.It follows from (45) that‖Bx− c‖ ≥ σ√2 lnC1.0039n. (48)Typically, the sampling decoder will find the closest vectorpoint Bx if the distance from c to the lattice is less than theright-hand-side (RHS) of (48). In this regard, the RHS of (48)can be interpreted as the decoding radius of bounded distancedecoding (BDD), that isd = ‖Bx− c‖ , σ√2 lnC1.0039n. (49)When σ = mini ‖b̂i‖/(2√π), we haved =√12π· ln C1.0039n· min1≤i≤n‖b̂i‖. (50)Clearly, the decoding radius d monotonically increases withC, implying a flexible trade-off between the decoding radiusand complexity.V. CONCLUSIONSIn this paper, the independent MHK algorithm was fur-ther studied. The exact spectral gap of the transition matrixwas derived, which led to a precise estimate of the con-vergence rate of the underlying Markov chain. Then, weapplied the independent MHK sampling to lattice decodingand obtained the decoding complexity for solving the CVPas Õ(e‖Bx−c‖2/mini ‖b̂i‖2). To evaluate the decoding perfor-mance, the decoding radius was further derived, revealing aflexible trade-off between performance and complexity.ACKNOWLEDGMENTThis work was supported in part by FP7 project PHYLAWS(EU FP7-ICT 317562).REFERENCES[1] P. Klein, “Finding the closest lattice vector when it is unusually close,”in ACM-SIAM Symp. Discr. Algorithms, 2000, pp. 937–941.[2] S. Liu, C. Ling, and D. Stehlé, “Decoding by sampling: A randomizedlattice algorithm for bounded distance decoding,” IEEE Trans. Inform.Theory, vol. 57, pp. 5933–5945, Sep. 2011.[3] Z. Wang, S. Liu, and C. Ling, “Decoding by sampling - Part II:Derandomization and soft-output decoding,” IEEE Trans. Commun.,vol. 61, no. 11, pp. 4630–4639, Nov. 2013.[4] B. Hassibi, M. Hansen, A. Dimakis, H. Alshamary, and W. Xu,“Optimized Markov Chain Monte Carlo for signal detection in MIMOsystems: An analysis of the stationary distribution and mixing time,”IEEE Transactions on Signal Processing, vol. 62, no. 17, pp. 4436–4450, Sep. 2014.[5] T. Datta, N. Kumar, A. Chockalingam, and B. Rajan, “A novel MonteCarlo sampling based receiver for large-scale uplink multiuser MIMOsystems,” IEEE Transactions on Vehicular Technology,, vol. 62, no. 7,pp. 3019–3038, Sep. 2013.[6] B. Farhang-Boroujeny, H. Zhu, and Z. Shi, “Markov chain MonteCarlo algorithms for CDMA and MIMO communication systems,” IEEETrans. Signal Process., vol. 54, no. 5, pp. 1896–1909, 2006.[7] R. Chen, J. Liu, and X. Wang, “Convergence analysis and comparisonsof Markov chain Monte Carlo algorithms in digital communications,”IEEE Trans. on Signal Process., vol. 50, no. 2, pp. 255–270, 2002.[8] P. Aggarwal and X. Wang, “Multilevel sequential Monte Carlo al-gorithms for MIMO demodulation,” IEEE Transactions on WirelessCommunications, vol. 6, no. 2, pp. 750–758, Feb. 2007.[9] H. Zhu, B. Farhang-Boroujeny, and R.-R. Chen, “On performance ofsphere decoding and Markov chain Monte Carlo detection methods,”IEEE Signal Processing Letters, vol. 12, no. 10, pp. 669–672, 2005.[10] Z. Wang and C. Ling, “Independent Metropolis-Hastings-Klein al-gorithm for lattice Gaussian sampling,” in Proc. IEEE InternationalSymposium on Information Theory (ISIT), Hong Kong, China, June2015, pp. 2470–2474.[11] J. S. Liu, “Metropolized independent sampling with comparisons torejection sampling and importance sampling,” Statistics and Computing,vol. 6, pp. 113–119, 1996.[12] D. A. Levin, Y. Peres, and E. L. Wilmer, Markov Chains and MixingTime, American Mathematical Society, 2008.[13] Z. Wang and C. Ling, “On the geometric ergodicity of Metropolis-Hastings algorithms for lattice Gaussian sampling,” Submitted toIEEE Transactions on Information Theory, 2015., [Online] Avail-able:http://arxiv.org/pdf/1501.05757v2.pdf.[14] D. Micciancio and O. Regev, “Worst-case to average-case reductionsbased on Gaussian measures,” in Proc. Ann. Symp. Found. ComputerScience, Rome, Italy, Oct. 2004, pp. 372–381.[15] J. H. Conway and N. A. Sloane, Sphere Packings, Lattices and Groups.New York: Springer-Verlag, 1998.[16] W. Banaszczyk, “New bounds in some transference theorems in thegeometry of numbers,” Math. Ann., vol. 296, pp. 625–635, 1993.

http://spiral.imperial.ac.uk/bitstream/10044/1/40322/2/bare_jrnl.pdf

Further results on independent Metropolis-Hastings-Klein sampling

Abstract

Similar works

Full text

Available Versions

Crossref

Spiral - Imperial College Digital Repository