Estimating the Partition Function of 2-D Fields and the Capacity of Constrained Noiseless 2-D Channels

Hans-Andrea Loeliger, Mehdi Molkaraie

Description

62009 IEEE Information Theory WorkshopEstimating the Partition Function of 2-D Fields andthe Capacity of Constrained Noiseless 2-DChannels Using Tree-Based Gibbs SamplingHans-Andrea Loeliger and Mehdi MolkaraieETH ZurichDept. of Information Technology & Electrical Engineering8092 Zurich,¨ SwitzerlandEmail:floeliger, [e-mail]= = = =Abstract—Tree-basedGibbssampling(proposedbyHamzeandde Freitas) is used to compute a Monte-Carlo estimate of thepartition function of factor graphs with cycles. The proposedmethod can be used, in particular, to compute the capacity ofnoiseless constrained 2-D channels.= = = =I. INTRODUCTIONLetX ;X ;:::;X be finite sets, letX be the Cartesian1 2 N4productX =X X :::X , and letf be a nonnegative1 2 N= = = =function f :X!R. We are interested in computing (exactlyor approximately) the quantityX4Z = f(x) (1)x2X= = = =1(or, equivalently, logZ) for cases whereN Fig. 1. Forney-style factor graph for Example 1. The unlabeled boxes X ;:::X are “small” sets (e.g.,jXj =jXj =::: = 2), represent factors as in (4).1 N 1 2 N is large, and f has a “useful” factorization (as will be detailedvalue 1. Letf be the indicator function of this constraint, which canbelow).be factored into factors of the formNote that14 0; if x =x = 1k ‘(x ;x ) = (4)p(x) = f(x) (2) k ‘1; else,Zis a probability mass function onX . We will also need the set with one such factor for each adjacent pair (x ;x ).k ‘The ...

Informations

• Publié par: Ebbi
• Nombre de lectures: 23
• Langue: English

Extrait

2009 IEEE Information TheoryWorkshop

Estimating the Partition Function of 2-D Fields and the Capacity of Constrained Noiseless 2-D Channels Using Tree-Based Gibbs Sampling

Hans-Andrea Loeliger and Mehdi Molkaraie ETH Zurich Dept. of Information Technology & Electrical Engineering 8092Z  urich,Switzerland Email:{loeliger, molkarai}@isi.ee.ethz.ch

= = = = Abstract—Tree-based Gibbs sampling (proposed by Hamze and de Freitas) is used to compute a Monte-Carlo estimate of the partition function of factor graphs with cycles. The proposed method can be used, in particular, to compute the capacity of noiseless constrained 2-D channels. = = = = I. I NTRODUCTION LetX 1 ,X 2 , . . . ,X N be finite sets, letXbe the Cartesian 4 productX=X 1 × X 2 ×. . .× X N , and letfbe a nonnegative = = = = functionf:X →R. We are interested in computing (exactly or approximately) the quantity X 4 Z=f(x)(1) x∈X = = = = 1 (or, equivalently,logZ) for cases where N Fig. 1.Forney-style factor graph for Example 1. The unlabeled boxes X 1 , . . .X N are “small” sets (e.g.,|X 1 |=|X 2 |=. . .= 2), represent factors as in (4). • • Nis large, • andfhas a “useful” factorization (as will be detailed below). value1. Letfbe the indicator function of this constraint, which can be factored into factors of the form Note that  4 1 0,ifx k =x ` = 1 p(x) =f(x)(2) κ(x k , x ` ) =(4) 1,else, Z is a probability mass function onX. We will also need the set with one such factor for each adjacent pair(x k , x ` ). The corresponding Forney-style factor graph offis shown in 4 + Fig. 1, where the boxes labeled “=” are equality constraints [5]. X f ={x∈ X:f(x)6= 0}.(3) (Fig. 1 may also be viewed as a factor graph as in [4] where the boxes labeled “=” are the variable nodes.) The quantity (1) is known as the “partition function” in This example is known as the 2-D(1,∞)constrained channel [6]. statistical physics (where it is considered as a function of a “temperature” parameter that is of no concern to us here). The Note that, in this example,Z=|X f |. For this particular + 1 example 2 log computation of (1) is also the key to computing information ,lim M→∞ M 2 (Z)≈0.5879is known to nine rates of source/ channelmodels with memory [1]–[3]. decimal digits [7], [8]. However, the method proposed in this Iffhas a cycle-free factor graph with not too many states, paper works also for various generalizations of this example then the sum (1) can be computed by sum-product message for which this limit is not known to any useful accuracy [6]. passing [1], [4]. In this paper, however, we consider the case A number of Monte-Carlo methods to estimateZhave where no such cycle-free factor graph exists. In particular, we been proposed, see [9], [10]. However, these methods assume are interested in examples of the following type. thatfis strictly positive, which excludes applications as in Example 1; more about this will be said in Section II. Example 1 (Simple 2-D Constraint).Consider a grid ofN= In this paper, we propose a new Monte-Carlo method for M×Mbinary (i.e.,{0,1}-valued) variables with the constraint that no two (horizontally or vertically) adjacent variables have both the the computation ofZthat works also for Example 1. In

9781424449835/09/$25.00 © 2009 IEEE228

✬ ✩ ✬✩ constrast to the method of [3], the proposed method converges 1 to the exact value oflog(Z)in the limit of infinitely many N = = = = samples. II. E STIMATING 1/ZU SING G IBBS S AMPLING One method to estimate1/Z(and thusZitself) goes as follows. = = = = (1) (2)(K) 1) Drawsamplesx,x, .. . ,xfromX f according + top(x)defined in (2). 2) Compute K X 4 1 1 = = = = Γ =(5) (k) | K|X f f(x) + k=1 It is easily verified thatE[Γ] = 1/Z. This method was proposed in [11], see also [9]. = = = = However, there are two major issues with this method. First, it is usually assumed (as in [9], [11]) thatfis strictly positive. In this case,X f =Xand|X f |=|X |is + + ✫ ✪ ✫✪ known. However, this assumption excludes applications as in (k) Example 1. (Indeed, in Example 1, we would havef(x) = 1 Fig. 2.Partition of Fig. 1 into two cycle-free parts (one part inside the two (k) ovals, the other part outside the ovals). for all samplesx, and|X f |=Zis the desired unknown + quantity.) We will see how this issue is resolved by an idea from [12]. (1) B. Tree-BasedEstimation of1/Z[12] Second, there is the problem of generating the samplesx, (2) (K) x, . . . ,xaccording top(x). A standard general method is Let X 4 Gibbs sampling [9], [13], which, however, produces strongly f A (x A ) =f(x A , x B )(8) dependent samples. In consequence, the required number of x B samplesKis likely to exceed the limits of practicality. We and will see how this issue is eased by tree-based Gibbs sampling X 4 f B (x B ) =f(x A , x B ).(9) as proposed by Hamze and de Freitas [14]. x A III. P ROPOSED N EW M ETHOD Since The proposed method combines tree-based Gibbs sampling X X X f(x A ) =f(x B ) =f(x) =Z,(10) from [14] with an idea from [12]. x A x B x Let(A, B)be a partition of the index set{1, . . . , N}such that, (i) for fixedx A , the factor graph off(x) =f(x A , x B ) we can estimateZby applying the algorithm of Section II to is a tree and (ii) for fixedx B , the factor graph off(x) = f A or tof B (as noted in [12]). Specifically, an estimateΓ A f(x A , x B )is also a tree. An example of such a partition is of1/Zis formed as follows: shown in Fig. 2. ((1) (2)K) 1) Drawsamplesx,x. . ,, .xfrom(X A )ac- + A AA f A 4 P A. Tree-BasedGibbs Sampling [14] cording top(x A ) =p(x A , x B ) =f A (x A )/Z. x B (0) (0) (0) 2) Starting from some initial configurationx= (xx ,), A B K (k) (k) X (k) the samplesx= (x ,x),k= 1,2, . . ., are created as 4 1 1 A B Γ A =(11) (k) (k) follows. First,xis sampled according to K|(X A )| + A f f(x) A k=1 A (k−1) (k−1) p(x A |x B =x)∝f(x A , x);(6) where B B 4 (k) (X A ) ={x A :f A (x A )6= 0}.(12) + f thenxis sampled according to A B (k) (k) By symmetry, we also have an analogous estimateΓ B . The p(x B |x A =x)∝f(x ,x B ).(7) A A computation of X The point is that the sampling can be done very efficiently (k) (k) f(x) =f(x ,x B ),(13) A A in both cases since the corresponding factor graphs are cycle- x B free; see the appendix for details. Tree-based Gibbs sampling mixes much faster than naivewhich is required in (11), is easy since the corresponding Gibbs sampling [14].factor graph is a tree.

229

Fig. 3.Estimated capacity (in bits per symbol) vs. the number of samplesK for a10×10grid with a(1,∞)constraint (Example 1). The plot shows 10 independent sample paths, each with two estimates, one fromΓ A and one fromΓ B .

The quantity|(X A )|in (12) may be easy to determine + f A even iffis not strictly positive. This applies, in particular, to Example 1 (and many similar examples) where (X A ) ={x A :f(x A ,0)6= 0}.(14) + f A P In this case,|(X A )|=f(x A ,0)is easily computed by + f x A A sum-product message passing in the (cycle-free) factor graph off(x A ,0). C. AHappy Combination (1) (2) It is now obvious to create the required samplesx,x, A A (K) . . . ,xin (11) by means of tree-based Gibbs sampling as A in Section III-A. The marginals (13) may then be obtained as a by-product of the tree-based sampling (see the appendix). We thus obtain two estimates,Γ A andΓ B , as a by- product of tree-based Gibbs sampling with virtually no extra computations. IV. N UMERICAL E XPERIMENTS Some experimental results with the proposed method are shown in Figures 3 through 6. All figures refer tofas in 1 Example 1 and show the quantity (the “capacity”)log (Z). 2 N Figures 3 and 4 use a factor graph partition as in Fig. 2. In Fig. 3, we haveN= 10×10and the estimated capacity is about0.6082. In Fig. 4, we haveN= 60×60; for this size of grid there are issues with slow convergence. To improve the convergence and to speed up the mixing, we can partition the factor graph (the extension of Fig. 1 toN= 60×60) into “thicker” vertical strips. Such thick strips have cycles, but exact sum-product computation is still possible, e.g., by converting the strip into an equivalent cycle-free factor graph. The computation time is exponential in the thickness of the strip, but the faster mixing (as shown in Figures 5 and 6) results in a substantial reduction of total computation time for strips of moderate width. From Fig. 6, the estimated capacity is about0.5914.

Fig. 4.Estimated capacity (in bits per symbol) vs. the number of samplesK for a60×60grid with a(1,∞)constraint (Example 1). The plot shows 10 independent sample paths, each with two estimates, one fromΓ A and one fromΓ B .

230

Fig. 5.Same conditions as in Fig. 4, but with strips of width two.

Fig. 6.Same conditions as in Fig. 4, but with strips of width three.

Also shown in the figures (as a horizontal dotted line) is g k−1 g k g k+1 1 the infinite-grid limitlim M→∞ 2 log 2 (Z)≈0.5879, which X k−2 X k−1 X k X k+1 M ✛ ✛ is known for this simple example (see Section I). All figures show the estimates fromΓ A and fromΓ B for ←− Fig. 7.Forney-style factor graph of (16) with messagesµ X (17). k several independent experiments. ←− V. B OUNDS FOR I NFINITE G RID in Fig. 7 (cf. [5]). We then haveµ X n (x n ) = 1and X 4 4 ←−←− 1 LetC M = 2 log (Z)be the capacity of a constraint as µ(x) =g(xx ,)µ(x)(17) M 2 X k k k+1k k+1X k+1 k+1 in Example 1 for anM×Mgrid. It is clear (from tiling the x k+1 n whole plane withM×Msquares) thatC ∞ ≤C M for any X Y finiteM. =g m (x m−1 , x m )(18) x k+1 ,...,x n m=k+1 On the other hand, by tiling the plane withM×Msquares separated by all-zero guard rows and all-zero guard columns, fork=n−1, n−2, . . . ,1. Then M 2 we obtainC ∞ ≥C M ( ). X M+1 p(x 1 ) =p(x 1 , . . . , x n )(19) In the example of Figures 4–6 (withM= 60), we thus x 2 ,...,x n obtain0.5721≤C ∞ ≤0.5914. ←− ∝µ X 1 (x 1 )(20) VI. C ONCLUDING R EMARKS and ←− g k (x k−1 , x k )µ X k (x k ) p(x k |x k−1 ) =(21) ←− We have shown that tree-based Gibbs sampling (as proposed µ X k−1 (x k−1 ) by Hamze and de Freitas) can be used to compute an estimate fork= 2, . . . , n. The proof of (21) follows from noting that of the partition function with virtually no extra computational −→←− cost. The proposed method can be used, in particular, to p(x) =γ µ(x)µ(x)(22) k−1X k−1 k−1X k−1 k−1 compute (a Monte Carlo estimate of) the capacity of noiseless and constrained 2-D channels. Our preliminary numerical experi- −→←− ments are encouraging. p(x k−1 , x k ) =γ µ X k−1 (x k−1 )g k (x k−1 , x k )µ X k (x k )(23) −→ whereµ X k−1 is the forward sum-product message along the A PPENDIX : S AMPLING FROM M ARKOV C HAINS edgeX k−1 and whereγis the missing scale factor in (16). We recall some pertinent facts about the simulation of We also note that X X Markov chains and cycle-free factor graphs. Letp(x) = ←− µ X 1 (x 1 ) =g(x)(24) p(x 1 , . . . , x n )be the probability mass function of a Markov x 1 x chain. Ifp(x)is given in the form whereg(x)is defined as the right-hand side of (16). In this n Y paper, this fact is used to compute the marginals (13) as a p(x) =p(x 1 )p(x k |x k−1 ),(15) by-product of the sampling. k=2 The generalization of all this to arbitrary factor graphs without cycles is straightforward. then it is obvious how to create i.i.d. samples according to p(x). Now consider the case wherep(x)is not given in the A CKNOWLEDGEMENT form (15), but in the more general form The authors wish to thank David MacKay and Iain Murray n Y for pointing out to us the work by Hamze and de Freitas [14]. p(x)∝g k (x k−1 , x k )(16) k=2 R EFERENCES with general factorsg k . It is then still easy to create i.i.d. samples according top(x), which may be seen as follows. , ˇ vicc  a,dn.WeZgn.O.Vontobel,A.Ka,PerigelLoA..-,HdlonrA.D]1[ “Simulation-based computation of information rates for channels with First, a probability mass function of the form (16) can be memory,”IEEE Trans. Inform. Theory, vol. 52, no. 8, pp. 3498–3508, rewritten in the form (15) (which allows efficient simulation). August 2006. Second, this reparameterization ofp(x)may be efficiently [2] P. H. Siegel, “Information-theoretic limits of two-dimensional optical recording channels,” inOptical Data Storage 2006(Proc. of SPIE, Vol. carried out by backward sum-product message passing, as 6282, Eds. Ryuichi Katayama and Tuviah E. Schlesinger), Montreal, will be detailed below. The resulting algorithm is know as Quebec, Canada, April 23–26, 2006, pp. 62820W-1 – 2820W-13. “backward-filtering forward-sampling” (or, in a time-reversed [3] O. Shental, N. Shental, S. Shamai (Shitz), I. Kanter, A. J. Weiss, and Y. Weiss, “Discrete-input two-dimensional Gaussian channels with version, as “forward-filtering backward-sampling”) [15]. memory: estimation and information rates via graphical models and ←− Specifically, letµ X k be the backward sum-product message statistical mechanics,”IEEE Trans. Inform. Theory,vol. 54, pp. 1500– along the edgeX k in the factor graph of (16), as is illustrated 1513, April 2008.

231

[4] F.R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,”IEEE Trans. Inform. Theory,vol. 47, pp. 498– 519, Feb. 2001. [5] H.-A.Loeliger, “An introduction to factor graphs,”IEEE Signal Proc. Mag.,Jan. 2004, pp. 28–41. [6] K.Kato and K. Zeger, “On the capacity of two-dimensional run-length constrained channels,”IEEE Trans. Inform. Theory,vol. 45, pp. 1527– 1540, July 1999. [7] N.J. Calkin and H. S. Wilf, “The number of independent sets in a grid graph,”SIAM J. Discr. Math.,vol. 11, pp. 54–60, Feb. 1998. [8] W.Weeks and R. E. Blahut, “The capacity and coding gain of certain checkerboard codes,”IEEE Trans. Inform. Theory,vol. 44, pp. 1193– 1203, May 1998. [9] G.Potamianos and J. Goutsias, “Stochastic approximation algorithms for partition function estimation of Gibbs random fields,”IEEE Trans. Inform. Theory,vol. 43, pp. 1984–1965, Nov. 1997.

[10] F.Huang and Y. Ogata, “Comparison of two methods for calculating the partition functions of various spatial statistical models,”Australian & New Zealand Journal of Statistics,vol. 43 (1), pp. 47–65, 2001. [11] Y. Ogata and M. Tanemura, “Estimation of interaction potentials of spatial point patterns through the maximum likelihood procedure,”Ann. Inst. Statist. Math.,vol. 33, pp. 315–338, 1981. [12] H.-A.Loeliger and M. Molkaraie, “Simulation-based estimation of the partition function and the information rate of two-dimensional models,” Proc. 2008 IEEE Int. Symp. on Information Theory,Toronto, Canada, July 6–11, 2008, pp. 1113–1117. [13] D.J. C. MacKay, “Introduction to Monte Carlo methods,” inLearning in Graphical Models,M. I. Jordan, ed., Kluwer Academic Press, 1998, pp. 175–204. [14] F.Hamze and N. de Freitas, “From fields to trees,”Proc. Conf. on Uncertainty in Artificial Intelligence,Banff, July 2004. [15] D.Gamerman and H. F. Lopes,Markov Chain Monte Carlo: Stochastic Simulation for Bayesian Inference.2nd ed., CRC Press, 2006.

232