Generalized Belief Propagation for the Noiseless Capacity and Information Rates of Run-Length Limited Constraints

1 Generalized Belief Propagation for the Noiseless Capacity and Information Rates of Run-Length Limited Constraints Giov anni Sabato, Member , IEEE and Mehdi Molkaraie, Member , IEEE Abstract —The performance of the generalized belief propa- gation algorithm to compute the noiseless capacity and mu- tual information rates of finite-size two-dimensional and three- dimensional run-length limited constraints is inv estigated. In both cases, the problem is reduced to estimating the partition function of graphical models with cycles. The partition function is then estimated using the region-based free energy approxi- mation technique. For each constraint, a method is proposed to choose the basic r egions and to construct the region graph which provides the graphical framework to run the generalized belief propagation algorithm. Simulation r esults for the noiseless capacity of different constraints as a function of the size of the channel are reported. In the cases that tight lower and upper bounds on the Shannon capacity exist, conv ergence to the Shannon capacity is discussed. For noisy constrained channels, simulation results are reported for mutual information rates as a function of signal-to-noise ratio. Index T erms —Generalized belief propagation algorithm, run- length limited constraints, partition function, factor graphs, region graphs, noiseless capacity , Shannon capacity , mutual information rate. I . I N T R O D U C T I O N Run-length limited (RLL) constraints are widely used in magnetic and optical recording systems. Such constraints reduce the effect of inter-symbol interference and help in timing control. In track-oriented storage systems constraints are defined in one dimension. W e say a binary one-dimensional (1-D) sequence satisfies the ( d, k ) -RLL constraint if the runs of 0 ’ s have length at most k and the runs of 0 ’ s between successive 1 ’ s have length at least d . W e suppose that 0 ≤ d < k ≤ ∞ . The Shannon capacity of a 1-D ( d, k ) -RLL constraint is defined as C ( d,k ) 1 D 4 = lim m →∞ log 2 Z ( m ) m , (1) Giov anni Sabato is with P ARALLEL Informatik AG, CH-6005 Luzern, Switzerland. Mehdi Molkaraie is with the Dept. of Information T echnol- ogy and Electrical Engineering, ETH Z ¨ urich, CH-8092 Z ¨ urich, Switzerland. Emails: giovanni.sabato@parallel.ch, molkaraie@isi.ee.ethz.ch. where Z ( m ) denotes the number of binary 1-D sequences of length m that satisfy the ( d, k ) -RLL constraint, see [1], [2]. W ith the rise in demand for larger storage in smaller size and with recent developments in page-oriented storage systems, such as holographic data storage, two-dimensional (2-D) constraints have become more of interest [3]. In these systems, data is org anized on a surface and constraints are defined in two dimensions. A 2-D binary array satisfies the ( d 1 , k 1 , d 2 , k 2 ) -RLL con- straint if it satisfies a ( d 1 , k 1 ) -RLL constraint horizontally and a ( d 2 , k 2 ) -RLL constraint vertically . If a 2-D binary array satisfies a 1-D ( d, k ) -RLL constraint both horizontally and vertically , we simply say that it satisfies a 2-D ( d, k ) -RLL constraint. Example: 2-D (2 , ∞ ) -RLL constraint: The 2-D (2 , ∞ ) -RLL constraint is satisfied in the following 2-D binary array segment. In words, in every row and ev ery column of the array there are at least two 0 ’ s between succes- siv e 1 ’ s; but the runs of 0 ’ s can be of any length (howe ver , 1 ’ s can be diagonally adjacent). . . . 0100100001001000100000100010 . . . . . . 1000010000100010000100000100 . . . . . . 0001000010000001000000010001 . . . . . . 0100100100010000001000100000 . . . The Shannon capacity of a 2-D ( d 1 , k 1 , d 2 , k 2 ) -RLL con- straint is defined as C ( d 1 ,k 1 ,d 2 ,k 2 ) 2 D 4 = lim m,n →∞ log 2 Z ( m, n ) mn , (2) where Z ( m, n ) denotes the number of 2-D binary arrays of size m × n that satisfy the ( d 1 , k 1 , d 2 , k 2 ) -RLL constraint. Similarly , the Shannon capacity can be defined for higher dimensional constrained channels. For example, the Shannon capacity in three dimensions C ( d 1 ,k 1 ,d 2 ,k 2 ,d 3 ,k 3 ) 3 D depends on Z ( m, n, q ) , the number of three-dimensional (3-D) binary 2 arrays of size m × n × q that satisfy a ( d 1 , k 1 , d 2 , k 2 , d 3 , k 3 ) - RLL constraint. The noiseless capacity of a constrained channel is an impor- tant quantity that provides an upper bound to the information rate of any encoder that maps arbitrary binary input into binary data that satisfies a given constraint. There are a number of techniques to compute the 1-D Shannon capacity (for example combinatorial or algebraic approaches) [1]. In contrast to the 1-D capacity , except for a few cases, e xact v alues of two and higher dimensional (positiv e) Shannon capacities are not known, see [4]–[9]. For noisy 1-D constrained channels, simulation-based tech- niques proposed in [16], [17] can be used to compute mutual information rates. Howe ver , computing mutual information rates of noisy 2-D RLL constraints has been an unsolved problem. In this paper , the goal is to apply the generalized belief propagation (GBP) algorithm [10] for the above-mentioned problems, namely , to compute an estimate of the capacity of noiseless 2-D and 3-D RLL constrained channels and mutual information rates of noisy 2-D constrained channels. F or both problems GBP turns out to yield very good approximate results. Preliminary versions of the material of this paper have appeared in [11] and [12]. In [11], we applied GBP to compute the noiseless capacity of 2-D and 3-D RLL constrained chan- nels. In [12], GBP was applied to compute mutual information rates of a 2-D (1 , ∞ ) -RLL constrained channel with relatively small size and only at high signal-to-noise ratio (SNR). In this paper , we show that both problems reduce to estimating the partition function of graphical models with cycles. W e then apply GBP to both problems and consider new constraints and larger sizes of grid. Our main motiv ations for this research were the successful application of GBP for information rates of 2-D finite-state channels with memory in [13], Kikuchi approximation for decoding of LDPC codes and partial-response channels in [14], and tree-based Gibbs sampling for the noiseless capacity and information rates of 2-D constrained channels in [12], [15]. The outline of the paper is as follo ws. In Section II, we consider the problem of computing the partition function and discuss how this problem is related to computing the noiseless capacity and information rates of constrained channels. Region graphs, GBP , and re gion-based free energy are outlined in Section III. Section IV discusses the capacity of noiseless 2-D constraints. Numerical v alues and simulation results for the capacity of noiseless 2-D and 3-D RLL constraints are reported in Section IV -A. In Section V, we apply GBP to compute mutual information rates of noisy 2-D RLL constraints and report numerical experiments for mutual information rates in Section V -A. I I . P R O B L E M S E T - U P Consider a 2-D channel of size N = m × m with a set of X = { X 1 , X 2 , . . . , X N } random variables. Let x i denote a realization of X i and let x denote { x 1 , x 2 , . . . , x N } . W e assume that each X i takes values in a finite set X i . Also let X be the Cartesian product X 4 = X 1 × X 2 × . . . × X N . In constrained channels, not all sequences of symbols from the channel alphabet X are admissible. Let S X ⊂ X be the set of admissible input sequences. W e define the indicator function f ( x ) 4 = ( 1 , x ∈ S X 0 , x / ∈ S X (3) The partition function Z is defined as Z 4 = X x ∈X f ( x ) . (4) W ith the abov e definitions, Z = |S X | is the number of sequences that satisfy a gi ven constraint. Therefore, computing the capacity of constrained channels as expressed in (2), is closely related to computing the partition function as in (4). Also note that with the above definitions p ( x ) = f ( x ) Z (5) is a probability mass function on X . For a noisy 2-D channel, let X be the input and Y = { Y 1 , Y 2 , . . . , Y N } be the output of the channel. The mutual information rate is 1 N I ( X ; Y ) = 1 N  H ( Y ) − H ( Y | X )  . (6) Let us suppose that H ( Y | X ) is analytically available. In this case, the problem of estimating the mutual information rate reduces to estimating the entropy of the channel output, which is H ( Y ) = − E  log p ( Y )  . (7) As in [16], we can approximate the expectation in (7) by drawing L samples y (1) , y (2) , . . . , y ( L ) according to p ( y ) and use the empirical av erage as H ( Y ) ≈ − 1 L L X ` =1 log( p ( y ( ` ) )) . (8) 3 Therefore, the problem of estimating the mutual information rate reduces to computing p ( y ( ` ) ) for ` = 1 , 2 , . . . , L . W e will compute p ( y ( ` ) ) based on p ( y ( ` ) ) = X x ∈X p ( x ) p ( y ( ` ) | x ) , (9) which for a fixed y ( ` ) has also the form (4) and therefore requires the computation of a partition function. RLL constraints impose restrictions on the values of vari- ables that can be v erified locally . F or example, in a 2-D (1 , ∞ ) - RLL constraint no two (horizontally or vertically) adjacent variables can both ha ve the value 1 . The indicator function of this constraint factors into a product of kernels of the form κ a ( x i , x j ) = ( 0 , if x i = x j = 1 1 , else, (10) with one such kernel for each adjacent pair ( x i , x j ) . The factorization with kernels as in (10) can be represented with a graphical model. In this paper, we focus on graphical models defined in terms of F orne y factor graphs . Fig. 1 shows the Forney factor graph of a 2-D (1 , ∞ ) -RLL constraint where the box es labeled “ = ” are equality constraints [19]. (Fig. 1 may also be vie wed as a f actor graph as in [18] where the boxes labeled “ = ” are the variable nodes). In general, we suppose that the indicator function f ( x ) of an RLL constraint factors into a product of non-negati ve local kernels each having some subset of x as arguments; i.e. f ( x ) = Y a f a ( x a ) , (11) where x a is a subset of x and each kernel f a ( x a ) has elements of x a as arguments. In this case, the partition function in (4) can be written as Z = X x ∈X Y a f a ( x a ) . (12) If the factorization in (11) yields a cycle-free factor graph (with not too many states), the sum in (12), or equiv alently the sum in (4), can be computed efficiently by the sum-product message passing algorithm [18]. Howe ver , for the examples we study in this paper, like the Forney factor graph of a 2-D (1 , ∞ ) -RLL constraint in Fig. 1, factor graphs contain (many short) cycles. In such cases computing Z requires a sum with an exponential number of terms and therefore we are interested in applying approximate methods. Due to the presence of many short cycles in the factor graph representation of 2-D and 3-D RLL constraints, loopy belief propagation often fails to conv erge. As a result, we apply GBP = = = = = = = = = = = = = = = = Fig. 1. Forney factor graph for a 2-D (1 , ∞ ) -RLL constraint. to estimate Z , which then leads to estimating the noiseless capacity and mutual information rates of RLL constraints. I I I . G B P A N D T H E R E G I O N G R A P H M E T H O D In statistical physics, Z defined in (4) is known as the partition function and the Helmholtz free ener gy is defined as F H 4 = − ln( Z ) . (13) The partition function and the Helmholtz free energy are important quantities in statistical physics since the y carry information about all thermodynamic properties of a system. A number of techniques have been dev eloped in statistical physics to approximate the free energy . The method that we apply in this paper is known as the region-based free energy approximation, in particular we use the cluster v ariation method to select a valid set of regions and counting numbers, see [10] and [20] for more details. W e start by introducing the region graph representation of our problem. Such a region graph will provide a graphical framew ork for GBP algorithm. For each RLL constraint, the size of the basic region is chosen based on the constraint parameters. For a 2-D ( d 1 , k 1 , d 2 , k 2 ) -RLL constraint with finite k 1 and k 2 , the width and the height of the basic region is chosen as W R = k 1 + 1 H R = k 2 + 1 , and for the infinite case, the size is chosen as W R = d 1 + 1 H R = d 2 + 1 . 4 = x 7 f L x 8 = x 8 f M x 9 = x 7 x 8 x 9 f H f I f K x 4 x 5 x 6 = x 4 f F x 5 = x 5 f G x 6 = x 4 x 5 x 6 f C f D f E x 1 x 2 x 3 = x 1 f A x 2 = x 2 f B x 3 = Fig. 2. Basic region of size 2 × 2 for a 2-D (1 , ∞ ) -RLL constraint. Such a choice for the basic regions seems plausible since the validity of a given array can be determined by verifying the constraints in each region and sliding the basic regions along the rows and along the columns of the array . F or a 2-D (1 , ∞ ) -RLL constraint, Fig. 2 sho ws the basic regions and Fig. 3 shows the region graph and the counting numbers associated with each region. After forming the region graph using the cluster variation method, we perform GBP on this graph by sending messages between the regions while performing e xact computations inside each region. W e will need the region-based free energy to estimate the number of arrays that satisfy a gi ven constraint. Therefore, we operate GBP on the corresponding region graph until con ver gence and use the obtained region beliefs { b R ( x R ) } to compute the region-based free energy ˆ F H (as an estimate of F H ). The region-based free energy ˆ F H can then be used to estimate the partition function Z using (13). W e compute ˆ F H as ˆ F H = min { b R } F R ( { b R ( x R ) } ) = X R ∈R c R X x R b R ( x R )  ln b R ( x R ) − ln Y a ∈ A R f a ( x a )  (14) Here R denotes the set of all regions, c R is the counting number , x R stands for the set of variables in region R , and A R is the set of factors in region R . See Fig. 3. I V . C A PAC I T Y O F N O I S E L E S S 2 - D R L L C O N S T R A I N T S For a 2-D RLL constrained channel of width m and of size N = m × m , we run GBP on the corresponding re gion graph to compute ˆ F H and an estimate of Z . W e can then compute C ( m, m ) = log 2 Z ( m, m ) m × m , (15) f A f C f D f F x 1 x 2 x 4 x 5 +1 f B f D f E f G x 2 x 3 x 5 x 6 +1 f F f H f I f L x 4 x 5 x 7 x 8 +1 f G f I f K f M x 5 x 6 x 8 x 9 +1 ? H H H H H H H j         H H H H H H H j         H H H H H H H j         ? f D x 2 x 5 − 1 f F x 4 x 5 − 1 f G x 5 x 6 − 1 f I x 5 x 8 − 1 H H H H H H H H j J J J J ^               x 5 +1 Fig. 3. The region graph for Forney factor graph in Fig. 2. where Z ( m, m ) denotes the number of 2-D binary arrays of size m × m that satisfy the constraint. In our numerical experiments in Section IV -A, for dif ferent RLL constraints we show con vergence of C ( m, m ) to the Shannon capacity as m increases. For example, let us consider a 2-D (1 , ∞ ) -RLL constraint with corresponding Forney factor graph in Fig. 1. For this constraint, we chose basic regions with size 2 × 2 in a sliding window manner ov er the factor graph, see Fig. 2. Starting from such basic regions, we applied the cluster variation method on the factor graph in Fig. 2 to obtain the corresponding region graph depicted in Fig. 3. The counting numbers { c R } are shown next to each region. A. Numerical Experiments Here we present the numerical results of applying GBP to estimate the finite-sized noiseless capacity of RLL constraints. T ight lower and upper bounds were gi ven for the Shannon capacity of a 2-D (1 , ∞ ) -RLL constraint in [4]. The bounds were further improved in [22] and [23], now known to nine decimal digits. 0 . 5878911617 ... ≤ C (1 , ∞ ) 2 D ≤ 0 . 5878911618 ... (16) For this constraint, Fig. 4 shows C ( m, m ) defined in (15) versus the channel width m ov er the interval [2 , 300] . The estimation was performed using the parent-to-child and two- way GBP algorithms. The two algorithms giv e almost identical results. The horizontal line in Fig. 4 shows the Shannon capacity for this channel in (16). For a channel of width 300 , the estimated noiseless capacity is about 0 . 5884 . 5 0.58 0.6 0.62 0.64 0.66 0.68 0.7 0.72 2 4 6 8 10 20 30 40 50 100 150 200 300 bits/symbol Channel width m Fig. 4. Estimated capacity (in bits per symbol) vs. channel width m for a 2-D (1 , ∞ ) -RLL constraint. The horizontal dotted line shows the Shannon capacity for this channel as in (16). 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 2 3 4 5 6 8 10 20 30 40 50 100 150 200 bits/symbol Channel width m 2D RLL Channels (1, ∞ ) (1, ∞ , 2, ∞ ) (1, ∞ , 3, ∞ ) (1, ∞ , 4, ∞ ) Fig. 5. Estimated capacities (in bits per symbol) vs. channel width m for a class of 2-D (1 , ∞ , d, ∞ ) -RLL constraints with d = (1 , 2 , 3 , 4) . Shown in Fig. 5 are plots of C ( m, m ) for 2-D (1 , ∞ , d, ∞ ) - RLL constraints with d = (1 , 2 , 3 , 4) from top to bot- tom, versus the channel width m over the interv al [2 , 200] . Fig. 6 sho ws the plots of C ( m, m ) for 2-D (1 , ∞ , 2 , 4) - RLL and (1 , ∞ , 2 , 3) -RLL constraints versus m over the interval [4 , 200] . From our simulation results, for a chan- nel of width 200 the estimated noiseless capacities for 2-D (1 , ∞ , d, ∞ ) -RLL constraints with d = (2 , 3 , 4) are about (0 . 4994 , 0 . 4346 , 0 . 3864) and the estimated noiseless capaci- ties for 2-D (1 , ∞ , 2 , 4) -RLL and (1 , ∞ , 2 , 3) -RLL are about (0 . 3106 , 0 . 2109) . T o the best of our knowledge, no theoretical upper or lo wer bounds e xist for these constraints. All plots are obtained using the parent-to-child algorithm. Note that 2-D (1 , ∞ , 1 , ∞ ) -RLL plot in Fig. 5 is the same as the plot in Fig. 4. Also shown in Fig. 7 is the plot of C ( m, m ) for a 2-D (2 , ∞ ) -RLL constraint versus m over the interval [3 , 400] . For 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 4 5 6 8 10 20 30 40 50 100 200 bits/symbol Channel width m 2D RLL Channels (1, ∞ , 2, 3) (1, ∞ , 2, 4) Fig. 6. Estimated capacity (in bits per symbol) vs. channel width m for a 2-D (1 , ∞ , 2 , 4) -RLL and (1 , ∞ , 2 , 3) -RLL constraints. 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 3 6 8 10 20 30 40 50 100 150 200 300 400 bits/symbol Channel width m Fig. 7. Estimated capacity (in bits per symbol) vs. channel width m for a 2-D (2 , ∞ ) -RLL constraint. 0.5 0.52 0.54 0.56 0.58 0.6 0.62 0.64 2 4 6 8 10 20 30 40 bits/symbol Channel width m Fig. 8. Estimated capacity (in bits per symbol) vs. channel width m for a 3-D (1 , ∞ ) -RLL constraint. The horizontal dotted lines show upper and lower bounds on the Shannon capacity for this channel as in (18). a channel of width 400 , the estimated noiseless capacity is about 0 . 4462 . Best known lo wer and upper bounds for the 6 Shannon capacity of a 2-D (2 , ∞ ) -RLL constraint are giv en in [8] and [9] respecti vely , as 0 . 4453 ≤ C (2 , ∞ ) 2 D ≤ 0 . 4457 (17) Our proposed method can be generalized to compute the noiseless capacity of 3-D and higher dimensional RLL con- straints. For a 3-D (1 , ∞ ) -RLL constraint the following lower and upper bounds were introduced in [23] 0 . 5225017418 ... ≤ C (1 , ∞ ) 3 D ≤ 0 . 5268808478 ... (18) Fig. 8 shows the noiseless capacity estimates of a 3- D (1 , ∞ ) -RLL constraint, obtained using the parent-to-child algorithm, versus the channel width m . The horizontal dotted lines sho w the upper and lower bounds for the Shannon capacity . For a channel of width m = 40 the GBP estimated capacity is about 0 . 5267 which falls within these bounds. Simulation results and numerical values for the noiseless capacity of many other 2-D RLL constraints are reported in [21]. B. Bounds for the Shannon Capacity For any finite m , it is possible to compute lower and upper bounds on the Shannon (infinite-size) capacity using C ( m, m ) the capacity of a 2-D RLL constrained channel of width m . For example, consider a 2-D (1 , ∞ ) -RLL constraint with local kernels as in (10). From tiling the whole plane with m × m squares, it is clear that C ( m, m ) is an upper bound for the Shannon capacity C (1 , ∞ ) 2 D . On the other hand, by tiling the plane with m × m squares separated by all-zero guard rows and all-zero guard columns, we obtain ( m m +1 ) 2 C ( m, m ) ≤ C (1 , ∞ ) 2 D . From Fig. 4, the estimated capacity at m = 300 is about C (300 , 300) = 0 . 5884 , we thus obtain the following lo wer and upper bounds for the Shannon capacity 0 . 5844 ≤ C (1 , ∞ ) 2 D ≤ 0 . 5884 . Note that although GBP performs remarkably well for 2- D constrained channels, it is an approximate algorithm which yields approximations to the lower and upper bounds to the Shannon capacity . Howe ver in order to achieve a desired precision, the bounds could provide a criterion for choosing the value of m . V . I N F O R M AT I O N R A T E S O F N O I S Y 2 - D R L L C O N S T R A I N T S As e xplained in Section II, the problem of computing mutual information rates reduces to computing the output probability . = Z Z Z Z y 1 x 1 = Z Z Z Z y 2 x 2 = Z Z Z Z y 3 x 3 = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z = Z Z Z Z Fig. 9. Extension of Fig. 1 to a Forney factor graph of p ( x , y ) with p ( y | x ) as in (22). Therefore, the remaining tasks are 1) Drawing input samples x (1) , x (2) , . . . , x ( L ) from S X ac- cording to p ( x ) and therefrom creating output samples y (1) , y (2) , . . . , y ( L ) . 2) Computing p ( y ( ` ) ) for each ` = 1 , 2 , . . . , L . W e will compute p ( y ( ` ) ) based on p ( y ( ` ) ) = X x ∈S X p ( x ) p ( y ( ` ) | x ) , (19) where p ( x ) is a probability mass function on S X . Let us assume uniform distribution over the admissible channel input configurations. Therefore we have p ( x ) = |S X | − 1 (20) = 2 − N C 2 D , (21) we also assume the channel is memoryless and p ( y | x ) factors as p ( y | x ) = N Y i =1 p ( y i | x i ) . (22) For such a noisy 2-D constrained channel, the corresponding Forne y factor graph, as an extension of Fig. 1, is sho wn in Fig. 9. Using (21) and (22), we can rewrite (19) as p ( y ( ` ) ) = 2 − N C 2 D X x ∈S X N Y i =1 p ( y ( ` ) i | x i ) , (23) = 2 − N C 2 D Z ( y ( ` ) ) , (24) where Z ( y ( ` ) ) has the same form as the sum in (12). The input samples x (1) , x (2) , . . . , x ( L ) are generated as fol- lows. W e run GBP on Fig. 3 until con ver gence to compute 7 the region beliefs { b R ( x R ) } at each region R . The region be- liefs are GBP approximations to the corresponding marginals { p R ( x R ) } . In our numerical experiments, each sample x ( ` ) is then generated piecewise sequentially according to the beliefs b R ( x R ) in basic regions. For example, in the region graph of Fig. 3, after computing b R ( x 1 , x 2 , x 4 , x 5 ) , sample x 1 is drawn according to b R ( x 1 ) , sample x 2 is drawn according to b R ( x 2 | x 1 ) , etc. The input samples x (1) , x (2) , . . . , x ( L ) are then used to create output y (1) , y (2) , . . . , y ( L ) using (22). The beliefs are directly proportional to the factor nodes in volv ed in each region, which guarantees that the samples are drawn from S X . Moreov er, since beliefs are good approx- imations to the marginal probabilities, one expects that the samples are drawn from a distribution close to p ( x ) , see [10]. In order to compute Z ( y ( ` ) ) , as in Section III, we start from the factor graph in Fig. 9 to build the region graph representing the problem and run GBP on this region graph. Finally , the estimated p ( y (1) ) , p ( y (2) ) , . . . , p ( y ( L ) ) are used to compute an estimate of H ( Y ) as in (8). A. Numerical Experiments In our numerical experiments we consider (1 , ∞ ) -RLL and (2 , ∞ ) -RLL constrained channels with size N = 30 × 30 and input alphabet X = {− 1 , +1 } N . Noise is assumed to be i.i.d. zero mean Gaussian with variance σ 2 and independent of the input. W e thus hav e H ( Y | X ) = N 2 log(2 π eσ 2 ) , (25) and p ( y | x ) in (22) has kernels of the form p ( y i | x i ) = 1 √ 2 π σ 2 exp  − 1 2 σ 2  y i − x i  2  . (26) SNR is defined as SNR 4 = 10 log 10 1 σ 2 (27) Shown in Fig. 10 is the estimated information rate vs. SNR ov er the interv al [ − 10 , 10] dB for a noisy 2-D (1 , ∞ ) -RLL constraint. The horizontal dotted line sho ws the estimated noiseless capacity which can be read from Fig. 4 and is about 0 . 5943 for this size of channel. Illustrated in Fig. 11 is the estimated information rate vs. SNR over the interv al [ − 10 , 10] dB for a noisy 2-D (2 , ∞ ) -RLL channel. The horizontal dotted line sho ws the estimated noiseless capacity which can be read from Fig. 7 and is about 0 . 4552 for this size of channel. The simulation results were obtained by averaging ov er L = 1000 realizations of the channel output. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 -10 -8 -6 -4 -2 0 2 4 6 8 10 bits/symbol dB Fig. 10. Estimated information rate (in bits per symbol) vs. SNR (in dB) for a 30 × 30 channel with a (1 , ∞ ) -RLL constraint and additi ve white Gaussian noise. 0 0.1 0.2 0.3 0.4 0.5 -10 -8 -6 -4 -2 0 2 4 6 8 10 bits/symbol dB Fig. 11. Estimated information rate (in bits per symbol) vs. SNR (in dB) for a 30 × 30 channel with a (2 , ∞ ) -RLL constraint and additi ve white Gaussian noise. Simulation results and numerical v alues for mutual infor - mation rates of many other 2-D RLL constraints are reported in [21]. V I . C O N C L U D I N G R E M A R K S W e proposed a GBP-based method to estimate the noiseless capacity and mutual information rates of RLL constraints in two and three dimensions. For noiseless RLL constraints, the method was applied to estimate the finite-size capacity of different constraints and to show conv ergence to the Shannon capacity as the size of the channel increases. In particular, the proposed method can be used to estimate the noiseless capacity of RLL constraints in the cases that the capacity is not kno wn to a useful accuracy . The method was also applied to estimate mutual information rates of noisy RLL constraints with additiv e white Gaussian noise and with a uniform distri- bution o ver the admissible input configurations. Our simulation 8 results show mutual information rates of different constraints as a function of SNR. A C K N O W L E D G E M E N T S The authors gratefully acknowledge the support of Prof. H.-A. Loeliger . The first author wishes to thank Ori Shental for his helpful comments on GBP implementation. W e would also like to thank the revie wers for their many helpful suggestions that helped to improve the presentation of our paper . R E F E R E N C E S [1] K. A. Schouhamer Immink, Codes for Mass Data Storag e Systems. Eindhoven: Shannon Foundation Publishers, 2004. [2] C. E. Shannon, “ A mathematical theory of communications, ” Bell Sys. T ech. Journal, vol. 27, pp. 379–423, July 1948. [3] P . H. Sie gel, “Information-theoretic limits of two-dimensional optical recording channels, ” in Optical Data Storage (Proc. of SPIE, V ol. 6282, Eds. Ryuichi Katayama and T uviah E. Schlesinger), Montreal, Quebec, Canada, April 2006. [4] N. J. Calkin and H.S. W ilf, ‘The number of independent sets in a grid graph, ” SIAM J. Discr . Math., vol. 11, pp. 54–60, Feb. 1998. [5] S. Forchhammer and T .V . Laursen, “Entropy of bit-stuffing-induced measures for two-dimensional checkerboard constraints, ” IEEE T rans. Inform. Theory , vol. 53, pp. 1537–1546, April 2007. [6] H. Ito, A. Kato, Z. Nagy , and K. Zeger , ‘Zero capacity region of multidimensional run length constraints, ” The Electr onic Journal of Combinatorics, vol. 6(1), 1999. [7] K. Kato and K. Ze ger, ‘On the capacity of two-dimensional run-length constrained channels, ” IEEE T rans. Inform. Theory , vol. 45, pp. 1527– 1540, July 1999. [8] E. Ordentlich and R. M. Roth, “ Approximate enumerative coding for 2- D constraints through ratios of matrix products, ” Pr oc. 2009 IEEE Int. Symp. on Information Theory, Seoul, Korea, pp. 1050–1054. [9] I. T al and R. M. Roth, “Concave programming upper bounds on the capacity of 2-D constraints, ” IEEE T rans. Inform. Theory , v ol. 57, pp. 381–391, Jan. 2011. [10] J. S. Y edidia, W . T . Freeman, and Y . W eiss. “Constructing free energy approximations and generalized belief propagation algorithms, ” IEEE T rans. Inform. Theory, vol. 51, pp. 2282–2312, July 2005. [11] G. Sabato and M. Molkaraie, “Generalized belief propagation algorithm for the capacity of multi-dimensional run-length limited constraints, ” Pr oc. 2010 IEEE Int. Symp. on Information Theory , Austin, USA, June 13–18, pp. 1213–1217. [12] M. Molkaraie and H.-A. Loeliger , “Estimating the information rate of noisy constrained 2-D channels, ” Pr oc. 2010 IEEE Int. Symp. on Information Theory , Austin, USA, June 13–18, pp. 1678–1682. [13] O. Shental, N. Shental, S. Shamai (Shitz), I. Kanter, A. J. W eiss, and Y . W eiss, “Discrete-input two-dimensional Gaussian channels with memory: estimation and information rates via graphical models and statistical mechanics, ” IEEE T rans. Inform. Theory , vol. 54, pp. 1500– 1513, April 2008. [14] P . Pakzad and V . Anantharam, “Kikuchi approximation method for joint decoding of LDPC codes and partial-response channels, ” IEEE T rans. Communications, vol. 54, pp. 1149–1153, July 2006. [15] H.-A. Loeliger and M. Molkaraie, “Estimating the partition function of 2-D fields and the capacity of constrained noiseless 2-D channels using tree-based Gibbs sampling, ” Proc. 2009 IEEE Information Theory W orkshop, T aormina, Italy , Oct. 11–16, pp. 228–232. [16] D. Arnold, H.-A. Loeliger , P . O. V ontobel, A. Kav ˇ ci ´ c, and W . Zeng, “Simulation-based computation of information rates for channels with memory , ” IEEE Tr ans. Inform. Theory , vol. 52, no. 8, pp. 3498–3508, August 2006. [17] H. D. Pfister, J.-B. Soriaga, and P . H. Siegel, “On the achie vable information rates of finite-state ISI channels, ” in Proc. 2001 IEEE Globecom , San Antonio, USA, Nov . 25–29, pp. 2992–2996. [18] F . R. Kschischang, B. J. Frey , and H.-A. Loeliger, “F actor graphs and the sum-product algorithm, ” IEEE T rans. Inform. Theory, vol. 47, pp. 498– 519, Feb. 2001. [19] H.-A. Loeliger , “ An introduction to factor graphs, ” IEEE Signal Proc. Mag., Jan. 2004, pp. 28–41. [20] M. W elling, “On the choice of regions for generalized belief propa- gation, ” Proc. 2004 confer ence on Uncertainty in artificial intelligence, Banff, Canada, pp. 585–592. [21] G. Sabato, Simulation-based techniques to study two-dimensional ISI channels and constrained systems. Master thesis, Dept. Inform. T echn. & Electr . Eng, ETH Z ¨ urich, Switzerland, 2009. [22] W . W eeks and R. E. Blahut, “The capacity and coding gain of certain checkerboard codes, ” IEEE T rans. Inform. Theory , vol. 44, pp. 1193– 1203, May 1998. [23] Z. Nagy and K. Zeger, ‘Capacity bounds for the three-dimensional (0 , 1) run-length limited channels, ” IEEE Tr ans. Inform. Theory, vol. 46, pp. 1030–1033, May 2000.