Eric

Jules Julien

Marie-Capucine Guillaume

Lorine Brice

Juliette Lucie

Maxime Germain

Keywords: factor graph, LDPC codes, variational free energy, region-based approximation, trapping sets, chaos. xi

En dehors du laboratoire, j'ai été soutenu par mes parents, mes beauxparents, mes frères et mes belles-soeurs que je remercie tous chaleureusement. Je remercie avec grande joie mes cousins Nathan, Anna, Tom et Maïa avec qui j'ai beaucoup parlé et relativisé les obstacles de la thèse, mes cousines Aurélie et Sarah pour leurs encouragements, et tout spécialement mon cousin Philippe qui, depuis les études supérieures me suit, me conseille et me soutient avec beaucoup d'affection, qui a passé beaucoup de temps avec moi pour discuter des solutions possibles aux impasses que je rencontrais, qui a su lire et critiquer avec précision les écrits que je produisais ainsi que le présent manuscrit.

posé cette thèse et qui m'ont laissé libre dans mes choix d'axes de recherche, Jean-Pierre Tillich, qui m'a accueilli à l'INRIA ainsi que dans le projet COCQ, et qui a toujours montré de l'intérêt pour mes travaux. Je remercie également les autres membres du laboratoire avec qui j'ai interagi: Frédéric Precioso qui a été très accueillant et de bons conseils, Michel Leclerc et Laurent Protois qui m'ont beaucoup aidé pour que je puisse utiliser au mieux les ressources informatiques, Annick Bertinotti qui malgré les multiples affaires administratives à régler savait prendre du temps pour nous les doctorants. Je remercie Emmanuel Bourdel et Inbar Fijalkow, ainsi que Sylvain Reynal et Michel Leclerc, pour m'avoir donné l'opportunité d'enseigner dans divers domaines, ce que j'ai beaucoup apprécié.

Les doctorants du laboratoire ont également participé, de près ou de loin, à mon travail. Je remercie Antoine Rolland de Rengervé qui savait m'écouter et faire ressortir dans mes problèmes les points importants, Alexis Lechervy qui a passé beaucoup de temps à m'aider à améliorer mes programmes et à structurer mon travail. J'ai une pensée toute particulière pour mes collègues et amis Jean-Emmanuel (au sens de l'humour si personnel mais aigu), Romain Tajan, Gaël Rigaud, et Ludovic Danjean qui ont travaillé avec moi pendant des demi-journées (voire journées) entières sur mes problèmes, qui ont lu, critiqué, corrigé et considérablement enrichi mes articles et mon travail, et qui ont été des oreilles attentives à mes détresses scientifiques. Je remercie aussi mes collègues de bureau Alan Julé pour les déjeuners sympathiques et les gâteaux bretons, Rémi Régnier pour la profonde connaissance et admiration que nous partagions à propos de Kaamelott, ainsi que Liang Zhou, Raouia Masmoudi, Shuji Zhao et Sonia Khatchadourian. Je remercie Erbao Li avec qui je n'ai malheureusement pas eu le temps d'écrire un article, mais qui m'a aidé à développer des idées et concepts, Fouad Sahraoui avec qui j'ai beaucoup parlé de spiritualité mais aussi de la compréhension de l'informatique et ses défis. Je remercie les autres doctorants et docteurs du laboratoire avec qui j'ai passé de bons moments : Laurent Rodriguez, Rodrigue Imad, Babar Aziz, David Gorisse, Thomas Lefebvre, Leïla Meziou, Mathilde Brandon.

iii

Version franc ¸aise

Dans cette thèse, nous étudions le problème de l'inférence bayésienne dans les graphes factoriels, en particulier les codes LDPC, quasiment résolus par les algorithmes de message-passing. Nous réalisons en particulier une étude approfondie du Belief Propagation (BP), un algorithme de message-passing dont nous démontrons l'équivalence avec l'approximation de Bethe utilisée pour étudier les verres de spins en physique statistique. La question de la sous-optimalité est soulevée dans le cas o ù le graphe factoriel présente des boucles. Ensuite est détaillé le Generalized Belief Propagation (GBP), un autre algorithme de message-passing élaboré à partir d'une généralisation de l'approximation de Bethe, l'approximation basée régions. Fondé sur un découpage en clusters du graphe factoriel, nous montrons par des expériences que le GBP peut être significativement plus performant que le BP dans les cas o ù le découpage donne naissance à des clusters très liés au structures topologiques, en particulier les trapping sets, qui empêchent le BP de bien décoder. Nous ne confrontons pas les performances des deux algorithmes uniquement par rapport à leurs taux d'erreur binaire mais aussi par rapport à leurs comportements dynamiques face à des événements d'erreur non triviaux. En utilisant des estimateurs inédits, nous montrons que l'algorithme du GBP peut dominer celui du BP en termes de précision et de convergence.

v

Résumé

Dans cette thèse, nous nous attaquons aux problèmes liés à l'inférence statistique, en particulier l'approximation des distributions marginales sur les noeuds d'un graphe factoriel à boucles. L'algorithme type message-passing du Belief Propagation (BP) créé par Pearl en 1988 ne parvient pas à approcher ces grandeurs de manière optimale, malgré de très bonnes performances. En effet, dans le cas des codes LDPC, le décodage par le BP présente des oscillations voire des comportements non classiques, comme du chaos. Après avoir présenté le fonctionnement de cet algorithme, nous montrons qu'il est équivalent à l'approximation de Bethe utilisée en physique statistique comme une approche de champ moyen pour estimer la fonction de partition des verres de spins par une minimisation de l'énergie libre variationnelle.

Les problèmes d'oscillations dans les graphes factoriels, dits aussi de frustration pour les verres de spins, peuvent être résolus grâce à l'approximation basée régions qui consiste à créer un réseau bayésien associé, appelé graphe des régions, dont les noeuds, ou régions, sont des regroupements des noeuds qui composent le graphe factoriel. Par ce moyen, il est attendu que la précision de l'approximation soit meilleure que celle donnée par l'approximation de Bethe. Par la même équivalence qui lie cette dernière au BP, nous démontrons que l'approximation basée régions peut être traduite comme un algorithme de message-passing appelé Generalized Belief Propagation (GBP). Les performances de décodage du GBP peuvent s'avérer meilleures que celles du BP pourvu que le découpage en clusters ayant donné naissance au graphe des régions, soit bien réalisé. En effet, aucun critère de choix permettant de justifier de l'utilisation d'un graphe des régions ou d'un autre n'a émergé jusque-là. Après avoir détaillé les équations de mise à jour du GBP, nous détaillons la construction du graphe des régions et nous insistons sur le fait que, conformément à son existence, il doit gérer les structures topologiques du graphe factoriel qui posent le plus de problèmes. Dans le cadre du décodage des codes LDPC, modélisables graphiquement par un graphe dit de Tanner, de telles structures ont été répertoriées dans la littérature. Il s'agit des trapping sets. Responsables de la dégradation du taux d'erreur binaire à très forts rapports signal-à-bruit (RSB), nous utilisons en conséquence ces structures pour réaliser un découpage original et pertinent du graphe de Tanner. Nous précisons que le GBP étant un algorithme naturellement très instable, nous introduisons un amortissement dans toutes nos simulations dont la loi d'évolution force l'algorithme à converger. vii viii Par des résultats expérimentaux sur un code LDPC particulier entièrement couvert par des trapping sets, nous montrons que le GBP muni d'un amortissement a de meilleures performances que le BP pour de très forts RSB.

La comparaison du BP et du GBP pour des RSB moyens en termes de taux d'erreur binaire n'est pas concluante, les deux algorithmes étant équivalents. Nous utilisons donc les outils de la théorie des systèmes dynamiques pour tenter de les départager. Le champ de recherche lié à l'étude de la dynamique des décodeurs itératifs pour les codes LDPC est récente, par conséquent la littérature à ce sujet n'est pas très fournie. Les articles publiés sur la dynamique du BP traitent de l'utilisation des fonctions classiques liées aux systèmes chaotiques: le diagramme de bifurcation pour obtenir les RSB correspondant à de brusques changements de comportement et l'exposant de Lyapunov pour juger de la stabilité, i.e. de la sensibilité aux conditions initiales. Le GBP n'a en revanche jamais été traité de ce point de vue, faisant de notre étude un cas nouveau. A l'aide de ces fonctions mais également d'outils nouveaux que nous introduisons, nous réalisons une étude expérimentale afin de connaître les types d'attracteurs auxquels sont soumis le BP et le GBP ainsi que leur relation sur leurs performances de décodage. Il apparaît que les deux décodeurs ont tous deux des comportements typiques similaires: convergence vers un point fixe, attraction vers un cycle limite, attraction vers un attracteur chaotique. Cependant, les propriétés de ces attracteurs nous permettent d'affirmer que le GBP a une tendance plus prononcée que le BP à converger vers un bon décodage.

Résultats principaux et contributions:

• utilisation des trapping sets pour la construction du graphe des régions,

• introduction d'un coefficient d'amortissement,

• mise en place d'une méthode pour l'étude des décodeurs à très forts RSB,

• étude de la dynamique du GBP,

• introduction d'outils originaux pour la description de la dynamique des décodeurs en lien avec leurs performances de décodage.

Cette thèse est également un travail original au sens o ù elle est à la croisée entre divers champs de recherche: les codes correcteurs d'erreurs, la physique statistique, l'inférence statistique, l'optimisation lagrangienne et la théorie des sytèmes dynamiques. C'est pourquoi la première partie du manuscrit est dédiée à une revue détaillée des méthodes et outils issus de ces champs de recherche, tandis que la seconde partie traite de l'explication détaillée des deux algorithmes BP et GBP ainsi que de leur étude d'un point de vue dynamique.

Mots-clés : graphe factoriel, codes LDPC, énergie libre variationnelle, approximation basée régions, trapping sets, chaos.

Abstract

In this thesis, we tackle statistical inference on loopy factor graphs to approximate marginal probability distributions on nodes. In 1988, Pearl introduced a message-passing algorithm, called the Belief Propagation algorithm (BP), to solve this problem, but in spite of very good performance, BP did not offer optimal approximations. As a matter of fact, when applied on LDPC codes decoding, BP may present oscillations or even chaotic behaviors. Following the presentation of BP equations, we show that it is equivalent to Bethe approximation, used in statistical physics as a mean field approach, to estimate partition function of spin glasses by minimizing the variational free energy.

Oscillations of factor graphs, also called frustrations for spin glasses, may be solved by region-based approximation. It consists in creating an associated Bayesian network, called region-graph, which nodes, or regions, are gathering of factor graph nodes. This way, accuracy of this new approximation is hoped to be better than Bethe approximation's. According to equivalence between Bethe approximation and BP, we demonstrate that region-based approximation may be derived as a message-passing algorithm, called Generalized Belief Propagation (GBP). GBP decoding performance may be better than BP's provided that the clustering used to gather factor graph nodes is well suited. In fact, no choice criterion has been introduced to discriminate between all possible clusterings. After demonstrating GBP update equations, we present in detail the region-graph construction, and we emphasize that, in accordance with its origin, it has to deal with harmful topological structures of factor graph. LDPC codes, represented by particular factor graphs called Tanner graphs, present such structures, that have been published in the literature, called trapping sets. They are responsible for degradation of the Bit-Error Rate (BER) for high values of the Signal-to-Noise Ratio (SNR). We accordingly use these structures to cluster, in an original and relevant way, Tanner graphs. We specify that the GBP algorithm is naturally instable, which makes us introduce damping factors, of various evolutions, that help decoder converge. We show with experimental results on a specific LDPC code, entirely covered with trapping sets, that a dedicated construction coupled with a damped GBP offer better performance than BP for high SNR values. ix x Comparing BP and GBP for middle SNR values only using BER is not conclusive as both decoders exhibit equivalent results. Tools of dynamical systems theory open opportunities to distinguish them. Research area on dynamics of iterative LDPC decoders has been recently introduced making the literature on this topic not very abundant. Published papers on BP dynamics deal with usual functions that describe chaotic systems: bifurcation diagram offers SNR values where the decoder abruptly changes its behavior, Lyapunov exponent informs about BP stability, i.e. sensitivity to initial conditions. GBP has never been studied from this point of view, making this work innovative. By means of usual functions but also of new tools, we carry out an experimental investigation to bring out BP and GBP attractors and their links with decoding performance. Both algorithms similarly behave as they present close attractors for same SNR values: fixed point, limit cycle, chaos. However, attractors properties allow us to affirm that GBP tends to converge to correct decoding better than BP.

Main results and contributions:

• use of trapping sets construct region-graphs,

• introduction of a general damping factor for GBP,

• development of an experimental protocol to study decoders for high SNR values,

• study of GBP dynamics,

• introduction of innovative tools to describe decoders dynamics and to link them to their decoding performance.

The originality of this thesis remains in the fact that is is at the crossroads of several research areas: error-correcting codes theory, statistical physics, statistical inference, Lagrange optimization and dynamical systems theory. Therefore, first part of the manuscript is dedicated to reviews on methods and tools from most of these fields, and second part deals with foundations of the BP and GBP algorithms as well as the study of their dynamics.

Introduction

History and motivations

In this thesis we expose the analysis and the developments of methods to construct and solve graphical models. Many engineering and research problems can be formulated as a statistical inference wording, provided that they are modeled as Bayesian networks, Markov random fields or factor graphs. These problems mostly take part in activities about image processing, neural networks, spin glasses, error-correcting codes, etc. All of them are aimed at extracting either the joint probability distribution of their graphical model or the marginal distributions of the variables that build the graph. This way, it can be extracted a large number of properties that make sense for any of the concerned research areas, that helps construct strong connections between them.

The spin glasses are physical networks naturally represented by factor graphs. Statistical physics has the great challenge to obtain thermodynamical functions at equilibrium, as average energy, free energy, entropy, pressure etc. It has been demonstrated that all these functions are analytically dependent on a single quantity that is the partition function Z. The partition function has its origin in the Boltzmann's law that provides the joint distribution of any spin glass. Indeed, Z is a constant that ensures that the Boltzmann's law is normalized, i.e. it has to be a true distribution. It has been proved that computing Z consists in spanning of all possible states of the spin glass, that represents a considerable number of computations. Partition function is accordingly intractable, therefore it became a central topic of research in statistical physics. Physicists created original methods to approximate its value in particular cases: Sherrington-Kirkpatric model, Edwards-Anderson model, etc. The main method still studied nowadays is the Mean-Field (MF) approximation. The core of this method is to consider that any spin of a given spin glass is not singularly correlated to each of its neighbors anymore but to a mean field. The energy function, that summarizes all spins couplings, is then modified. Boltzmann's law assigns a probability weight to each state of the spin glass according to its energy, therefore the MF approximation entails a modification of the joint distribution. It turns out to be the product of all marginal distributions of probability. This mean field distribution is the basis to approximate the partition function by minimizing the mean-field variational free energy, whose minimum is the log-partition function.

By means of the statistical physics, the joint distribution of any graphical model could be approximated thanks to the MF approach. In addition, the notions of energy and entropy became familiar in other research areas. For instance, researchers connected parity-check equations of error-correcting codes to the spin glass energy function, and they succeed in relating the information and thermodynamical entropies. Low-Density Parity-Check (LDPC) codes, created by R. Gallager in 1963 and known to approach Shannon's limit, were then given a physical meaning. In the research area on error-correcting codes, LDPC codes are one of the main topics of research. Their links with spin glasses still spark much interest, as shown by the numerous papers published each year on this issue. The graphical model of any LDPC code is a Tanner graph, a specific factor graph. Edges of any factor graph are used as supports to solve inference by the use of message-passing algorithms, e.g. the Belief Propagation algorithm ( BP) invented by J. [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF]. This algorithm has proved to be the best compromise between computation time and accuracy. Beyond energy, entropy and other physical quantities that are shared with statistical physics, it turned out that specific methods as BP are also available in both research fields. Indeed, fixed-points of BP are strictly equivalent to stationary points of the Bethe approximation, an extension of the MF approach. This equivalence helped researchers define BP as an iterative procedure to minimize the Bethe variational free energy F β of any factor graph. In other words, the general inference problem was demonstrated to be identical to an optimization problem solved by the Lagrangian formalism.

At the origin, J. Pearl specified that BP was a perfect inference algorithm for any graphical model whose topology is tree-like. In addition, it was demonstrated in papers of T. Heskes in the 2000s that the optimization of F β of any factor graph was always possible if it contains at most one loop. Unfortunately, this remarkable work and the related results based on the convexification of F β could not be extended to any loopy factor graph. At this time, the challenge was to modify the BP algorithm to overcome the undesired loopy effects. Several methods were proposed, e.g. the ConCave-Convex Procedure (CCCP) from A. Yuille, but they implied a dramatic increase in the computation time. Once again, connections with spin glasses turned out to be conclusive. As a matter of fact, the problem of loops especially in BP decoding of the LDPC codes was perfectly understandable as the problem of frustration in statistical physics. Actually, not only the quantities and the methods but also the drawbacks are shared with the statistical physics, that helps to solve few issues. Physicists could make use of the cluster variation method or region-based approximation invented by Kikuchi in 1951 to deal with the frustration in spin glasses. This method, also called Kikuchi approximation after its designer, was eventually a MF approach that generalized the Bethe approximation by clustering any spin glass. Rather than factorizing the joint distribution on the single spins marginal distributions, it is approximated by the product of estimates of the pseudo-marginal distributions, namely the beliefs, on subsets of spins, called regions. It resulted in a new variational free energy F K , which minimization by Lagrangian formalism provided a message-passing algorithm, fairly detailed in 2001 by J. Yedidia et al: the Generalized Belief Propagation algorithm (GBP). The support of this algorithm is the region-graph: a Bayesian network constituted by regions.

The relevancy of the region-based approximation relies on the significant degree of freedom to cluster any factor graph. Therefore, it appeared that the GBP algorithm would be an efficient substitute to the BP algorithm to annihilate the frustration, i.e. to cancel the harmful loops effect in the decoding of the LDPC codes.

At this point comes the main issue of the region-based approximation. Indeed, no criterion was given by Kikuchi to discriminate the numerous clusterings of any factor graph. Thus, given a factor graph, no one is able to judge the relevancy of an arbitrary construction of the region-graph. J. Yedidia et al introduced two principles that any region-graph should fulfill to present challenging results. The first states that the energy function of the region-graph should strictly equal the energy function of the basic factor graph. Accordingly, the number of occurrences of a given spin, or node of the factor graph, has to be identical in both energy functions. The second principle relies on the fact that the Kikuchi approximation evaluates the entropy by a regionbased entropy, that is a linear decomposition over the regions entropies of the region-graph. It can be shown that even in the case of a tree-like topology, this decomposition remains an approximation. The best that can be done is to verify that making all the beliefs uniform distributions implies that the region-based entropy is maximum. Despite the power of this principle, it is again intractable to check for any region-graph. Thus the only way to evaluate the advantage of a region-graph is to run the GBP algorithm and to analyze the corresponding results versus BP's. This is the core of the work presented in this thesis.

A booming topic is the study of the dynamical behaviors of the inference algorithms, in particular in the case of the LDPC codes decoding. Few papers emerged in the last decade to deal with the bifurcations and the chaos in the decoding algorithms, especially the BP algorithm. Given that the dynamical system related to the BP algorithm is a multi-dimensional space, the mathematical developments are very hard tasks. Therefore, the study is oriented toward a use of numerous experiments of specific quantifiers to highlight typical behaviors. It appears that the BP algorithm is subject to non-trivial behaviors for particular signal-to-noise ratio range. However, the literature on the dynamics of the LDPC decoding is not substantial, given that this research domain is only in its early age. We demonstrate in this thesis that the study of the dynamics of the decoders helps compare and confront them. We introduce a toolbox to analyze the behaviors of the BP and the GBP. In this way, we emphasize the power of a region-graph construction, that provides a way to judge its relevancy for the region-based approximation.

Organization of the manuscript

The investigation we present in this manuscript lies within the areas of errorcorrecting code theory, statistical physics and dynamical systems theory. From these three fields are extracted tools and methods to develop and analyze message-passing algorithms with an extensive focus on BP and GBP. Thus we organize the manuscript as follows.

Part I is dedicated to the introducing the context in which the thesis began, digital communications. Brief notions on channel coding are presented the Belief Propagation algorithm and its derivatives, together with the optimization problem it is equivalent to. Then, we show preliminary results on the performance of BP that tend to investigate its dynamical properties. This constitutes Chap.1.

Chap.2 deals with statistical physics. We expose a few methods that are aimed at providing the partition function, as MF approach, replica method and cavity approach. About the latter, we accurately derive the corresponding equations to highlight their similarity with the BP equations.

In Chap.3 are exposed the main notions concerning the dynamical systems. By the use of famous examples, such as the logistic map and the Lorenz system, we present the classical quantifiers utilized to numerically describe the non-trivial behaviors like chaos.

Part II deals with the detailed development of the message-passing algorithms. Chap.4 contains the foundations of the BP algorithm. We introduce the Bayesian inference on factor graphs, then is presented the Bethe approximation and the equivalence between its stationary points and the fixed-points of the BP algorithm. We end this chapter by the open problems linked to the loop-like topology of the Tanner graph.

In Chap.5 is presented the region-based approximation. We demonstrate the equations of the GBP algorithm by the minimization of the region-based free energy. After describing the introduction of a damping factor in the GBP equations to maintain its convergence, we expose the basic construction rules of the region-graph with few examples on small factor graphs. After that is detailed our original constructions based on the specific topology of a selected LDPC code, then we exhibit encouraging experimental results.

In Chap.6 is described a toolbox to investigate the dynamics of the BP and the GBP algorithms for the LDPC codes. Presented as dynamical systems, we define their mathematical environment and the classical quantifiers, the bifurcation diagram and the Lyapunov exponents, among others. We then expose original tools to relate dynamics and decoding performance of the algorithms, allowing us to confront them.

Introduction

In this first part, we tackle a large overview of different algorithms that are aimed at analyzing and solving factor graphs. We describe the issue of this work according to several research areas, given that this thesis deals with a complex multifaceted investigation. As a matter of fact, inference problem on graphical models represents a challenge in error-correcting codes theory, statistical physics, Bayesian inference, nonlinear optimization and discrete dynamical systems theory. Each of these domains helps construct an exceptional research field where it is crucial and equally hard to merge the corresponding expertises. The work that we have carried out in this thesis needs an introduction in all these fields to well-understand the in-depth-study we will present. To this end, we describe in a few pages essential notions that helps the reader to get accustomed to these theories and the associated notations. It also allows us to focus on open problems that we will deal with all along the manuscript.

The methods presented in this part are not original, as many articles and theses have already dealt with their foundations and their properties. Thus, we do not exhibit many experimental results but the main papers and their significant contributions to the advances in the corresponding topics. In addition, we make an effort to introduce connections between different chapters in order to make the whole presentation consistent, given that problems and solutions emerging from any research field is not independent from the other ones. Moreover, we attempt to tackle each theory with the same global vision to enlighten noteworthy points, useful in our work. Accordingly, we do not describe all relevant aspects of each theory but only the ones that we will make use of in the manuscript, that still represent a quite substantial work.

This part is organized as follows. In first chapter, we introduce the general engineering context in which the thesis began, and the corresponding raised issues. We take the time to present a few decoding algorithms that attempt to solve inference on Tanner graphs. The second chapter is dedicated to an overview of statistical physics, especially methods that help compute the partition function -to some extent the free energy -of any spin glass. It includes the main equations of the mean-field approximation and the cavity method that both enable us to connect with the equations of iterative decoding algorithms. The third chapter deals with a description of the tools and methods that we generally need to describe any dynamical system. We explain the main notions about chaos, stability and bifurcations by use of typical examples, like the Lorenz system and the logistic map. We also exhibit calculations of the Lyapunov exponent that help numerically analyze chaos, which will prove to be very helpful in the second part of the manuscript.

-Chapter 1 -

Context: channel coding § 1.1 Introduction

In this chapter we present the engineering context and the motivations that inspired this thesis. To this end, we first provide few basic definitions about channel coding and Low-Density Parity-Check (LDPC) codes. Afterwards, we focus on the main decoding algorithms derived from the Belief Propagation (BP) that have been published in the literature. We exhibit brief results concerning their relative performance. We also relate the iterative decoding with the optimization of the free energy in statistical physics, that will help us focus on drawbacks of BP. We mention the main articles that made the research progress. We end this chapter considering BP as a dynamical system, and we give the main publications that have proposed first investigations on this topic. § 1.2 Digital communications model

We consider a couple made with an emitter and a receiver. The emitter sends a signal, mapped to a sequence u = [u 1 , . . . , u κ ] of κ bits, to the receiver. In the case where no disturbance affects the signal the receiver would obtain the very sequence u. However, in real life, the "space" between both users is non-empty in the sense that many other signals are traveling inside it. These signals interfere with the one we are working on. Furthermore, space itself is disruptive due to its non-homogeneous constitution, therefore it is as well a source of noise. The consequence is that the emitted signal is strongly disturbed. The raised issue is then how to faithfully recover it. We call the space between the emitter and the receiver the channel. A channel can be modeled as a function ch(u) that maps the input sequence u to an output sequence y = [y 1 , . . . , y κ ] which values depend on the nature of the channel, see Fig. The channel is assumed to produce random noise and makes the output a random reading. Therefore we can define a conditional probability distribution p(y|u) that is called the likelihood. We consider a memoryless channel therefore the likelihood is factorisable as:

p(y|u) = κ ∏ i=1 p(y i |u i ) (1.1)
One of the main channels investigated in the literature are:

• the Binary Symmetric Channel (BSC) of crossover probability p. Any input bit u i is flipped by the channel with a probability p. We represent this channel by the diagram in Fig. 1.2(a). The likelihood of the bit u i is:

p(y i = 0|u i = 0) = p(y i = 1|u i = 1) = 1 -p (1.2) p(y i = 0|u i = 1) = p(y i = 1|u i = 0) = p (1.3)
• the Additive White Gaussian Noise Channel (AWGNC) of power σ 2 . Any input bit u i is disturbed by a random element n i that follows a Gaussian law whose variance is σ 2 , i.e., n i ∼ N (0, σ). We represent this channel by the diagram in Fig. 1.2(b). The likelihood of any bit u i is given by The aim of the receiver is to estimate u by a sequence û using the likelihoods and metrics like the Hamming distance d H (u, û). Defined as the number of bits that differ between both sequences, it is aimed at being as low as possible. § 1. 3

p i (y i |u i ) = 1 √ 2πσ 2 e - (y i -u i ) 2 2σ 2 (1.4) 1 0 1 0 p 1 -p (a) Binary symmetric channel u ∈ F κ 2 n ∼ N (0, σ) y ∈ R + (b) Additive white Gaussian noise channel

Naive estimation

The input sequence u is estimated by scanning all possible states of u and extract the most likely one according to the a posteriori probability distribution p(u, y). This distribution is computed according to the Bayes rule:

p(u, y) = p(y|u)p(u) (1.5)
As every word u ∈ F κ 2 has the same probability to be sent, then û is:

û = arg max u p(y|u) (1.6)
The memoryless assumption allows us to reduce the search from 2 κ states to a search over two states for each bit, i.e. 2κ states. In other words, maximizing the joint probability distribution is equivalent to maximizing each marginal likelihood:

û = [ û1 , . . . , ûκ ] (1.7) ∀i ∈ {1, . . . , κ}, ûi = arg max u i p i (y i |u i ) (1.8)
We display in Fig. 1.3 the performance of this method with κ = 200. The commonly used estimator to judge about the accuracy of any transmission is the Bit-Error Rate (BER), i.e., the ratio between the number of wrong estimated bits and κ. This transmission method is not very efficient given that the BER is too far away from Shannon's limit. Practically we need the information to be more robust against the channel noise to ensure a reliable transmission. § 1. 4

Channel coding

To improve the performance of a transmission, error-correcting codes (ECC) are used. An ECC makes the input signal more robust against the channel noise basically by encoding it. The encoding consists in adding redundancy bits that the receiver uses to decode. The decoding process is aimed at correcting errors caused by the channel. The model is presented in Fig. 1 The dashed box frames the channel coding part, i.e., the part that we will tackle in this manuscript, particularly the box D. We do not focus on the design and the inversion of G that are part of a related but different research area in coding theory.

Encoding

In this work, we investigate the linear ECC, whose encoding blocks are binary matrices, denoted by G, that map u ∈ F κ 2 to x ∈ F N 2 with N > κ. Each coordinate x i is a linear combination of few elements of u. We define the ECC, called simply the code, as the set C = Im(G), each vector of this set is called a codeword. For example the Hamming encoding matrix [R.W50] is:

G =           1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 1 1 1 0 1 1 0 1 1 0 1 1           whose associated code is a set of 2 4 = 16 codewords C = {x = [x 1 . . . x 7 ]}
such that for any of these codewords:

x 1 = u 1 x 2 = u 2 x 3 = u 3 x 4 = u 4 x 5 = u 1 ⊕ u 2 ⊕ u 3 x 6 = u 1 ⊕ u 2 ⊕ u 4 x 7 = u 1 ⊕ u 3 ⊕ u 4 1.4. Channel coding
with ⊕ the binary addition. We define the parity-check matrix H of a code as Im(G) = Ker(H). The advantage of this matrix is that it maps the bits x = [x 1 , . . . , x N ] that are the elements we work on in the channel coding. It is easier to check if a vector x is a codeword by comparing H.x with 0 rather than searching for a vector u such that x = G.u. For example the Hamming parity-check matrix is:

H =   1 1 1 0 1 0 0 1 1 0 1 0 1 0 1 0 1 1 0 0 1  
thus any codeword x ∈ C verify the following equalities:

x 1 ⊕ x 2 ⊕ x 3 ⊕ x 5 = 0 x 1 ⊕ x 2 ⊕ x 4 ⊕ x 6 = 0 x 1 ⊕ x 3 ⊕ x 4 ⊕ x 7 = 0
The most famous codes are the Turbo-codes [C. 93] and the Low-Density Parity-Check (LDPC) codes [R.G63]. The work we will present here is oriented to the latter ones. Any row H i of H applied on a vector x is a binary sum of few elements of x:

H i .x = x j ⊕ x k ⊕ x l ⊕ . . . (1.9) with j, k, l, • • • ∈ {1, . . . , N}.
This sum equals zero if x is a codeword, and one otherwise. That is why the M rows of H are called parity-check equations.

In the example of the Hamming code, M is smaller than N which is actually a common property of ECC. Indeed, using an ECC makes the number N of bits to transmit x = [x 1 , . . . , x N ] larger than the number κ of original bits u = [u 1 , . . . , u κ ]. Adding parity-check equations provides a strong robustness against the channel noise, but it also dramatically reduces the information flow. Using more parity-check equations than the number of transmission bits is not reasonable, therefore we always consider that M < N, and most times, κ = N -M.

Maximum Likelihood Decoding

Using the ECC that encoded the information, the receiver is now able to recover the word x as faithfully as possible, which allows it to recover the information word u. To this end it uses a decoding algorithm. A decoding algorithm or decoder consists in finding the codeword x ∈ C that is the closest one to the vector x, according to the Hamming distance, given the knowledge of H and y. The Maximum Likelihood Decoder (MLD) is said to be optimal in the sense that it offers the lowest Hamming distance. This algorithm consists in spanning through the whole code to extract the codeword x that maximizes the likelihood: ∀x ∈ C\x, p(y|x) ≤ p(y|x) (1.10)

In Fig. 1.5 is presented the performance of the MLD applied on the Hamming code. Even though this algorithm is optimal, it involves a really large computation time given that the code contains 2 κ codewords. Furthermore, in practice the code C is not systematically available, i.e. the receiver only has H at its disposal. In this case, it is unavoidable to scan the 2 N sequences of F N 2 to find the optimal x, given that N can reach very large values. Finally this optimal decoding is not well-suited because of its computation time. § 1.5 Practical decoding algorithms Any error-correcting code has to be associated with a decoding algorithm, also called decoder, to retrieve the information bits that were corrupted by the channel noise. A few codes are associated with specific decoders, i.e. algorithms that are either built especially according to their own topology, e.g. the polar codes [E. 09], or built with dedicated algebraic rules e.g. the Reed-Solomon codes [I.S60] and the non-binary BCH codes [A. 59]

,[R.C60].
As mentioned in the previous paragraph, the optimal decoder for any LDPC code is the MLD but it has no practical interest therefore it is replaced by iterative decoders of dramatic lower complexity.

Graphical representation

Any ECC that is represented by its parity-check matrix has a graphical representation, called the Tanner graph [F.R01]. It is a bipartite graph whose vertices are:

• variable nodes {X i } i • (one per bit),
• check nodes {c a } a (one per parity-check equation).

To lighten notations, we denote by X i (resp. c a ) the corresponding bit and the variable node (resp. the corresponding parity-check equation and the check node). The value x i is said to be the state of X i . We draw an edge e ia between nodes X i and c a if and only if X i is an argument of the parity-check equations c a .

And we say that c a (resp. X i ) is a neighbor of X i (resp. c a ), i.e., c a ∈ N i (resp. X i ∈ X a ). Usually, we represent the Tanner graph in a 2-level mode, e.g. the Tanner graph of the Hamming code represented in Fig. 1 The G-B is a message-passing algorithm where each edge of the Tanner graph carries two directed messages. A message from a variable node to a check node is the state that this variable node should take according to all the check nodes it is related to, except the destination. A message from a check node to a variable node is the XOR between all the incoming messages except the one from the destination. Actually, any message is valued in {0, 1}. The BF appears to be the most natural algorithm, after the MLD: knowing that the aim of an LDPC code is to obtain all its parity-check equations verified, it consists in flipping the bits whose neighboring check nodes are mainly unsatisfied until all the check nodes are verified.

.6. X 1 X 2 X 3 X 4 X 5 X 6 X 7
Such an algorithm is not strictly speaking a message-passing algorithm, but thinking that the flip of a variable node X i results from the information collected on its neighborhood, it could be regarded as messages transmissions.

In the case of a BSC, the BF and the G-B have the clear advantage to be very light according to the hardware requirements, and of quite low complexity given that the update rules are very simple. We display in Fig. 1.7(a) the necessary number of iterations K to obtain a valid answer from both decoders, i.e. when the estimate of any of them is either verified by all the check nodes or identical to the estimate of the previous iteration. We observe that the G-B turns out to be significantly faster than the BF. In Fig. 1.7(a) we display the BER of the two algorithms. We clearly conclude that the G-B is also more accurate than the BF. As an example, for p ≈ 10 -2 , G-B yields a hundred times lower error rate than BF with only four times less computing iterations.

Belief Propagation

Diversity of transmission channels made research being more active about decoders given that some channels result in soft output values. Hard decoders are not able to deal with them, e.g. AWGNC, Rayleigh channel etc. In addition, even for hard valued channels, see results on Fig. 1.7, generalizations of the G-B worth studying in order to improve it. Extending the messages alphabet should provides significant improvements.

Considering no hardware implementation, we can make use of an infinite alphabet by assuming that any message from a node a to a node b in the Tanner graph is the probability distribution of b conditioned by a. This assumption makes the decoder being not a hard decoder anymore but a soft decoder. Such a model was proposed in [START_REF] Pearl | Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF] as a solution to the inference problem on Bayesian trees and polytrees. It provides marginal probability distributions, or marginals, of all random variables that build any tree. Extracting arguments of the maximum of all marginal distributions then results in the most likely state. This model is the Belief Propagation (BP) that constitutes the basis of many other iterative algorithms as it was pointed out in

[F.R98],[D. 95],[R.J98].
A decoder that is largely used in practical cases and that stems from BP is the Viterbi algorithm [A. 67]. It was created for decoding convolutional codes and Turbo-codes. It does not extract approximates of marginals but it provides the most likely sequence of bits according to a graphical representation of the code, the trellis, that is not a Bayesian network but a hidden Markov chain. The Viterbi algorithm is also known as the Max-Product version of BP closely related to dynamics programming [D.P01]. Also used to decode Turbo-codes and convolutional codes, the BCJR decoder [F.R01], often called Forward-Backward, focuses on the minimization of the bit-error rate whereas the Viterbi is aimed at minimizing the word-error rate.

The Sum-Product Algorithm (SPA) [F.R01] is the name of BP written especially for factor graphs e.g. Tanner graphs of LDPC codes. Originally built for pairwise Bayesian networks that are not bipartite graphs, we only find one update rule. The message from the variable node X i to the variable node X j , displayed on Fig. 1.8, is defined as:

m ij (x j ) = ∑ x i ψ ij (x i , x j )φ i (x i , y i ) ∏ X k ∈N i \X j m ki (x i ) (1.11)
where:

• ψ ij (x i , x j ), the compatibility function, indicates if the values x i and x j of the variable nodes X i and X j are valid. In the context of LDPC codes, ψ ij (x i , x j ) = 1 ⊕ x i ⊕ x j is the parity-check equations that links X i with X j ,

• φ i (x i , y i ), the evidence, provides the observation we obtain from the node X i . It is the likelihood p i (y i |x i ).

X i X j N i Figure 1.8: Message m ij in a Bayesian network
This update equation is a sum-product rule, given that it contains both basic operations. Qualitatively, sum operand allows or not to propagate information coming to X i (except information from X j ), and product operand provides resulting information. This update rule has been extended [F.R01],[J.S05] to graphs that are neither pairwise nor Bayesian networks but multi-body interaction undirected graphs as Markov Random Fields (MRF) or factor graphs G = (V, E). These models are bipartite, i.e. the set of nodes V is a union of two sets of nodes of different kinds X ∪ C. Multi-body interactions between nodes of X are not pairwise then they cannot be implicitly represented by edges, they need their own nodes C. An edge e ∈ E represents a contribution of a given variable of X to a given multi-body interaction in C. Two update rules, that gave birth to the SPA, were defined according to Fig. 1.9: from variable nodes to check nodes:

∀e ia ∈ E, ∀x i ∈ {0, 1}, m (k) ia (x i , y) = l i (x i , y) Z ia ∏ c b ∈N i \c a n (k-1) bi (x i , y) (1.12)
from check nodes to variable nodes:

∀e ia ∈ E, ∀x i ∈ {0, 1}, n (k) ai (x i , y) = ∑ x a ∪x i f a (x a ) ∏ X j ∈X a \X i m (k) ja (x j , y) (1.13)
where:

• k is the iteration index,

• l i (x i , y) = p(y i |x i ) is the likelihood of the variable X i in the state x i , • Z ia ensures that m (k) ia (0, y) + m (k) ia (1, y) = 1, • x a = X i ∈X a x i , • f a (x a ) = 1 ⊕ X i ∈X a x i . X i c a N i \c a n ia
(a) Message n ia from the variable node X i to the check node c a with ia (x i , y) = l i (x i , y). Statistically, the initialization choice does not wield much influence on the outcome of the algorithm, therefore it is up to the experimenter. BP provides the probability distribution {b i (x i , y)} x i ,y of each variable node X i , called the belief, according to the following equation:

N i = {c b ∈ C | (c b , X i ) ∈ E} X i c a X a \X i m ai (b) Message m ai from the check node c a to the variable node X i with X a = {X j ∈ X | (X j , c a ) ∈ E}
∀X i ∈ X, ∀x i ∈ {0, 1}, b (k) i (x i , y) = l i (x i , y) Z i ∏ c a ∈N i n (k) ai (x i , y) (1.14)
where

Z i = b (k) i (0, y) + b (k)
i (1, y) ensures that the belief builds a true distribution. Finally a codeword estimate is computed as the concatenation of the most likely states of these distributions:

x(k) = X i arg max x i b (k) i (x i , y) (1.15)
SPA is stopped whether H.x = 0, i.e. it results in a codeword, or for any edge e ia , the message m

(k) ia (x i , y) (or n (k)
ai (x i , y)) does not change from an iteration to the next one, i.e. SPA has converged.

Alternatives to the BP

Even though the complexity of the SPA was very low compared with the MLD one, it turned out necessary to simplify the update equations for specific decoding applications. One of the solutions consists in working in the logdomain. It leads us to the logarithm version of the SPA, partly introduced in [R.G63], where a vector message [m ai (0) m ai (1)] T becomes a scalar message M ai = log m ai (0) m ai (1) . The update rule is then: 16) where N ja = log n ja (0) n ja (1) . Computation complexity is dramatically decreased, and computation memory is twice less given that only one quantity is used for any message and for any belief. However the applications tanh and tanh -1 are not easily handled. To circumvent this undesired effect of the log-domain we could store an accurate tab of these functions but it would need a nonnegligible memory. The Min-Sum (MS) algorithm, an approximation to the log-SPA [J. 96], was defined such that a message from check node c a to variable node X i contains only basic operations:

M ai = 2 tanh -1   ∏ X j ∈X a \X i tanh N ja 2   (1.
M ai =   ∏ X j ∈X a \X i sign(N ja )   min X j ∈X a \X i {|N ja |} (1.17)
that substantially reduced computation time, unfortunately at the expense of the performance.

Last generic decoders we have to mention are the Finite Alphabet Iterative Decoders (FAID) also called multilevel decoders presented in [S.K10],[S.K11]. They consist in modifying the MS algorithm by changing the update rule of the messages from the variable nodes to the check nodes {n ia } i,a , each of these messages being computed according to nonlinear equations. Studies tend to prove that these decoders provide lower BER than the SPA on the BSC for very low values of the crossover probability of the channel, that is a non-negligible progress [D. 12]. Unfortunately, FAID are somehow hard decoders that are not easy to use on AWGNC, furthermore they cannot be used on any LDPC code. We display on Fig. 1.10 performance of hard and soft decoders on BSC and AWGNC. BP sparked a lot of interest among different research areas: Bayesian inference, channel coding and also nonlinear programming. In [J.S03a],[J.S02] is described a physical aspect of the decoding algorithms for spin glasses [START_REF] Nishimori | Statistical physics of spin glasses and information processing an introduction[END_REF] and more generally for factor graphs, e.g. the Tanner graphs. It turns out that the inference problem is equivalent to the minimization of the free energy F of the graph. The free energy is computed from the Boltzmann's joint distribution p(x):

p(x) = e -E(x) Z (1.18) F = -log Z (1.19)
where Z is the partition function that ensures the normalization of the distribution, Z = ∑ x e -E(x) . Given that computing Z consists in spanning through the 2 N possible states of the graph, we conclude that it is not tractable. Instead of the free energy, variational free energies [H.A35],[R. 51] computed by tractable factorized joint distributions prove to be relevant approximations.

Therefore the main issue is to find the factorized distribution that minimizes the variational free energy. This problem can be solved by BP. Indeed, it was shown in [J.S05],[P. 03] that fixed points of BP are equivalent to stationary points of the variational Bethe free energy F β . However, the update rules of BP are not exactly established according to the constrained Lagrangian optimization of F β , the equivalence is only true for steady states, that are not always encountered. In addition the loop-like topology of most factor graphs prevents the minimization from being optimal, i.e. either the result is only a local minimum, or the minimization completely fails in the sense that it oscillates between two or more states or even it falls into chaotic behaviors. Actually, the barrier is the fact that F β is very rarely a convex functions of the beliefs, making its minimization a very hard task i.e. lack of convergence of BP is an open issue.

The optimization of F β is dealt in [T. 03a],[T. 04],[T. 06a] by extracting conditions on the graph for the existence and the uniqueness of the stationary points, e.g. the BP yields exact minimum of F β for graphs that contain only one single loop. In spite of this encouraging property, this result makes the investigation on the BP more complex because many graphs contains much more than one single loop. An alternative version of the BP is proposed in [A.L03] to ensure the convergence. This ConCave-Convex Procedure (CCCP) consists in splitting F β into a convex function F vex and a concave function F cav in order to iteratively find, as in Newton's method, the point that cancels the gradient of the variational free energy:

∇F vex ({b (k+1) i } i ) = -∇F cav ({b (k) i } i ) =⇒ ∇F β ({b (k→K) i } i ) = 0 (1.20)
with K the large number of iterations. The algorithm is made of two nested loops:

• the inner loop on the Lagrange multipliers ensures that the linear constraints are verified (iteration index: k),

• the outer loop provides the values of the beliefs (iteration index: τ).

On one hand, update equations of the multipliers involve the beliefs, on the other hand, the update equations of the beliefs involve the multipliers. Then depending on the initialization of the algorithm, computation time can be very large compared to the original BP one. This property often makes this algorithm unpractical. However, messages update rules are strictly consistent with the Bethe free energy minimization, it even appears clear that the messages are none other than the Lagrange multipliers themselves. Thus the CCCP is of very high interest to understand a few subtleties of BP. We display on Fig. 1.11(a) the BER of the CCCP and the BP algorithms on the Tanner code. We observe that the accuracy of the CCCP is slightly improved compared with the BP, as the SNR is increased. The drawback of this method is depicted in Fig. 1.11(b) where it appears that the number K of loops (outer loops for the CCCP) is almost doubled by the use of the CCCP compared with the BP. Furthermore, according to Fig. 1.12, for low SNR values, the number of inner loops τ is non-negligible. For instance, for an SNR = 3.2 dB, we obtain τ = 77 that is very high. Even though this value decreases as k increases, computation time is highly lengthened by the CCCP method. Knowledge of the average BER is not sufficient to understand the whole BP algorithm, e.g. properties about stability, convergence, fluctuations around the BER, etc. For instance, we represent on Fig. 1.13 the convergence rate and the required number of iterations to converge, both averaged over ten thousand simulations. One would a priori think that, as the BER, these quantifiers should decrease according to the SNR. However, it appears an unexpected behavior of BP around SNR = 1.0dB: the convergence rate significantly decreases and the number of iterations required to have the BP converge considerably increases. Between 2.0dB and 3.0dB both quantities recover their expected behavior relatively to their initial values at 0.0dB. As mentioned in many papers, the loop-like topology of most codes could be responsible for this undesired phenomenon.

Another factor is the nonlinearity of the update rule that enables us to compute the messages from the check nodes to the variable nodes (1.13), (1. 16), (1.17). According to research results on dynamical systems, nonlinearities are responsible for non-trivial or even unexpected phenomena [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF] e.g. oscillations, aperiodicity or even chaos.

Studies of the nonlinearities of the BP has been well introduced for Turbocodes [D. 00],[L. 02], with few improvements that help the decoding to be more accurate. It was introduced measurements to quantify nonlinear effects. The bifurcation diagram revealed that SNR values could be shared into intervals according to the BP behavior, and the largest Lyapunov exponent helped confirm the chaotic behavior of the decoder for particular SNR values. Concerning LDPC codes, a less extensive work was presented in the sequel [L. 06]: only pseudo-trajectories without any bifurcation diagram or Lyapunov spectrum that could exhibit BP chaotic behaviors. In [X. 05] is presented a dense work on BP decoding using the mean square beliefs: Research on decoding algorithms is confronted to the problem of suboptimality of BP-based iterative decoders. Optimizing the variational free energy, taken as the cost function of any message-passing algorithm, is a very hard task given that this function is almost never convex. The CCCP appears a good candidate to overcome this problem, but we saw that its computation time makes it unrealistic. Thus, investigating on the free energy is necessary to understand its importance in the decoding algorithms. In this purpose, we describe, in the following of the manuscript, basic notions of statistical physics of spin glasses. Furthermore, preliminary results that we exhibited about BP convergence allows us to deeply investigate this decoder with the tools provided by the dynamical systems theory. Therefore we need to introduce basic notions and tools of this theory, in order to use them to highlight the dynamical properties of the decoder.

E(k) = 1 N N ∑ i=1 b (k) i (x i ) 2 (1.
-Chapter 2 -

Statistical physics of spin glasses § 2.1 Introduction

In this chapter we describe physics behind the work presented in this thesis. Decoding algorithms have deep relations with methods of statistical physics of spin glasses. In order to keep a reasonable manuscript size, we will not go deeply into the following theories. The purpose is to give an overview of the tools used in our research.

In section 2.2 are exposed the main principles of statistical physics that are necessary to develop the decoding algorithms for error-correcting codes. To this end, we explain in section 2.3 the notion of spin glasses and the associated free energy. In sections 2.4, 2.5, 2.6 and 2.7 we describe methods used to compute or estimate the free energy. § 2.2 Problems

In statistical physics, one of the main goals is to estimate thermodynamical functions of any system at equilibrium. Among them, the average energy U and the entropy S help us compute the Helmholtz free energy F:

F = U -TS (2.1)
where T is the temperature set arbitrarily to 1 in the whole study as it simply sets a scale for the units in which one measures energy. This quantity is crucial because it provides knowledge of many other state functions and it is notably related to the partition function Z as:

F = -ln Z (2.2)
where Z provides the exact probability distribution of a physical network S = {S i } i , e.g. spin glasses, given by the Boltzmann's law:

p(S = s) = e -E(s) Z (2.3)
where E(s) is the energy function. However, computing F is a very hard task as physical systems are often very large and interactions between elements of the system are very complex. More precisely there is no global solution but using equation (2.1), i.e researchers have no choice but to study all systems one by one.

Fortunately, physical systems are divided up into classes, as Sherrington-Kirkpatrick model, Edwards-Anderson model, etc. Each class has been studied to exhibit its own solving method. Among them we focus on the spin glasses because they are in a one-to-one match with graphs used in errorcorrecting coding theory. § 2.3 Spin glasses From [START_REF]Decoherence and Entropy in Complex Systems, Selected Lectures from DICE 2002 -Spin glasses: Still Complex After All These Years[END_REF],[M. 87],[B. 00], spin glasses are defined as systems of atoms with localized magnetic moments, namely spins, whose interactions are characterized by quenched randomness: a given pair of spins have a roughly equal a priori probability of having a ferromagnetic or an anti-ferromagnetic interaction, see Fig. The prototype material is a dilute magnetic alloy, with a small amount of magnetic impurity randomly substituted into the lattice of a nonmagnetic metallic host. The signature feature of spin glass behavior is the phase transition, i.e. a dramatical change in the physical magnitudes of the system. The main interesting property of this kind of network, for our work, is that spin glasses suffer from frustration: no spin configuration can simultaneously satisfy all couplings. Such a phenomenon is observed also in channel coding when not all parity-check equations of an error-correcting code are satisfied because the channel has changed the values of the bits.

According to [F. 12], the coupling J ij between two spins s i , s j around a magnetic impurity has random sign and intensity because the distance between two spins is a random variable. The simplest idealization of the interaction between two spins s i and s j is a coupling term -J ij s i s j in the energy function E(s) with J ij a random variable. The disorder of the material comes from the random values of the couplings {J ij } <i,j> . Due to the time independent valuation of these terms, the disorder is said quenched. The spin glass energy function without any external magnetic field is: As explained previously, frustration is a phenomenon that emerges because of the randomness of coupling constants: any spin is subjected to fields due to its neighbors that have different signs. Some want it to point up and others to point down then it cannot be set to a unique state. We represent this event in Fig. 2.2.

E(s) = -∑ <i,j> J ij s i s j (2.
The first figure shows a spin glass where one interaction has been modified such that it goes from anti-ferromagnetic to magnetic. Then on the other figures are depicted the route to the frustration: step by step, the spins flip (in red) to make the interactions, either ferromagnetic or anti-ferromagnetic, true. We then observe that the left blue spin should flip to be consistent with its left and right neighbors but in this case it prevents the interaction with its bottom neighbor from being true. Practically, if the spin flips then we again observe a step by step route to frustration that leads to the same situation. This is very close to the problem of loop-like topology in the case of the LDPC codes.

Frustration makes the search for the ground state not trivial, many degenerate ground states, or metastable states, emerge. Actually, a challenge of the physics is to find the state ŝ that provides the lowest value of the energy function E(s). If all the J ij 's were of same sign then the solution would appear very easily. However, in most cases, systems are paramagnetic the coupling constants have different signs. The naive though optimal method would consist in computing the energy function for all spins configurations, then we extract the lowest energy state. Unfortunately, if it could be done for low dimensional systems, it is completely intractable for high dimensional systems, given that 2 N states have to be computed (N being the number of spins). One would notice that this is actually closely related to the MLD that is also intractable for LDPC codes for the same reason. We here after present an alternative solution based on the Mean-Field (MF) approximation. Its relevance surely comes from the fact that it successfully attempts to represent some important features of the physical spin glass systems of great interest for their peculiar properties, at least at the level of the mean field approximation. A relevant summary of the MF with theory and results is exposed in [F. 06].

From [START_REF]Introduction to the Functional Renormalization Group[END_REF], [START_REF] Nishimori | Statistical physics of spin glasses and information processing an introduction[END_REF], the idea consists in developing the term J ∑ <i,j> s i s j in the expression of the energy function, where J is the coupling assumed to be invariant on the couples of spins. We consequently use the mean spin value:

m i =< s i > . (2.5)
Its value can be derived as:

m i = p i (S i = +1) × (+1) + p i (S i = -1) × (-1) (2.6) = 2p i (S i = +1) -1 (2.7)
Given that for each spin S i the state is s i = m i + δ i , where δ i is the fluctuation of the spin around its mean value, the cross product between two neighboring spins is:

s i s j = m i m j + m i δ j + m j δ i + δ i δ j (2.8)
The MF approximation is to neglect the quadratic fluctuations δ i δ j , considering small enough in particular for high dimensional systems, i.e.:

s i s j ≈ m i m j + m i δ j + m j δ i (2.9)
Furthermore the mean value m i of a spin s i is expected to be site-independent, i.e. for any spin s i , < s i >= m. This yields:

s i s j ≈ m 2 + m(δ j + δ i ) (2.10)
Using the definition of the fluctuation δ i = s im, it appears that:

s i s j ≈ -m 2 + m(s i + s j ) (2.11)
As the spins are valued in {-1, +1} their sum is either 0 or 2s i . Thus the remaining term is:

s i s j ≈ -m 2 + 2ms i (2.12)
The energy function is written such that all the couplings are considered equivalent:

E(s) = -J ∑ <i,j> s i s j (2.13)
Under the MF approximation it becomes:

E(s) = JN c m 2 -Jzm ∑ i s i (2.14)
where N c is the number of interacting couples and z is said to be the coordination number. Considering an external magnetic field interacting with the spins through the energy E ext = -h ∑ i s i we obtain the total MF approximation:

E(s) = JN c m 2 -(h + Jzm) ∑ i s i (2.15)
The whole system has then been decoupled into a sum of one-body system with an effective mean-field h e f f = h + Jzm, which is the sum of the external field h and of the mean-field induced by the neighboring spins. Using the Boltzmann's law it appears that the MF joint distribution is:

p MF (x) = 1 Z e -JN c m 2 N ∏ i=1 e h e f f s i (2.16)
Finally we can write that the distribution is a factorization along the spins:

p MF (x) = 1 Z MF N ∏ i=1 p MF,i (x i )
(2.17)

The MF approximation provides analytic expressions of physical quantities as the partition function and the free energy, i.e. it solves the spin glasses.

Message-Passing

The MF approximation provides an algorithm to obtain the {m i } i . Given (2. 16) and (2.17) we can write the MF marginal on the spin S i as:

p i (s i ) = 1 Z i e h+J ∑ S j ∈N i s j (2.18)
Then the magnetization of the spin S i is according to (2.6):

m i = tanh(h + ∑ S j ∈N i m j ) (2.19)
A simple approach to obtain {m i } i is iteratively updating them according to the following equation:

m (k+1) i = tanh(h + ∑ S j ∈N i m (k) j ) (2.20)
This message-passing results after K iteration in the computation of the marginals given that:

p i (s i = +1) = 1 + m (K) i 2 , p i (s i = -1) = 1 -m (K) i 2 (2.21)
Actually the MF expression of m i is not as accurate as we would expect. Furthermore when averaged over all the configurations of the spin glass, < m i > is not easily handled because of the nonlinearity of the tanh function

< tanh(h + ∑ S j ∈N i m j ) > = tanh(h + ∑ S j ∈N i < m j >)
that has been investigated, particularly in the TAP method proposed by Thouless, Anderson and Palmer [START_REF]Advanced Mean Field Methods[END_REF]. § 2.5 Replica method

As mentioned in the introduction of the chapter, statistical physics is aimed at computing the free energy F =ln Z. It turns out that computing integer moments of the partition function Z n with n ∈ N is much easier than computing the log-partition function. Fortunately, mathematics offer a way to use this alternative by the following identity:

ln Z = lim n→0 Z n -1 n (2.22)
Now we DESCRIBE, according to [M. 06], [START_REF]Advanced Mean Field Methods[END_REF], the basics of the replica method based on this identity. Z n is a product of n identical partition functions. We consider a system S of N spins, represented by an energy function E(s), that is independently duplicated into n replicas. The overall distribution of these n systems is:

p(S 1 = s 1 , . . . , S n = s n ) = p 1 (S 1 = s 1 ) . . . p n (S n = s n ) (2.23)
where the distribution of a system S i verifies the Boltzmann's law, i.e.:

p(S 1 = s 1 , . . . , S n = s n ) = 1 Z 1 e -E(s 1 ) . . . 1 Z n e -E(s n ) (2.24) = e -(E(s 1 )+•••+E(s n )) Z 1 . . . Z n (2.25)
We obtain a Boltzmann distribution with an energy function linearly decomposed over the different replicas of the system:

Ẽ(s 1 , . . . , s n ) = E(s 1 ) + • • • + E(s n ) (2.26)
and the partition function is:

Z = ∑ s 1 • • • ∑ s n e -(E(s 1 )+•••+E(s n )) (2.27)
For more details on this method that sparked off much interest, we refer the reader to the excellent reviews in [START_REF] Nishimori | Statistical physics of spin glasses and information processing an introduction[END_REF] and [M. 06]. § 2.6 Cavity approach

The cavity method [F. 11] relies on the computation of the partition function Z by an iterative process. We will see that such an algorithm is completely equivalent to the BP as it propagates messages between nodes of the graph.

That is why we will describe it in a more detailed way than the other methods. The idea is to compute the partition function recursively, by considering a dynamical graph: at each step is added a variable node to the graph, the graph being completely recovered at the last step. Equivalently one may think of taking one variable node out of the system and computing the change in the partition function. The name of the method comes exactly from this image: one digs a cavity in the system.

S 1 S 2 S 3 -→ h

Ising chain

As for BP, the original algorithm is defined on trees. Let's consider an Ising chain, i.e. a pairwise tree whose spins are binary variables and of degree ≤ 2, with a small number of spins to develop the equations, see Fig. 2.6. We now derive the equations step by step given that two spins S i , S j interact via the interaction -JS i S j and that a spin S i interacts with the magnetic field via the energy -h i S i .

Step 1: spin S 1

The tree is made with spin S 1 only and field -→ h .

S 1 -→ h Figure 2.4: N = 1
Then the partition function is:

Z (1) = ∑ s 1 Z 1 (s 1 ) (2.28)
with Z 1 (s 1 ) the partial partition function with the spin S 1 fixed to the value s 1 : Z 1 (s 1 ) = e h 1 s 1 (2.29)

Step 2: spins S 1 and S 2

The tree is made with the spin S 1 and the spin S 2 and the field -→ h .

S 1 S 2 -→ h Figure 2.5: N = 2
Then the partition function is:

Z (2) = ∑ s 2 Z 2 (s 2 ) (2.30)
with Z 2 (s 2 ) the partial partition function with the spin S 2 fixed to the value s 2 :

Z 2 (s 2 ) = e h 2 s 2 ∑ s 1 e Js 1 s 2 e h 1 s 1 (2.31)
In this equation we recognize Z 1 (s 1 ) defined earlier. However the context is different given that we added a variable to the system. Therefore Z 1 (s 1 ) is not anymore the partial partition function with S 1 fixed to the value s 1 . In order to make it clear we denote the former Z 1 (s 1 ) by Z 1→2 (s 1 ). Then:

Z 2 (s 2 ) = e h 2 s 2 ∑ s 1 e Js 1 s 2 Z 1→2 (s 1 ) (2.32)
Step 3: spins S 1 , S 2 and S 3

The tree is made with the spin S 1 and the spin S 2 and the field -→ h .

S 1 S 2 S 3 -→ h Figure 2.6: N = 3
The partition function is:

Z (3) = ∑ s 3 Z 3 (s 3 ) (2.33)
with Z 3 (s 3 ) the partial partition function with S 3 fixed to the value s 3 :

Z 3 (s 3 ) = e h 3 s 3 ∑ s 2 e Js 2 s 3 e h 2 s 2 ∑ s 1 e Js 1 s 2 e h 1 s 1 (2.34)
Extending the same previous argument, we now denote the former Z 2 (s 2 ) by Z 2→3 (s 2 ), that provides:

Z 3 (s 3 ) = e h 3 s 3 ∑ s 2 e Js 2 s 3 Z 2→3 (s 2 ) (2.35)

Interpretation

We observe for each step i, the partial partition function Z i (s i ) is a product of information coming from the the magnetic field e h i s i with the information coming from its neighbor S j revealed at the previous step ∑ s j e Js i s j Z j→i (s j )

Generally for any tree model, not necessarily a chain, the partial partition function of any spin S i fixed to the value s i is computed as:

Z i (s i ) = e h i s i ∏ S j ∈N i ∑ s j e Js i s j Z j→i (s j ) (2.36)
where N i is the set of spins connected to S i . The quantity Z j→i (s j ) is the partial partition function of the system without the spin S i with S j fixed to s j . One could verify this way that we obtain that ∑ s i Z i (s i ) is strictly identical to the total partition function Z.

Case of a factor graph

The previous claim seems to only address pairwise trees, but it also applies to any model such that the energy function is of the following form:

E(s) = -∑ i h i s i -∑ a J ∏ S i ∈N a s i (2.37)
where ∑ a is the sum over all the spins interactions, and N a is the set of spins that interact through a common interaction denoted by a. Let's consider the simple example on Fig. 2.7 to develop the cavity method. The partial partition function for the spin S 1 in the state s 1 is:

Z 1 (s 1 ) = e h 1 s 1 ∑ s 2 ∑ s 3 e Js 1 s 2 s 3 Z 2→ a 1 (s 2 )Z 3→ a 1 (s 3 ) (2.38)
where:

Z 2→ a 1 (s 2 ) = e h 2 s 2 (2.39) Z 3→ a 1 (s 3 ) = e h 3 s 3 (2.40) × ∑ s 4 e Js 3 s 4 Z 4→ b 3 (s 4 ) (2.41) × ∑ s 5 ∑ s 6 e Js 3 s 5 s 6 Z 5→ c 3 (s 5 )Z 6→ c 3 (s 6 ) (2.42)
and:

Z 5→ c 3 (s 5 ) = e h 5 s 5 (2.43) Z 6→ c 3 (s 6 ) = e h 6 s 6 (2.44)
The subscript a in Z 2→ a 1 helps understand "from S 2 to S 1 through a". One could easily verify that we get the exact total partition function by computing:

∑ s 1 Z 1 (s 1 )

Bipartite messages

In order to simplify the implementation, above equations could be sorted into two kinds. To understand how, one should have a look to the equations of Z i→ α j (s i ). Since such a quantity is a partial partition function for the system without spin S j , it should not depend on S j . In the previous example, one would check that:

Z 3→ a 1 (s 3 ) = Z 3→ a 2 (s 3 ) (2.45)
In this case, the spins S 1 and S 2 do not play any role in the computation, therefore both quantities should be computed as a single message from the spin S 3 to the interaction, called factor, a. Generally we have for two spins S i , S j interacting through the factor a:

Y i→a (s i ) Z i→ a j (s i ) (2.46) = e h i s i ∏ α∈N i \a ∑ s α 1 ∑ s α 2 • • • ∑ s αm e Js α 1 s α 1 ...s αm ∏ S j ∈N α Y j→α (s j ) (2.47)
where {S α 1 , . . . , S α m } is the set of spins interacting via factor α. The second kind of message that we define is the reverse of Y i→a (s i ), i.e. the message that comes from a to S i . The equations of the previous case provide that:

W a→i (s i ) ∑ s a 1 ∑ s a 2 • • • ∑ s am e Js a 1 s a 2 ...s am ∏ S j ∈N a Y j→a (s j ) (2.48)
Then we can rewrite Y i→a (s i ) as:

Y i→a (s i ) = e h i s i ∏ α∈N i \a W α→i (s i ) (2.49)
In this way it is possible to give the partial partition function for a spin S i in the state s i as:

Z i (s i ) = e h i s i ∏ a∈N i W a→i (s i ) (2.50)
Equations (2.48), (2.49) and (2.50) are message-passing equations, analogous to those BP equations we mentioned in the previous chapter such that:

• W a→i is analogous to the message m a→i from the check node c a to the variable node X i ,

• Y i→a is analogous to the message n i→a from the variable node X i to the check node c a ,

• Z i is analogous to the belief b i on the variable node X i .

It helps define the equivalence between the cavity method and BP. As a notice, the cavity method can also be expressed as a message-passing either between the magnetizations of the spins or the probability distributions of the spins, see [F. 11] for details about this. § 2.7 Mean-Field generalizations

In chapters 4 and 5 we will present improvements of the MF approach, namely the Bethe approximations and the region-based approximations [J.S02],[J.S05].

The common principle of both methods is to generalize the equation (2.17) to the following:

p(x) = 1 Z R ∏ r∈R p c r r (x r ) (2.51)
Instead of factorizing the joint distribution on single variable nodes, we factorize on subsets or regions of variable nodes. This way, we keep, to some extent, the couplings between the variables mentioned in the energy function. These methods are actually part of the global region-based free energy approximation detailed in the chapter 5. § 2.8 Conclusion

The MF approximation is a relevant way to provide tractable computation of the thermodynamical functions. We saw that it can also be expressed as a message-passing algorithm between the spins of any spin glass. Similarly, the cavity approach has been showed to present the same equations as the BP's, highlighting a strong connection between statistical physics and decoding algorithms. As such, this method is also subject to the loop-like harmful effects of the spin glasses, or any other physical factor graphs. To our knowledge, the replica method has not been extended to decoding algorithms, partly due to its quite high complex theoretical content. Other methods that we did not describe here are also investigated in this respect, e.g. the renormalization group [G. 10],[J. 03b].

We will see in the next part that these methods can surely be improved by considering less restrictive decoupling between spins. That is the rationale behind MF generalizations, also known as Bethe and region-based approximations.

-Chapter 3 -

Dynamical systems § 3.1 Introduction

In this chapter, we present a few tools that help study dynamical systems. We cannot describe all the extent of the theory, therefore we will focus on the ones that are to be used on the decoding algorithms that are nonlinear systems due to the nonlinearities involved by the update rules. Linear systems are not surprising as they always give linear response, i.e. their effects are proportional to the causes. On the contrary, evolutions of nonlinear systems turn out to be often unexpected, though their update equations are deterministic. Nonlinearities can make saturate a system or even make it diverge.

Study of so called dynamical systems results from the fact that most systems evolve in time and that they need to be reliably anticipated. The property of stability remains central in the study of any nonlinear system. It is qualitatively for a system to keep a regular and controlled evolution when something is slightly changed in its initial conditions.

Dynamical systems are defined by state variables which time evolutions, or state equations, provide the description of the systems. These evolutions are logically defined by differential equations when variables are continuous or iterative equations when variables are discrete in time. Complex and often unexpected behaviors of nonlinear systems do not appear obvious even though state equations are completely deterministic. As a matter of fact, one could think that the lack of control would come out if the state equations followed some probabilistic laws, because it somehow represents a lack of knowledge. Actually, we will see that strange behaviors do not need such a property but only a series of causes and effects submitted to particular initializations and parameters values. This is partly due to the absence of analytic solution to the state equations but only numerical ones. Consequently, a slight difference in the numerical values of few state variables or parameters can lead to massive effects in the evolutions of the system, such as chaos. In such a case there is no predictability and no control on the system, the evolution seems to be erratic and random. Such a sensitivity is inherent to the system itself. Although the purpose of the investigation is to predict the states of any nonlinear system at any time, chaos prevents us from completely mastering it. Therefore, a crucial study is to bring out the main properties of chaos, particularly the minimal conditions under which it occurs.

In section 3.1, we describe the notion of chaos with two well-known examples, the logistic map and the Lorenz system. Section 3.2 is dedicated to the different scenarios that lead any nonlinear system to a chaotic behavior. In section 3.3, we detail the computation of the Lyapunov exponent that numerically estimates whether any system is chaotic or not. Last section 3.4 deals with the situations where a dynamical system is incomplete, i.e. a few state variables are not available. We briefly present a method that re-constructs the state space. § 3.2 Chaotic attractors It is crucial to understand that chaos is neither random nor nonsense. More accurately, chaos is basically opposed to randomness as a chaotic system is a deterministic system whose time evolution is perfectly described by the state variables and the associate update rule(s). Practical systems are subject to chaos because of our incapacity to analytically solve them and to compute state variables with an infinite precision. More precisely, an infinite precision on data would allow us to forecast any state of a system given its initial state whereas a finite precision implies an uncertainty, or computation error, after some computation time.

From [K.T96], chaos is a specific behavior defined as follows: when the system evolution could be done through numerous motions that are initially quite close in the state space, such that after some time, any two of them seem to have nothing in common except they are enclosed in the same region of the given state space (they are dissipative). The separation between them turns out to be often exponential in time. This phenomenon is an intrinsic property of the chaotic attractors, also called strange attractors. To illustrate this peculiar behavior, we will use two well-known examples: the logistic map that models population growth and Lorenz attractor that models -to some extent -the atmospheric convection.

Logistic map

A common model that represents the evolution {x k } k of a resource-limited population is given by:

x k+1 = ax k (1 -x k ) (3.1)
with a a positive real number and x 0 ∈]0, 1[. It appears that a very small value of x k entails a quite linear growing because 1x k ≈ 1 whereas a large value of x k prevents from a proportional growing. In other words, the maps is nonlinear for non small values of x k . Such a nonlinearity is the source of non trivial behaviors presented hereafter. We restrict the study to the case a = 4 but anyone who is interested in extracting the whole information about this system should investigate on all values of a. The first step of the study is to visualize the sequence {x k } k . A common method is displayed in Fig.

map x k+1 = 4x k (1 -x k )
We observe that no particular pattern emerges, we cannot distinguish any fixed-point or repeated motion. This is also exhibited in Fig. 3.2(a) where a trajectory of the map is represented, with the time evolution for a slightly different initial condition x 0 = x 0 + 10 -6 . At k = 14 appears the beginning of the divergence between both trajectories and as the time increases they evolve completely differently. In Fig. 3.2(b) is represented the time evolution of the separation between these trajectories. For k ≤ 21 it is possible to model it by the function: -α) with:

d(k) = e λ(k
λ ≈ 0.9, α ≈ 21. Coefficient α is a delay and λ is a crucial quantity that provides the intensity of the divergence between trajectories. Given that the exponential function takes precedence over any polynomial function, an exponential evolution of the separation is the most divergent motion that can be encountered. Such a property is an evidence of chaos: two initially close trajectories eventually follow very non similar evolutions. Quantity λ is commonly called the largest Lyapunov exponent, or simply the Lyapunov exponent.

In Fig. 3.2(b) it appears that d(k) is not exponential for k > 21. First of all, it is impossible for the separation to acquire values that would be larger than the state space size itself, therefore d(k) ≤ 1 for all k. After that, the trajectories of any dissipative systems live in a restricted part of the state space, named the attractor. This property leads us to evaluate the attractor size by measuring the maximum of d(k). We will study this more deeply in the last chapter of the manuscript.

Lorenz attractor

In climate science, the Lorenz attractor helps predict long term atmospheric aspects. Once more, slight difference in numerical initial conditions leads, under particular assignments of external parameters, to very different behaviors. The Ordinary Differential Equations (ODE) that model the Lorenz system are:

dx(t) dt = σ(y(t) -x(t)), (3.2
)

dy(t) dt = rx(t) -y(t) -x(t)z(t), (3.3
)

dz(t) dt = x(t)y(t) -bz(t), (3.4)
where x(t), y(t), z(t) are the state variables. This set of values traps the system into a non trivial attractor as depicted in Fig. 3.3. This attractor is the chaotic Lorenz attractor. The sensitivity to initial conditions is very high, see Fig. 3.4(a). The second trajectory is initially 10 -6 away from the first one at the Euclidean sense.

-20 -10 0 10 20 -30 -15 0 15 0 15 30 45 The full circles are the initial points and the full squares are the final points. Even though circles are very close, squares are quite far away one from the other, divergence appears substantial. This assumption is confirmed by Fig. 3.4(b) where the Euclidean distance d(t) along the time t ≤ 2000 can be modeled by an exponential function: -α) with:

x( t) , u( t) y(t), v(t) z( t) , w (t ) (
d(t) = e λ(t
λ ≈ 0.008 α ≈ 1430
From the brief study on the logistic map and the current one on the Lorenz system, the common point is that the chaos appears as soon as the quantity λ is positive. In other words, the sign of the Lyapunov exponent is the signature of chaos. § 3.3 Route to chaos

An aspect that needs investigation is the scenario that is responsible for chaos, in order to find out the situations where it can appear. Actually the chaos theory lists several routes to chaos that are sequences of bifurcations, a bifurcation being any sudden change in the behavior of a system as some parameter is varied [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF].

Period-doubling

For a given value of an external parameter β, the system originally oscillates with a period T. Increasing β makes the oscillation period doubled i.e. T ← 2T. It is actually due to the insertion of a new pattern inside the former periodic one. A larger value of β again doubles the period T ← 2T etc. After a certain time we observe that the system is trapped into a completely non periodic motion. As an illustration is displayed the evolution of the period T of the logistic map according to the parameter a in Fig. 3.5. Even though we can pick a few inter values of the period, the main behavior is such that the period is doubling as a increases. The circles are typical values, also represented by stems in Fig. 3.6(b), such that:

a = 3.421 =⇒ T = 2 a = 3.499 =⇒ T = 4 a = 3.566 =⇒ T = 8 a = 3.604 =⇒ T = 16 a = 3.617 =⇒ T = 32
The period-doubling process may continue until the period becomes infinite, that is the trajectory never repeats itself, it is then chaotic. We display in Fig. 3.6(a) the bifurcation diagram, i.e. the final states x K (K = 100), of the logistic map along the values of a. We observe the first bifurcation that results in oscillations for a = 3 and the next bifurcations that double the period in Fig. 3.6(b).

Quasi-periodicity

We use the Tinkerbell map to illustrate the quasi-periodicity and the associated route to chaos. Such a map is defined by the following state equations:

x k+1 = x 2 k -y 2 k + ax k + by k (3.8) y k+1 = 2x k y k + cx k + dy k (3.9)
and we choose parameters to be:

b = -0.6013, c = 2, d = 0.5.
The value of the parameter a will play the role of a potentiometer that makes quasi-periodicity appear, among others.

Poincaré section Actually this two-dimensional discrete map can be thought of as a reduced three-dimensional continuous map. This assumption comes from the Poincaré section method. The Poincaré section is somehow a stroboscopic portrait of the continuous trajectory, as displayed in Fig. 3.7.

The Poincaré section provides a discrete representation P 1 , P 2 , . . . , P K of the continuous system whose state space dimension is smaller than the original one. Formally, a system like:

Ẋ = f (X, Y, Z) Ẏ = g(X, Y, Z) Ż = h(X, Y, Z)
is mapped to an iterative reduced system:

P k+1 = F(P k ) P • • P 1 • P 2 • P 3 X Y Z Figure 3.7: Poincaré section P
where P k = x k y k such that:

x k+1 = F x (x k , y k ), y k+1 = F y (x k , y k ).
The iterative functions F x , F y are defined by the position of the Poincaré section in the three-dimensional state space.

A periodic three-dimensional continuous system looks like a continuous line in its state space, therefore the projection on the Poincaré section is a finite number of points. For a = 0.5, the Tinkerbell map follows this rule, see Fig. 3.8(a). The power spectrum Fig. 3.8(b) of the sequence {x k } k reveals a single frequency f 0 and its harmonics 2 f 0 , 3 f 0 . Quasi-periodic attractor A quasi-periodic motion presents two or more different frequencies. These frequencies are such that the three-dimensional continuous system is a torus and the Poincaré section is generated by using a plane that cuts through the torus. If the ratio of any two frequencies can be expressed as a rational fraction then the Poincaré section will consist of a finite number of points. In other words, this is a periodic motion. However, if only one ratio of two frequencies cannot be expressed as a ratio of integers, then the ratio is irrational. For this case, the Poincaré map points will eventually fill in a continuous curve, the motion is said to be quasi-periodic because the motion never exactly repeats itself. The points of the Poincaré section drift around the curve forming what is called a drift ring. As an example, we consider the Tinkerbell map for a = 0.1, see Fig. 3.10(a). We observe two leading spectral lines on the spectrum Fig. 3.10(b) whose associated frequencies f 1 , f 2 are such that the ratio is not "strictly" rational. Actually, given that any computer has a finite precision, it is always possible to extract a rational ratio between two arbitrary values, that makes the quasi-periodicity quite hard to detect. For instance, the current values of the two frequencies can be linked such that f 2 ≈ 20 7 f 1 .

-0.4-0.3-0.2-0.1 0 0.1 0.2 0.3 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6
-0.2-0.15-0.1 -5 • 10 -2 0 5 • 10 -2 0.1 0.15 -0.4 A crucial property of the quasi-periodic case is that it is not chaotic. The confusion could be made considering that the chaotic motion also provides a quite continuous curve on the Poincaré section, as we can see with the Tinkerbell map with a = 0.9 in Fig. 3.10(a). The difference comes from the fact that a quasi-periodic behavior is composed of several periodic components then the spectrum contains different spectral lines whereas in a chaotic motion, the spectrum involves a continuum of frequencies, see Fig. 3.10(b).

-0.3 -0.2 -0.1 0 0.1 0.2 0.3

Ruelle-Takens scenario or quasi-periodicity route to chaos

The system begins with a limit cycle. As a control parameter is changed, a second periodicity appears in the behavior of the system. This bifurcation event is a Hopf bifurcation. If the ratio of the period of the second type of motion to the period of the first is not a rational ratio, then the motion becomes quasi-periodic. Under some circumstances, if the control parameter is change further the motion becomes chaotic. One would expect to see a long sequence of different frequencies come in as the control parameter is changed, much like the infinite sequence of period-doubling. However at least in some cases, the system becomes chaotic instead of introducing a third frequency for its motion. In [S.E78] the authors proved that if the trajectory of a system is confined to a three-dimensional torus, i.e. it is made of three frequencies with irrational ratios, then even a small perturbation of the motion will damage the motion and lead to chaos.

-1.4 -1.2-1-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 -2 -1.5 -1 -0.

Intermittency and crises

The intermittency route to chaos is characterized by dynamics with irregular occurring bursts of chaotic behavior interspersed with intervals of apparently periodic behavior. As a control parameter is changed, the chaotic bursts become longer and occur more frequently until, eventually the entire time record is chaotic.

A crisis is a bifurcation event in which a chaotic attractor and its basin of attraction suddenly disappear or suddenly change in size as some control parameter is changed. Alternatively is the parameter is changed in the opposite direction the chaotic attractor can suddenly appear or the size of the attractor can suddenly be reduced. A crisis event involves the interaction between a chaotic attractor and an unstable fixed point or an unstable limit cycle. § 3.4 Lyapunov exponent If the system is behaving chaotically, the divergence of nearby trajectories will manifest itself in the following way: if we select some value x i from the sequence {x k } k and then search the sequence for another value x j that is close to x i then the sequence of differences:

d 0 = |x i -x j | d 1 = |x i+1 -x j+1 | . . . d k = |x i+k -x j+k |
is assumed to increase exponentially from one step k to the next one k + 1, i.e. d k+1 = d k e λ that could be also written:

d k = d 0 e λk (3.10)
where λ is the Lyapunov exponent. According to this assumption we can extract two methods to compute λ.

Computation method 1

The equation ( 3.10) stipulates that the evolution of log d k follows an affine law:

log d k = λk + log d 0 (3.11)
where λ is the slope of the line. A well suited tool to estimate this slope is the method of least squares. According to the theory we can easily show that the estimate of the slope is:

λ = K ∑ k=1 k log d k - 1 K K ∑ k=1 k K ∑ k=1 log d k K ∑ k=1 k 2 - 1 K K ∑ k=1 k 2 (3.12)
where K is the index of the last sample. Obviously, it is crucial to average this value over a large set of initial conditions in order to obtain a quite faithful value that depends much more on the parameters of the system than on its initial states.

Computation method 2

The value d k is the separation between the k-th iterate of the map f from x i and the k-th iterate of f from x j . Given that x i and x j are very close, it is more convenient to consider x i x 0 and x j x 0 + where 1. Thus:

d k = | f k (x 0 + ) -f k (x 0 )| (3.13)
Therefore:

d k = | f k (x 0 + ) -f k (x 0 )| (3.14)
In this way, it appears that:

lim →0 d k = ( f k ) (x 0 ) (3.15)
First of all, the derivatives of any function f k can be decomposed as:

( f k ) (x 0 ) = f (x k ) f (x k-1 ) . . . f (x 0 )
After that, the assumption that the separation evolves exponentially implies that:

( f k ) (x 0 ) = e λk i.e. λ = 1 K K ∑ k=0 log | f (x k )| (3.16)
In other words, the Lyapunov exponent λ is the log-average of the slopes over all the samples x 0 , . . . , x K . On Fig. 3.11 are displayed the estimate of λ from both methods. Even though we can observe few differences, e.g. at a = 3.24 the least squares method does not provide a λ as low as the differential method does, both methods offer estimates of a global same shape. On one hand, the least squares method presents variance for a ∈ [1.24, 3.48] but it is possible to compute λ for any systems, whatever the complexity of its state equations. On the other hand, the differential method is smooth for the whole a interval however it is not always convenient given that the derivatives of the map, or even of the set of maps in the multidimensional case, is not an easy computation task. As an example, it is not well suited for the iterative decoding algorithms as their update rules are sometimes very hard to handle, see chapter 5.

Assumptions

Both computation methods assume two main hypothesis:

• the number of samples K should be large enough to obtain accurate averages,

• the number of samples K should be small enough such that d K can still be approximated by d 0 e λK . For any dissipative system, any attractor whether chaotic or periodic (or quasi-periodic), has a finite size inside the state space. The maximum value of the separation is the attractor size itself. Waiting for too large K iterations would lead to a wrong estimation given that the exponential law of d K is not valid anymore. On Fig. 3.12(a) is displayed the evolution of log d k for the logistic map with a = 3.75. The k values are shared into two sets: for the interval [0, 28] indexed with the blue arrow log d k is linear then the approximation is right, for the interval [29,80] indexed with the red arrow log d k is constant to the size of the chaotic attractor then the approximation does not make sense anymore. The maximum iteration to use is then K = 28. However, the ideal value K is not always obvious, especially for multidimensional systems. The initial conditions coupled with external parameters largely wield influence on the behavior of the system, e.g. the size of the attractor or even its nature (fixed point, periodic, quasi-periodic or chaotic). Therefore finding the value K turns out to be a very hard task. Furthermore, the system can mimic a chaotic behavior for few iterations and then it collapses to a very small sized attractor or even a fixed point. We display such a phenomenon, the transient chaos [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF], in Fig. 3.12(b) for the logistic map with a = 3.56. In such a case, if we still consider that K = 28 then we conclude that the logistic map is chaotic whereas it is absolutely not. Waiting for larger K would lead to the better conclusion that the logistic map has suffered from a temporary almost chaotic behavior. The previous observation about the versatility of K and the current one about the temporary chaos lead us to consider that neither computation method is perfect.

Computation method 3

The issue is to model by a single value the whole behavior of the system along the iterations for any parameter value. The most convenient solution that can be computed for any dynamical system is to use the least squares method for a quite large K. The main consequence is that the output λ will not be the exact Lyapunov exponent because it cannot describe the exponential expansion. Nevertheless, this λ can describe the whole movement of the system. If the system encounters a chaotic attractor, λ is positive even though it is lower than the true Lyapunov exponent. If the behavior is temporary chaotic then the Lyapunov exponent would exhibit a largely positive value whereas the current one would be considerably smaller or even negative, which is very close to the reality. As an example, the computation of Lyapunov exponent using both methods on the logistic map with a = 3.56 in Fig. 3.13 results in:

• Lyapunov exponent: λ = 0.0022,
• alternative Lyapunov exponent: λ = -0.0694 (3.16). Given that the whole evolution of d k is used, not only the transient, therefore the least squares approximation results in least values. That is why around a = 2 the exponent exhibits less negative values than the Lyapunov exponent ones. Furthermore the fact it does not deal with the transient prevents it from being deceiving as we explained previously. Even though its shape is quite similar to the Lyapunov exponent one, it provides a more faithful explanation of the behavior of the logistic map. We do not display any curves about other systems as the Lorenz one. Nevertheless extensive experiments confirmed that the use of our alternative exponent to describe not only the transient but the whole evolution is relevant. The tools of the chaos theory that we have just presented provide information about the behavior of the studied system. The visualization of the associated trajectory in its own state space allows either to confirm the information or to improve it, e.g., the famous Hénon map and the Lorenz map represented in their state space in Fig. 3.16(a) and Fig. 3.15(a).

However, multidimensional systems cannot be so easily represented, given that our human skills prevent us from observing more than three dimensions (or four if we consider time as a dimension). Furthermore, multidimensional dynamical systems are described by so many variables that it appears often very hard to acquire them all. Even in reasonable state space, as the Hénon map one, a few state variables are unavailable. To keep a reliable representation of a system, we need an alternative state space built by the State Space Reconstruction (SSR) [START_REF] Kugiumtzis | State space reconstruction parameters in the analysis of chaotic time series -the role of the time window length[END_REF]. The SSR aims to put up an approximate state space by the only use of a subset of the state variables. Most time, the evolution of a single variable is enough to recover a faithful alternative trajectory, or more generally a one dimensional map of the state variables. Practically, we get back a time series {E(1), . . . , E(K)} of K points that we will map by the SSR to a matrix such that each column is a state variable of the alternative state space:

T =           T 1 T 2 . . . T n . . . T L           =           E(1) E(1 + τ) . . . E(1 + (m -1)τ) E(2) E(2 + τ) . . . E(2 + (m -1)τ) . . . . . . . . . E(n) E(n + τ) . . . E(n + (m -1)τ) . . . . . . . . . E(L) E(L + τ) . . . E(L + (m -1)τ)          
The quantities [T 1 . . . T L ] are L points of the alternative trajectory inside the new state space that represents the system evolution. The value m is the target dimension of the approximate space, called the embedding dimension [H. 06].

The term τ is a delay and L is chosen such that:

L + (m -1)τ ≤ K (3.17)
The SSR is a relevant tool as we are about to see considering the following situations:

• Hénon map: we use only one of the two state variables as the basic sequence to construct the new state space,

• Lorenz map: we use only one of the three state variables as the basic sequence to construct the new state space.

- We display the result in Fig. 3.15(b) and Fig. 3.16(b), with τ = 1, m = 2 for the Hénon map, and τ = 5, m = 3 for the Lorenz map. These parameters are arbitrarily chosen given the embedding dimension value m = 3. We carried out numerous experiments to find the better suited τ, i.e. the delays that provide the most faithful trajectory compared with the original one. The shape of the Hénon map is perfectly restored although the x-axis suffers from a contraction comparing with Fig. 3.15(a). Concerning the Lorenz map, the reconstruction is not perfect. However, the global shape is recovered, which is the goal of the SSR. Indeed, the alternative trajectory provides more qualitative descriptions of our systems than quantitative descriptions. Finally, it appears that the SSR is a good candidate for the representation of a dynamical system. Therefore the SSR will be a useful tool when we will expose the behaviors of the decoding algorithms that are multidimensional systems. We now have the main tools to describe the dynamics of the message-passing algorithms. One of the main task will be to define the state space, that basically conditions the quantifiers to use. Afterwards, it will be useful to redescribe their computation given that the decoders are less simple than the examples we made use of in this chapter. The goal will be to find the parameters that are responsible for non trivial behaviors of the algorithms, and to evaluate the nature of the associated attractors. We will also introduce other quantifiers dedicated to the decoder themselves.

1 0 1 -0.4 -0.2 0 0.2 0.4 X 1 X 2 (a) Hénon map: 2D -0.4 -0.2 0 0.2 0.4 -0.4 -0.2 0 0.2 0.4 X 1 X 2 (b) Hénon map reconstruction

Conclusion

We now have all the necessary tools to tackle the core of the problem, that is the explanation and the analysis of the Belief Propagation algorithm and its generalization, the Generalized Belief Propagation. The statistical physics and the optimization will help together to construct the message-passing equations of both algorithms. By the use of the probability theory on graphical models that we will also expose, we are able to make merge the different research areas in a single problem.

The chapter about the dynamical systems will be utilized in the last chapter, where we will try to confront the decoding algorithms on their stability and their behaviors according to the SNR of the AWGNC. Thanks to the presentation that we made, it will be possible to construct a relevant and consistent analysis of the dynamics of the Belief Propagation and the Generalized Belief Propagation.

57

Introduction

In this second part of the manuscript, we explain in details foundations of the Belief Propagation (BP) and the Generalized Belief Propagation (GBP) algorithms. The previous part helps get accustomed with the different theories needed to tackle this issue. In particular, statistical physics plays a major role in understanding the whole development that leads to the very equations of the decoding algorithms. This is essentially due to the root of this work and many others, that is to say Boltzmann's law. This equation allows us to roll out the whole procedure that concludes on the formulation of the inference problem as a constrained optimization of the free energy.

The BP and the GBP are constructed according to the stationary points of specific formulations of the free energy. As such, they are comparable algorithms, that is why it appears natural to confront them, even though the BP is much more tackled than the GBP in the literature. The relevancy and the originality of our investigation remains in the particular construction of the region-based approximation that we propose to obtain an efficient GBP algorithm. As a matter of fact, we study the topology of a given Tanner graph, that is not usual in the related papers, and we make use of the trapping sets to construct a powerful mean-field generalization, given that they represent one of the most challenging problems in channel coding. In this way, we obtain a relevant GBP algorithm especially for situations where the trapping sets are particularly harmful. We then propose a method to examine and confront the BP and the GBP in these situations.

We organize this part such that the first chapter is a deep investigation on the BP algorithm and its deep connections with the Bethe approximation in statistical physics. We also expose by the use of examples the suboptimality that it is subjected to, due to the loop-like topology of Tanner graphs. In the second chapter, we explain the origin and the functioning of the GBP algorithm based on the region-based approximation. We introduce an original method to factorize the joint distribution that implicitly leads to a specific form of the free energy and then to dedicated message-passing rules. We compare BP and GBP for LDPC codes where trapping sets are responsible for the BER degradation. In the last chapter, we consider the BP and the GBP algorithms as dynamical systems to bring out their properties about their behaviors face to non trivial situations. To this end, we use tools presented in first part, e.g. bifurcation diagram and Lyapunov exponent. We also introduce new tools to connect the dynamical behavior to the decoding performance of 61 both decoders. Thus, we determine sizes and positions of attractors, according to their natures and specific SNR values, provided by the previous dynamical quantifiers.

- In this chapter we show how error correcting codes are closely related to networks handled by BP. We give a few properties of LDPC codes, then we describe the rationale behind the Tanner graph. After that, we introduce a demonstration of the BP equations according to statistical physics. Finally, we present problems linked to BP that will open onto the next chapters. § 4.2 LDPC codes and factor graphs

Error correction capability of LDPC codes

In 1963, Gallager created the Low-Density Parity-Check codes during his PhD [R.G63]. Impractical to implement at this time, they were forgotten until 1996 when McKay and Neal re-discovered them [D. 95]. These codes are aimed at providing fast and relevant decoding algorithms for recovering the input codeword x. Tanner graphs of these codes do not contain many edges, they are sparse, making any message-passing algorithm, such as BP, low-complexity algorithms. To evaluate the error correction capability of a given code, we commonly use d min the minimum distance between any two of its codewords. Coding theory stipulates that the Maximum Likelihood Decoder (MLD) is able to correct a number of errors that equals: Decoders, such as BP, do not behave as MLD, i.e. nothing ensures that errors would be practically corrected. This theoretical aspect makes d min a good indicator to describe LDPC codes performance, whereas the BER is a fair indicator of decoders performance. The work we present in the following does not refer to d min as we focus on decoders constructions, but anyone has to keep in his mind that global performance of a transmission partly rely on the couple code+decoder.

= E d min -1 2 (4.

Tanner graph

Tanner graph is not only a support for message-passing algorithms, it has a deep meaning in the field of graphical models. Previously we pointed out that any channel could be modeled as a random map on the input codeword x. It means that we cannot obtain any deterministic information concerning x from the data y but only a probabilistic information. From the receiver point of view, bits x 1 , . . . , x N are mapped to a probabilistic set made of the random variables X = {X 1 , . . . , X N }. The only feasible process for the receiver when no knowledge of the LDPC code use to encode information is available, is to estimate the a posteriori probability p(X = x, y)1 according to Bayes rule, using the likelihoods, as we mentioned in the previous chapter concerning the naive estimation.

An LDPC code C transforms the probabilistic set X into a probabilistic network (X, C). Any parity-check equation c a is then a correlation function between its argument bits X a . In other words, any conditional probability distribution p(X i |X 1 , . . . , X i-1 , X i+1 , . . . , X N ) is reduced to:

p(X i |X 1 , . . . , X i-1 , X i+1 , . . . , X N ) = p(X i |X N i ), (4.2)
where X N i ⊆ X is the set of variables that are linked to X i by any paritycheck equation (distinct from N i , the set of check nodes linked to X i ), see the example in Fig. 4.1.

X 1 X 2 X 3 X 4 X 5 X 6 X 7 c a c b c c Figure 4.1: Tanner graph -X N 7 = {X 1 , . . . , X 6 }, N 7 = {c a , c b , c c }
Thus, the network is said to fulfill the Markov property. This makes the probabilistic set a Markov Random Field (MRF). We can represent this MRF by a bipartite graph G = (X ∪ C, E) where the vertices are the variables X = {X 1 , . . . , X N } and the check nodes

C = {c 1 , . . . , c M }. Any edge e = (c a , X i ) = (X i , c a ) ∈ E
is the dependence of the variables on the parity-check equations.

Factor graph

The Hammersley-Clifford theorem [P. 90], [J.M71] stipulates that the joint probability distribution of any MRF can be factorized over cliques of the graph G, a clique being defined as a fully connected subgraph. In the current study, any check node c a accompanied by its neighborhood X a constitute a clique. Therefore we obtain the following form:

p(x) p(X = x) = 1 Z M ∏ a=1 f a (x a ) (4.3)
with Z the normalization constant that ensures that p(x) is a genuine probability distribution, i.e. ∑ x p(x) = 1. Functions { f a (x a )} a,x a are symmetric functions called potentials, or kernels, depending on the context:

∀c a ∈ C, ∀X a i , X a j ∈ X a , ∀x a i , x a j , f a (x a 1 , . . . , x a i , . . . , x a j , . . . ) = f a (x a 1 , . . . , x a j , . . . , x a i , . . . ) (4.4)
Such a factorization (4.3) makes any MRF a factor graph. Generally, the factor graph distribution is a product of positive potential functions, but for LDPC codes, potentials are only valued in {0, 1} i.e. the state of any clique c a ∪ X a is strictly valid ( f a (x a ) = 1) or not ( f a (x a ) = 0). In other words, any potential function f a is the XOR of the associated parity-check equations c a :

∀c a ∈ C, ∀x a , f a (x a ) = 1 ⊕ X i ∈X a x i (4.5)
The factor graph of any LDPC codes is called a Tanner graph [F.R01]. Equation (4.3) is the a priori distribution on the variables X 1 , . . . , X N from the receiver point of view. This expression is quite natural in the sense that:

• for any vector that is not a codeword, there is at least one parity-check equation that is not verified, i.e.,

∀x ∈ F N 2 \C, ∃c a , f a (x a ) = 0 =⇒ p(x) = 0
• only the codewords have non-null probabilities to appear,

∀x ∈ C, ∀c a , f a (x a ) = 0 =⇒ p(x) = 0
• the distribution over the codewords is uniform and can be computed knowing that C is built with 2 K codewords (K is the number of information bits encoded by the LDPC code),

∀x ∈ C, p(x) = 1 2 K .
According to Bayes rule, we can finally assert that the a posteriori probability distribution is:

p(x, y) = 1 Z M ∏ a=1 f a (x a ) N ∏ i=1 p(y i |x i ) (4.6)
From the receiver point of view, the estimate of the input sequence x corresponds to: x = arg max x p(x, y) (4.7)

Unfortunately, as it was pointed out about MLD, this distribution is not tractable given that we would need to span through 2 N states to reach x, and N is typically very large (e.g. N = 64800 in LDPC codes used for DVB-S2 standard2 ) making it not practical. To circumvent this problem of computation complexity, we use a trial and tractable probability distribution {b(x)} x , called the joint belief, that we want to be as close as possible to the true distribution {p(x)} x .

In the following, we develop questions about the constraints it is subject to, and how to evaluate its accuracy. § 4.

Pairwise Tanner graph

First of all, we study the case of pairwise Tanner graphs, because they form the simplest examples to understand and handle the factorization of the joint distribution. In such graphs, any edge implicitly carries a check node, because any check node only has two variable nodes as neighbors. We start with a tree, i.e. a loop-free graph, and we continue with a square lattice that is a graph which drawing forms a regular square tiling.

Tree case

We consider the tree represented on Fig. 4.2. Bayes rule and Markov property tell us that the joint distribution is:

X 1 X 2 X 3 X 4
p(x) = p 1|23 (x 1 |x 2 , x 3 )p 234 (x 2 , x 3 , x 4 ) (4.8) p(x) = p 123 (x 1 , x 2 , x 3 ) p 23 (x 2 , x 3 ) p 234 (x 2 , x 3 , x 4 ) (4.9)
As {X 3 , X 1 , X 2 } is a chain, i.e. only one path connects these three nodes, we can easily factorize the term p 123 :

p(x) = p 12 (x 1 , x 2 )p 13 (x 1 , x 3 )p 24 (x 2 , x 4 ) p 1 (x 1 )p 2 (x 2 ) (4.10)
Such a formulation is a product of the couples distributions that are p 12 , p 13 , p 24 divided by the single variable intersections distributions that are p 1 , p 2 :

p(x) = ∏ <X i ,X j > p ij (x i , x j ) ∏ X i p d i -1 i (x i ) (4.11)
where d i is the number of edges incoming on the variable node X i .

Loopy case

We consider the graph represented on Fig. 4.3. The joint distribution is the same as (4.9). However, equation (4.10) cannot be exact on this graph given that X 2 and X 3 are independent conditionally to X 1 and X 4 whereas in the tree case it was only conditioned to X 1 .

X 1 X 2 X 3 X 4
In other words, if we can find more than one path between two nodes in the graph, i.e. if it contains loops, then equation (4.11) holds as an approximation. Most graphs are not trees, therefore factorization (4.11) remains an approximation in most cases. As a famous example in spin glasses theory or image processing, the square lattice of N = 9 variable nodes, represented in Fig. 4.4, can be approximated with the formula (4.11). 4.5(a). We attempt to extend the factorization (4.11) to non-pairwise Tanner graphs. First, if we omit the edges ( f a , X 1 ), ( f b , X 4 ), ( f c , X 3 ), see Fig. 4.5(b), we obtain a loop-free Tanner graph whose joint distribution is exactly the extension of (4.11):

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
X 1 X 2 X 4 X 3 X 5 X 6 X 7 f a f b f c (a) Tanner graph of the Hamming code X 1 X 2 X 4 X 3 X 5 X 6 X 7 f a f b f c (b)
p(x) = p 5|123 (x 5 |x 1 , x 2 , x 3 )p 6|124 (x 6 |x 1 , x 2 , x 4 )p 7|134 (x 7 |x 1 , x 3 , x 4 )(4.12) p(x) = p a (x a )p b (x b )p c (x c ) p 2 1 (x 1 )p 2 (x 2 )p 3 (x 3 )p 4 (x 4 ) (4.13)
Unfortunately, factorization for the true Hamming code cannot be reduced to (4.13):

p(x) = p 5|123 (x 5 |x 1 , x 2 , x 3 )p 6|124 (x 6 |x 1 , x 2 , x 4 )p 7|134 (x 7 |x 1 , x 3 , x 4 )(4.14) = p a (x a )p b (x b )p c (x c )p 1234 (x 1 , x 2 , x 3 , x 4 ) p 123 (x 1 , x 2 , x 3 )p 124 (x 1 , x 2 , x 4 )p 134 (x 1 , x 3 , x 4 ) (4.15)
It is due to the loop-like topology of the graph: several paths link any couple of variable nodes. Therefore the general rule of factorization is only an approximate of the joint distribution:

p(x) ≈ ∏ M a=1 p a (x a ) ∏ N i=1 p d i -1 i (x i ) (4.16)
and so does the belief:

b(x, y) = M ∏ a=1 p a (x a ) ∏ X i ∈X a l i (x i , y i ) N ∏ i=1 p i (x i )l i (x i , y i ) d i -1 (4.17)
It is a quite hard task to evaluate the influence of a given loop on the previous approximation, mainly because loops are very often interleaved altogether.

In addition, we will see that the degradation of the approximation strongly depends on the transmission channel noise, therefore decoders, that we will describe in the following, are the only method to obtain a fair insight.

The factorized form is not the only necessary practical approximation. Indeed, in such a formula, we assume that marginal distributions of potentials {p a (x a )} a,x a and single variable nodes {p i (x i )} i,x i are available. We recall that for any potential f a , for any state x a :

p a (x a ) = ∑ x∪x a p(x) (4.18)
and that for any variable node X i , for any state x i :

p i (x i ) = ∑ x∪x i p(x) (4.19)
These quantities cannot be practically obtained by these marginalizations as they depend on the joint distribution p(x) we are trying to approximate. Therefore, marginal distributions have to be estimated by trial distributions, the beliefs: to any potential f a is associated a belief {b a (x a )} x a and to any single variable X i is associated a belief {b i (x i )} i . § 4.5 Energy

To estimate marginal distributions, we use the form of the joint distribution given by statistical physics. Previously we highlighted that any Tanner graph was a factor graph which probability distribution is:

p(x) = 1 Z M ∏ a=1 f a (x a ) (4.20)
The Hammersley-Clifford theorem stipulates that any factor graph has the distribution of a Boltzmann's field [L. 64]:

p(x) = 1 Z e -E I (x) (4.21)
where E I (x) is an energy function. According to statistical physics, E I (x) provides links or interactions between variable nodes of the network X in the state x, e.g. the magnetic couplings {J ij } i,j between the spins {S i , S j } i,j of a spin glass S. The function E I is linearly decomposed as:

E I (x) = M ∑ a=1 E a (x a ) (4.22)
where each function E a is the interaction energy between variable nodes of the subset X a ⊂ X. This linear decomposition allows us to factorize the joint distribution as:

p(x) = 1 Z M ∏ a=1 e -E a (x a ) (4.23)
that provides, by a simple identification with (4.20), analytical links between potentials { f a (x a )} a,x a and energy functions {E a (x a )} a,x a :

∀c a ∈ C, ∀x a , E a (x a ) = -log f a (x a ) (4.24)
This formulation carries a physical interpretation of LDPC codes. A configuration, or state, x such that at least one parity-check equation c a is not verified, i.e. f a (x a ) = 0, is an infinite energy state. In other words, a forbidden state is made unreachable. In contrast, a state that verifies all parity-check equations is a ground state (state of minimum energy). Thus, to search for the most likely state of the distribution, i.e. x = arg max x p(x), is equivalent to search for the state that minimizes the energy function, x = arg min x E I (x).

Any additive channel noise, e.g. the noise from the Additive White Gaussian Noise Channel (AWGNC), is equivalent to an external magnetic field H in which a spin glass is embedded in. As the AWGNC disturbs independently each bit of a codeword, H interacts independently with each spin. These magnetic interactions are modeled by another energy function:

E H (x) = N ∑ i=1 E i (x i ) (4.25)
which expression is a scalar product between the magnetic field and the spin vectors 3 . The total energy function is then:

E(x, H) = E I (x) + E H (x) (4.26) = M ∑ a=1 E a (x a ) + N ∑ i=1 E i (x i ) (4.27)
Denoting any local interaction of a "spin" X i with the magnetic field H by

E i (x i ) = -log l i (x i , y i )
, where y i = x i + n i is the channel output and n i is a noisy element, we obtain:

E(x, y) = - M ∑ a=1 log f a (x a ) - N ∑ i=1 log l i (x i , y i ) (4.28)
Then, joint distribution becomes:

p(x, y) = 1 Z M ∏ a=1 e -E a (x a ) N ∏ i=1 e -E i (x i ) (4.29) or: p(x, y) = 1 Z M ∏ a=1 f a (x a ) N ∏ i=1 l i (x i , y i ) (4.30)
that equals (4.6). This helps demonstrate that an LDPC code is somehow equivalent to a spin glass, considering that:

• any coupled spins are bits of a common parity-check equation,

• any spin disturbed by an external magnetic field is a bit disturbed by a channel noise. § 4.6 Variational approach Equations (4.17) and (4.28) allow us to construct a variational approach to solve the inference problem. The variational approach principle is to use a trial distribution, the belief {b(x, y)}, that we modify to fit the true joint distribution {p(x, y)} as accurately as possible. Accuracy of any variational approach is evaluated according to a specific metric, as we will see in the following.

We do not take into account the a posteriori probability distribution, but the a priori one. Indeed, as the belief is a probability distribution, it satisfies the Bayes rule:

b(x, y) = b(x)b(y|x) (4.31)
As likelihoods {l i (x i , y i ) = p i (y i |x i )} i,x i are completely independent from factor graphs4 , the probability law of Y|X does not depend on the inference method. As a matter of fact, it helps to discriminate valid states, or codewords, brought out by {b(x)} x , according to a given observation.

This way, we focus on the computation of the a priori probability distribution, see (4.16):

b(x) = M ∏ a=1 b a (x a ) N ∏ i=1 b d i -1 i (x i ) (4.32)
where {b a (x a )} a,x a and {b i (x i )} i,x i are beliefs on potentials and variable nodes, respectively. To evaluate the relevancy of (4.32), a commonly used estimator is the Kullback-Liebler divergence [START_REF] Kullback | Information Theory and Statistics[END_REF]:

KL(b||p) = ∑ x b(x) log b(x) p(x) (4.33)
Reasons to select this estimator are that it offers a physical interpretation of the variational approach as explained in the following, and also that KL(b||p) has a meaning in information theory: it is the relative entropy between distributions b{(x)} x and p{(x)} x . This way, Kullback-Liebler divergence is a measure of the information lost when distribution b{(x)} x is used to approximate distribution p{(x)} x . Main properties of this function are:

• ∀{p(x)} x , ∀{b(x)} x , KL(b||p) ≥ 0, • KL(b||p) = 0 if and only if ∀x, b(x) = p(x),
• KL(•||p) and KL(b||•) are convex functions.

We will not detail the fact that it does not define a metric or a distance because we do not need this piece of information but the interested reader could refer to [START_REF] Kullback | Information Theory and Statistics[END_REF] for more details.

A physical interpretation of the variational approach consists in developing KL(b||p) considering p(x) as a Boltzmann's field:

KL(b||p) = ∑ x b(x) log b(x) -∑ x b(x) log e -E I (x) Z (4.34) = U b -S b + log Z (4.35)
where:

• U b = ∑ x b(x)E I (x)
is the expectation value of the energy function called the variational mean energy provided that the network is described by the law {b(x)} x ,

• S b = -∑ x b(x) log b(x)
is the entropy relative to the distribution {b(x)} x , called the variational entropy.

Properties of the Kullback-Liebler divergence provide that:

U b -S b ≥ -log Z (4.36) U b -S b = -log Z if and only if ∀x, b(x) = p(x) (4.37)
We obtain a dual formulation of the inference problem: the variational approach consists in varying the joint belief to minimize the cost function

F b = U b -S b (4.38)
called the variational free energy. Therefore it is possible to reach analytic expressions of the marginal distributions replacing the trial distribution {b(x)} x by its factorization (4.32). The expression of the variational entropy is:

S b = - M ∑ a=1 ∑ x a b a (x a ) log b a (x a ) - N ∑ i=1 (1 -d i ) ∑ x i b i (x i ) log b i (x i ) (4.39)
The variational mean energy is computed thanks to the linear decomposition of the energy function (4.22):

U b = M ∑ a=1 ∑ x a b a (x a )E a (x a ) (4.40) = - M ∑ a=1 ∑ x a b a (x a ) log f a (x a ) (4.41)
On one hand, as shown in [J.S05], the variational mean energy is not an approximation of the mean energy U p = ∑ x p(x)E(x), the only difference is due to the values of the beliefs5 . On the other hand, even though all beliefs {b a (x a )} a,x a , {b i (x i )} i,x i equal true marginals, the variational entropy is an approximation of the true entropy S p = -∑ x p(x) log p(x) because the factorization of the joint distribution is only an approximation for most factor graphs (the loopy ones). We obtain the equation of the so-called variational free energy of the Bethe approximation:

F b = M ∑ a=1 ∑ x a b a (x a ) log b a (x a ) f a (x a ) + N ∑ i=1 (1 -d i )b i (x i ) log b i (x i ) (4.42) § 4.

Optimization

To search for the most accurate belief is equivalent to optimize the cost function F b i.e. to search for stationary points of the variational free energy. These points are distributions {b(x)} x such that the gradient of F b is zero i.e.:

∂F b ∂b(x) = 0 (4.43)
In the following we present constraints under which a function {b(x)} x is acceptable to minimize F b .

Constraints

Factorization of belief (4.32) is constituted of variational functions {b a (x a )} a,x a and {b i (x i )} i,x i assumed to be distributions, then they have to satisfy local marginalizations:

• (C1) ∀c a ∈ C, ∑ x a b a (x a ) = 1, • (C2) ∀X i ∈ X, ∑ x i b i (x i ) = 1.
In addition, for any check node, distributions of neighboring variable nodes have to equal corresponding marginalizations:

• (C3) ∀c a ∈ C, ∀X i ∈ X a , ∀x i , b i (x i ) = ∑ x a ∪x i b a (x a )
Additional conditions ensure that these functions are non-negative:

• (C4) ∀c a ∈ C, ∀x a , b a (x a ) ≥ 0, • (C5) ∀X i ∈ X, ∀x i , b i (x i ) ≥ 0.
We obtain a constrained problem modeled, and sometimes solved, by the Lagrangian formalism [D.P01]. It is accordingly convenient to remove conditions that are linear combinations of other conditions, which helps lighten writings and computations:

• (C1) and (C3) imply constraint (C2) then (C2) is removable,

• (C3) and (C4) entail condition (C5) thus (C5) is voluntary omitted.

One would notice that the update equation of any belief {b i (x i )} i,x i (1.14), exhibited in previous chapter, included a normalization constant, whereas condition (C1) guarantees such a property. As a matter of fact, all conditions are guaranteed to be satisfied for any optimum while any other points is not assumed to fulfill such conditions. The BP algorithm is an iterative process that could be interrupted at any iteration k either for time's sake or for any other reason. Thus {b

(k) i (x i )} i,x i , {b (k)
a (x a )} a,x a have to define acceptable solutions, even approximative.

Normalization

It is crucial to understand the normalization constraint because failing to do so could lead to a completely wrong optimization. We present here two formulations of the minimization problem that we encountered in the literature. In order to make the explanation quite simple, we reduce the current problem to the following minimization:

f : R n + -→ Y ⊂ R x -→ f (x) such that x ∈ X = {x ∈ R n + , ||x|| = 1}, where ||x|| = ∑ n i=1 x i .
The first formulation of the optimization is:

(P1) :    f (x * ) = min x∈X f (x) x * = arg min x∈X f (x) (4.44)
A modified formulation consists in relaxing the constraint by making it inherent to the function. We accordingly obtain the second formulation of the optimization:

(P2) :

       f x * ||x * || = min x∈R n + f x ||x|| x * = arg min x∈R n + f x ||x|| (4.45)
As we will now explain, (P1) and (P2) are not strictly equivalent. To this end, we consider the example on Fig. 4.6 with n = 2, where function f (x 1 , x 2 ) is assumed to be a paraboloid.

X • • A • • M x 1 f (x 1 , x 2 ) x 2 1 2 1 2 0 1 2 (a) Constrained minimization (P1) X • • A • • • B B • • • C C • • • D D x 1 f (x 1 , x 2 ) x 2 1 2 1 2 0 1 2 (b) Relaxed minimization (P2)

Without any constraint, the minimum of

f is reached in M = [x 1 (M), x 2 (M)]. The normalization constraint is the vertical plan X = {(x 1 , x 2 )|x 1 + x 2 = 1}.
Points of the paraboloid that match with the constraint constitute the parabola. Problem (P1) consists in spanning through X to find the couple

(x * 1 , x * 2 ) such that: ∀(x 1 , x 2 ) ∈ X \(x * 1 , x * 2 ), f (x 1 , x 2 ) ≥ f (x * 1 , x *
2 ) Thanks to the convexity of the parabola, it is guaranteed that we find the point A = [x 1 (A), x 2 (A)], see Fig. 4.6(a).

Considering problem (P2) Fig. 4.6(b), searching for

(x * 1 , x * 2 ) is not restricted to X . For each point (x 1 , x 2 ) ∈ R 2 + , we first compute its normalized version (z 1 = x 1 x 1 +x 2 , z 2 = x 2
x 1 +x 2 ) that lives in X , then we check if f (z 1 , z 2 ) is the minimum. Using this method, the minimum f (A) is always found but the argument of this minimum is not reduced to a unique point. All points on the line that goes from the origin through (x 1 (A), x 2 (A)) are arguments of the minimum because their normalized version equals point A, i.e.:

arg min x∈X f (x) ⊂ arg min x∈R 2 + f x ||x|| (4.46)
Finally, this toy example makes reveal that it is not a mistake not to consider the constraint on the normalization, given that the solution of the optimization is not unique in this case. The current problem of the minimization of the variational free energy F b can then be formulated as searching for:

arg min

{b a (x a )} a,xa ∈R 2 da + {b i (x i )} i,x i ∈R 2 + F b   b a (x a ) ∑ x α b a (x α ) a,x a , b i (x i ) ∑ x j b i (x j ) i,x i   (4.47)
This way, enforcing normalization at every iteration does not prevent the BP from searching for and reaching the optimal solution. Another formulation of the marginals normalizations that provides wrong result is the following one. First of all, we extract {b a (x a )} a,x a , {b i (x i )} i,x i , unnormalized functions that minimize F b . As a post-treatment, we then normalize these functions. In the former example on the paraboloid, this method leads M to be projected on X as is exhibited in Fig. 4.7.

X • • A • A • • M x 1 f (x 1 , x 2 )
x We observe that point A differs from point A, the result is wrong. This error is noticeable when inserting post-normalization quantities {Z a } a , {Z i } i in the equation of the variational free energy. We obtain another function:

Fb = M ∑ a=1 1 Z a ∑ x a b a (x a ) log b a (x a ) f a (x a ) - M ∑ a=1 log Z a + N ∑ i=1 1 -d i Z i b i (x i ) log b i (x i ) - N ∑ i=1 (1 -d i ) log Z i . (4.48)
Signs and values of the new terms strongly affect the variational free energy. Then, there is absolutely no reason that Fb = F b , except if all postnormalization quantities equal 1, that is rarely true without any constraint on F b . Finally, normalizing marginals cannot be done as a post-treatment.

Omitted constraints

We now point out constraints that are not mentioned in the literature, to our knowledge. First of all, the global marginalizations:

∀c a ∈ C, ∀x a , b a (x a ) = ∑ x∪x a b(x) (4.49) ∀X i ∈ X, ∀x i , b i (x i ) = ∑ x∪x i b(x) (4.50)
According to [T. 06a], the reason not to take into account the global marginalization is that the Bethe approximation does not consist in factorizing the joint belief as (4.17) but in decomposing the variational free energy according to (4.42) using pseudo-marginals on check nodes and variable nodes. Indeed, equation of F b can be deduced according to the following principle: splitting variational mean energy and variational entropy over variable nodes and potentials. This principle is original whereas derivation of the joint distribution is somehow a consequence of this principle that was not originally presented in [H.A35]. The factorization helps introduce definitions to link with the inference problem. This is only an approximation and the values of {b(x)} x are never used in the optimization process then it turns out not to be necessary to consider any joint belief in any inference algorithm.

The second constraint that is not mentioned is the global normalization:

∑ x b(x) = 1 (4.51)
As we mentioned previously, the joint belief is somehow a way to create the pseudo-marginals {b a (x a )} a,x a , {b i (x i )} i,x i . But it turns out possible to define these pseudo-marginals without any joint belief, by asserting that we approximate both the mean energy and the variational entropy by linear decompositions. Then any constraint on the joint belief is irrelevant because this quantity does not truly exist. Furthermore conditioning the minimization to a normalized joint belief does not carry any relevant information, it is equivalent to consider that its true form is:

B(x) = 1 Z b(x) (4.52) B(x) = 1 Z ∏ M a=1 b a (x a ) ∏ N i=1 b d i -1 i (x i ) (4.53)
where Z is the normalization such that ∑ x B(x) = 1. The variational free energy becomes:

F B = 1 Z F b -log Z (4.54)
Thus, we observe that Z only implies an affine mapping of F b to F B , that does not affect the optimization problem. Finally we can assert that the use of Z is not necessary. Given that it also includes a non-negligible computation time, we do not take it into account in the optimization. Then we do not consider the optimization of F B but of F b .

Karush-Kuhn-Tucker conditions

Any problem with constraints of inequalities is generally solved using the Lagrangian formalism [D.P01]. To this end, it is necessary to check that all the constraints fulfill the Karush-Kuhn-Tucker conditions that are:

• continuous differentiability: any constraint and its differentiate are continuous,

• regularity: we can find a solution {b * a (x a )} a,x a , {b * i (x i )} i,x i such that none of the differentiates of the constraints with these solutions is zero,

• we can find a set of non-infinite Lagrange multipliers {γ * a } a , {λ ai (x i )} a,i,x i such that:

primal feasibility: the constraints are all verified, -stationarity: the gradient of the Lagrange function is zero.

Given that constraints on the Bethe approximation satisfy KKT conditions, we define Lagrange multipliers {γ a } a , {λ ai (x i )} a,i,x i for all constraints such that the Lagrange function is:

L = F b + ∑ c a ∈C γ a ∑ x a b a (x a ) -1 + ∑ c a ∈C ∑ X i ∈X a ∑ x i λ ai (x i ) ∑ x a b a (x a ) -b i (x i ) (4.55)
KKT conditions are used to define analytic expressions of marginal beliefs. However, one should keep in mind that these KKT conditions are only necessary conditions. The equivalence is not trivial, especially because the variational free energy is rarely convex. Proving equivalence is a very hard task that is not tackled a lot in the literature. In [T. 06a] are introduced a few assumptions that lead to the equivalence but there is no general proof, e.g. it is demonstrated for Tanner graphs that only contain one loop, which leads to a perfect optimization in that case. However, many loops are encountered for most LDPC codes. To our knowledge, no proof ensures the equivalence in these cases, and it is very likely that the equivalence is not verified. Despite this drawback, we previously observed that BP has very decent performance, meaning that this omission is not significantly damageable, even though it is still an open problem.

Solution

We now search for positive marginal functions that minimize the Lagrange functions, therefore we compute the gradient of the Lagrange function ac-cording to the marginal beliefs on check nodes and variable nodes:

∀c a ∈ C, ∀x a , ∂L ∂b a (x a ) = 1 + log b a (x a ) f a (x a ) + γ a + ∑ X i ∈X a λ ai (x i ) ∀X i ∈ X, ∀x i , ∂L ∂b i (x i ) = (1 -d i ) 1 + log b i (x i ) -∑ c a ∈N i λ ai (x i )
Arguments that nullify the gradient are called stationary points of Bethe approximation. We obtain their equations6 :

∀c a ∈ C, ∀x a , b * a (x a ) = f a (x a )e -1-γ a ∏ X i ∈X a e -λ ai (x i ) (4.56) ∀X i ∈ X, ∀x i , b * i (x i ) = e -1 ∏ c a ∈N i e λ ai (x i ) 1-d i (4.57)
We observe that these functions are positive then it is not useful to consider constraints (C4) and (C5), that is why they do not appear in the equation of the Lagrange function (4.55). As for any variable node Using constraints, it is possible to obtain a system of equations that does not depend on beliefs anymore but only on Lagrange multipliers, e.g. constraint (C3) is equivalent to:

X i , b i (x i , y i ) = b i (x i )l i (x i , y i ) and for any check node c a , b a (x a , y a ) = b a (x a )l a (x a , y a ) with l a (x a , y a ) = ∏ X i ∈X a l i (x i , y i ) we obtain: ∀c a ∈ C, ∀x a , b * a (x a , y a ) = l a (x a , y a ) f a (x a )e -1-γ a ∏ X i ∈X a e -λ ai (x i ) (4.58) ∀X i ∈ X, ∀x i , b * i (x i , y i ) = l i (x i , y i )e -1 ∏ c a ∈N i e λ ai (x i ) 1-d i
e λ ai (x i ) 1-d i = e -γ a ∑ x a ∪x i f a (x a ) ∏ X j ∈X a \X i e -λ aj (x j ) e λ ai (x i ) ∏ X b X i , f b = f a e λ bi (x i ) 1-d i (4.60)
However, the system is not linear and no method in the literature indicates how to handle equations to reach a solution. Partly tackled in [T. 06a] and [A.L03], computation of multipliers is very hard to handle but it emerged an alternative algorithm, the CCCP 1.6 to solve it. In spite of encouraging results [J.

-13a], CCCP does not success for any LDPC code, especially due their looplike topologies.

Another approach to this dual problem is offered by the BP algorithm. First of all, we assume that:

e -λ ai (x i ) = ∏ X b X i ,c b =c a e λ bi (x i ) 1-d i (4.61)
Afterwards, we can reduce (4.60):

e λ ai (x i ) 1-d i = e -γ a ∑ x a ∪x i f a (x a ) ∏ X j ∈X a \X i e -λ aj (x j ) (4.62)
These two equations are messages equations of the BP algorithm provided that messages and Lagrange multipliers are linked such that:

m ai (x i ) = e λ ai (x i ) 1-d i (4.63) n ia (x i ) = e -λ ai (x i ) (4.64)
Using an iterative process to solve these Lagrange multipliers is performing the message-passing of BP (1.12) and (1.13). We then obtain beliefs equations of BP:

b a (x a , y a ) = 1 Z a l a (x a , y a ) f a (x a ) ∏ X i ∈X a n ia (x i ) (4.65) b i (x i , y i ) = 1 Z i l i (x i , y i ) ∏ c a ∈N i m ai (x i ) (4.66)
where Z i , Z a are normalization constants. Identifying messages to Lagrange multipliers provides the crucial result that stationary points of the variational free energy of Bethe are fixed points of BP. We summarize three equivalences that we have encountered until now:

• searching for the most likely state of Tanner graphs is searching for the ground state of energy functions,

• searching for the most accurate belief is minimizing the variational free energy,

• searching for BP fixed points is searching for stationary points of the variational free energy of Bethe. § 4.8 BP assumption

Many experiments have shown in the literature that BP is quite efficient to decode LDPC codes, and more generally to infer any factor graph. However, its equations are based on a particular assumption that prevents it from being optimal. To understand it well, we use two very simple examples in which we extract the belief on a given variable node.

Tree-like Tanner graph

A Tanner graph has a tree-like topology if and only if any two nodes are linked by a unique path, e.g. Tanner graph in Fig. 4.8. In this case, messages and beliefs of BP reach canonical forms, that is to say expressions that only depend on parity-check equations and likelihoods. In previous figure, after few iterations, the marginal belief of variable node X 1 becomes:

X 1 X 2 X 3 X 4 X 5 X 6 X 7
b (k) 1 (x 1 , y 1 ) = l 1 (x 1 , y 1 ) Z ∑
x 2 ,x 3 ,x 4 ,x 5 ,x 6 ,x 7 l 2 (x 2 , y 2 )l 3 (x 3 , y 3 )l 4 (x 4 , y 4 ) ×l 5 (x 5 , y 5 )l 6 (x 6 , y 6 )l 7 (x 7 , y 7 )

× f a (x 1 , x 2 , x 3 ) f b (x 2 , x 4 , x 5 ) f c (x 3 , x 6 , x 7 ) (4.67)
We observe that it does not depend on messages or initialization quantities due to the iterative process but only on intrinsic parameters of the problem. As a matter of fact, anyone can verify that for any other variable node X i :

∀x i , b i (x i , y i ) = ∑ x∪x i 1 Z M ∏ a=1 f a (x a ) N ∏ i=1 l i (x i , y i ) (4.68)
which is the exact marginalization of the true joint distribution (4.6). It finally confirms that BP reaches optimal performance for tree-like LDPC codes.

Non tree-like Tanner graph

We use the previous Tanner graph considering variable nodes X 5 and X 6 merged, see Fig. 4.9. Follow BP update equations, we obtain at any iteration k:

X 1 X 2 X 3 X 4 X 5 X 6
(k) 2a (x 2 , y 2 ) = l 2 (x 2 , y 2 ) Z ∑ x 1 ,x 2 ,x 3 ,x 4 ,x 5 ,x 6 f b (x 2 , x 4 , x 5 ) f c (x 3 , x 5 , x 6 )
×l 3 (x 3 , y 3 )l 4 (x 4 , y 4 )l 5 (x 5 , y 5 )l 6 (x 6 , y 6 )

×m (k-3) 1a (x 1 , y 1 )m (k-3) 2a (x 2 , y 2 ) (4.69) Message m (k)
2a depends on its value three iterations before m

(k-3) 2a , due to the loop (c a ) -(X 2 ) -(c b ) -(X 5 ) -(c c ) -(X 3 ) -(c a ).

Loops

We show in Fig. 4.10 local Tanner graphs used to compute BP equations (1.12) and (1.13). Crosses indicate that BP assumes that there is a unique path between two nodes of the Tanner graph, we call this the BP assumption. Therefore, this assumption makes us temporarily ignore few edges when computing messages. Anyone can show that this prevents us from obtaining canonical equations of messages and then beliefs, as we previously saw (4.69). Thus, tree-like graphs result in optimal performances, contrary to loopy graphs.

X i X j X k X l c a (a) Message from c a to X i c a c b c c c d X i (b) Message from X i to c a
In addition, iterative process involves that errors due to loops are spread out in the graph depending on loops sizes: the shorter a loop is, the more harmful it is for the approximation as its effect appears in a number of iterations equal to the number of variable nodes it contains.

Designing LDPC codes of a high minimum distance needs sophisticated links between variable nodes to make them robust against channel noise. Regularly, robustness is relevant when any bit takes part in more than two paritycheck equations. The other side of the coin is that it creates loops in Tanner graphs. Challenge of coding theory is then to design codes with a high minimum distance and with loops as large as possible. Optimal solution then appear unreachable by the BP algorithm.

Physical interpretation

Considering the equivalence with statistical physics, such a result is not surprising. Indeed, Bethe approximation is a MF approach, see in first chapters 2.4, which is obvious when confronting both fundamental equations:

(MF) : b(x) = N ∏ i=1 b i (x i ), (Bethe) : b(x) = ∏ M a=1 b a (x a ) ∏ N i=1 b d i -1 i (x i ) (4.70)
It means that instead of ignoring couplings when computing marginal distributions on single variable nodes, as originally done in MF equations, Bethe approximation is aimed at neglecting longer range correlations between spins, that are loops. In other words, Bethe approximation extends MF approach. First of all, we will see in next chapter that this note considerably help improving again MF approaches. In addition, this also helps not to consider the BP assumption as a theoretical mistake but as a limit that goes further than the MF one.

Oscillation

An effect due to the loop-like topology of Tanner graphs is the oscillation of messages, and then of the beliefs. A simple example helps understand well this phenomenon. We consider the Tanner graph in Fig. 4.11. The only acceptable states are x ∈ {0, 1}. We now assume that likelihoods are:

X 1 X 2 X 3 X 4
l 1 (x 1 , y 1 ) = 0.75 0.25 , l 2 (x 2 , y 2 ) = 0.25 0.75 , l 3 (x 3 , y 3 ) = 0.25 0.75 , l 4 (x 4 , y 4 ) = 0.75 0.25 .

In this case, we can show, according to BP equations, that beliefs oscillate such that for any iteration k:

b (k) 1 (x 1 , y 1 ) = 0.25 0.75 , b (k+1) 1 (x 1 , y 1 ) = 0.75 0.25 , b (k) 2 (x 2 , y 2 ) = 0.75 0.25 , b (k+1) 2 (x 2 , y 2 ) = 0.25 0.75 , b (k) 3 (x 3 , y 3 ) = 0.75 0.25 , b (k+1) 3 (x 3 , y 3 ) = 0.25 0.75 , b (k) 4 (x 4 , y 4 ) = 0.25 0.75 , b (k+1) 4
(x 4 , y 4 ) = 0.75 0.25 .

The estimate of the input codeword then also oscillates:

x (k) = [1001], x (k+1) = [0110]
Both states 1001 and 0110 are unacceptable because they do not satisfy all parity-check equations. More precisely, for a given variable node, neighboring check nodes carry conflicting information about its most likely state. It is then impossible to simultaneously satisfy all parity-check equations. Such a behavior prevents BP from reaching any steady state, which clearly damages its performance. As we saw in chapter 2.3, oscillations also appear in statistical physics for spin glasses and are referred to as frustration.

Trapping-sets

The BP assumption can seriously damage loopy LDPC codes decoding. Worse than loops are combinations of loops, also known as trapping sets [T.J03]. A trapping set is a nonempty set of variable nodes in a Tanner graph that are not eventually corrected by the BP decoder7 . An (a, b) trapping set, shortly written TS(a, b), is a set of a variable nodes such that the induced subgraph has b odd-degree, or unverified, check nodes. A trapping sets ontology has been created for LDPC codes which variable nodes have three checks neighbors, and which shortest loops contain eight nodes [B. 09]. In this work, we focus on TS (5,3) represented in Fig. 4.12.

Trapping sets are especially investigated on discrete channels, particularly BSC, because in these cases, they match low-weight errors events, i.e. sequences of erroneous bits that BP cannot eventually decode. They are responsible for the error floor [T.J03], an abrupt degradation of the error-rate slope that occurs for low crossover probabilities of BSC, see Fig. 4.13. Researchers deeply investigate trapping sets and their relations to failure of iterative decoders. Two points of view remain to circumvent the error-floor:

X 1 X 2 X 3 X 4 X 5
• either we construct LDPC codes that contain neither small loops nor combinations of them [D.V12],

• or we design decoders able to manage trapping sets.

First point of view has already been mentioned when we introduced loops. In the second one, Finite Alphabet Iterative Decoders (FAID), strongly studied in [S.K11] and [D. 12], can significantly improve error-rate for very low values of the crossover probability. Two drawbacks of these algorithms in our study are that they are not defined on LDPC codes with arbitrary node degrees and they do not deal with any continuous channels, as the AWGNC. In our study, we choose not to discriminate between any of these channels, therefore we focus on a generalization of BP, described in next chapter, studied on BSC and AWGNC. This decoder would help improve the BP decoding performance, annihilating loops and often trapping sets effects.

BP and Bethe approximation are demonstrated to be equivalent, therefore they are two different points of view on a common problem. This result help interpret any Tanner graph as a physical network, e.g. spin glasses, which can be studied according to a variational approach. Indeed, minimizing the variational free energy of Bethe using trial distributions is searching for marginal beliefs as close as possible to exact marginals. We saw that BP -to some extent Bethe approximation too -was based on a crucial assumption telling that loops must be ignored to get messages update equations. According to statistical physics, this is somehow an extension of the MF approach that even ignores couplings between spins (only when factorizing joint distribution). Literature on coding theory indicates that most LDPC codes are full of loops and trapping sets which then prevent BP from being an optimal decoder. Thus, searching for a more general MF approach turns out to be necessary to overcome this very important problem.

-Chapter 5 -

Region-based approximation: beyond the Belief Propagation § 5.1 Introduction

In this chapter, we present an extension of the BP algorithm, the Generalized Belief Propagation (GBP) algorithm. Stemming from the cluster variation method that gave birth to the Kikuchi approximation, also known as regionbased approximation, in statistical physics, this method provides a more general approach to factorize the joint distribution. This is actually the general form that the Bethe approximation and the Mean-Field are particular cases of. It helps deal with the loops problem inside Tanner graphs of LDPC code extending nodes to the concept of regions, where regions are defined as subgraphs, or gathering of nodes. This allows us to build a mapping of any Tanner graph to a region-graph where harmful loop-like structures are absorbed in regions.

We begin this chapter introducing the region-based formalism for any factor graph. The region-based approximation is then tackled from the point of view of the inference and statistical physics, that makes emerge a new variational free energy called the region-based free energy. This stands for a new cost function which stationary points provide equations of a new messagepassing, the GBP algorithm. We also describe general damping rules that help the decoder converge, contrary to its original version. Afterward, we present our original rules of the region-graph construction based on trapping sets. In addition, we exhibit results that reflect significant performance improvements compared with the BP, according to an original method that emphasizes on error-events based on trapping sets.

A part of the following results are published in [J.-12c] and [J.-12a]. § 5.2 The region-based formalism Bethe approximation, presented in previous chapter, mostly relies on an approximate factorization of the joint distribution:

p(x) ≈ ∏ M a=1 p a (x a ) ∏ N i=1 p d i -1 i (x i ) (5.1)
which is exact if and only if the corresponding factor graph is a tree. We recall that this formulation stems from Bayes rule and Markov condition 1 . As a matter of fact, Markov condition does not restrict the factorization of the joint distribution to be only (5.1). More precisely, we are allowed to extend Markov condition asserting that the distribution of any set of variable nodes only depends on the neighborhoods of all these nodes. This way, we find other factorizations that can be more accurate or even exact. To illustrate such tellings, we consider the factor graph in Fig. 5.1. Bayes rule allows us to formulate the joint distribution as: (5.3) Both factorizations, either exact or approximative, are also distributions of other graphical representations, see Fig. 5.2 and Fig. 5.3. For these factorizations, marginal distributions are computed over subgraphs of the factor graph. These subgraphs are called regions, a region r being made of:

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
p(x) = p 124578 (x 1 , x 2 ,
• a set of variable nodes X r ⊆ X,

• a set of check nodes C r ⊆ C 3 that can be an empty set, 1 that stipulates that the distribution of any variable node only depends on its neighborhood 2 any two nodes are only linked by a unique path 3 provided that any check node c a of C r is accompanied by its neighborhood X a ⊆ X r because a region has to be a subgraph (associated edges are implicitly included in r)

The regions set used to write down the factorization is denoted by R. Joint distribution of a factor graph is then approximated by: where c r denotes what is called the counting number of r, that is the weight of its distribution inside the joint distribution, e.g. in (5.3) region {X 1 , X 2 , X 4 , X 5 } has a unit weight, and region {X 2 , X 5 } has a weight of -1. This formulation is a region-based formulation of the distribution.

p(x) ≈ ∏ r∈R p c r r (x r ) (5.4) X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 1 , X 2 , X 4 , X 5 , X 7 , X 8 X 2 , X 5 , X 8 X 2 , X 3 , X 5 , X 6 , X 8 , X 9

Links with MF approach

At this point appears a deep link with the mean-field approach we exposed in chapter 2.4. We recall that MF consists in approximating the joint distribution by the following factorization:

p MF (x) = N ∏ i=1 p i (x i ) (5.5)
This corresponds to (5.4) with any region being a single variable node with a unit weight. The rationale behind the MF is that any single variable node considers all the other variable nodes as a unique global field or node. The Bethe approximation, linked to the factorization (5.1), is a first extension of the MF. In this case, cliques of variable nodes consider the other cliques, and to some extent the other variable nodes, as a unique global node. In this way, the joint belief could approach the exact form of the joint distribution, and therefore it could be more accurate than the MF. As an example, the Bethe approximation is exact if the factor graph is a tree whereas the MF cannot be.

In the current case, we are given more freedom to write the factorization. It allows us to deal with bigger subsets of variable nodes, and to some extent to deal with harmful topological structures e.g. loops and trapping sets, as we will see in the following.

Counting numbers

In Bethe approximation, factorization of the joint distribution is:

b(x) = M ∏ a=1 b a (x a ) N ∏ i=1 b 1-d i i (x i ) (5.6)
Comparing with the region-based factorization, we observe that any potential has a unit weight, and any variable node X i has a weight of 1d i . These weights directly stem from Bayes rule and Markov condition, such that any single variable node -and to some extent any potential -has to contribute only once to the factorization. According to the Bethe approximation, considering an arbitrary variable node X i , we observe that:

• X i appears in d i marginal functions on potentials,

• X i appears once in marginal {b i (x i )} x i of multiplicity 1d i .

As d i + (1d i ) = 1, we are ensured that X i only "participates" once to the factorization. It is somehow as chemical equations where each element appears the same number of times on both equation sides, that we can model:

b(x 1 , x 2 , . . . , x N ) -→ M ∏ a=1 b a (x a ) N ∏ i=1 b 1-d i i (x i ) (5.7)
In [J.S05], {c r } r are called counting numbers, and authors introduced the concept of validity of a region-based approximation when it satisfies this principle.

In the following, all region-based approximations will fulfill this rule.

Computation rule

For any couple of regions (r, s) in a given regions set R, s ⊆ r means that X s ⊆ X r and C s ⊆ C r . We previously mentioned that regions distributions have to be weighted not to over count inner variable nodes and potentials in the factorization of the joint distribution (5.4). This is modeled according to following equations:

∀X i ∈ X, ∑ r∈R X i ∈X r c r = 1, (5.8) ∀c a ∈ C, ∑ r∈R c a ∈C r c r = 1 (5.9)
These equations can be rephrased for regions, such that any region only contributes once to the region-based factorization:

∀r ∈ R, ∑ s⊇r c s = 1 (5.10) which can be rephrased as:

∀r ∈ R, c r = 1 -∑ s⊃r c s = 1 (5.11)
This helps compute counting numbers. First of all, any region r that is not included in any other region of the same regions set R has a unit counting number according to the previous equation. Afterward, we can compute any region that is included in other regions of unit counting numbers, and so on. We iteratively compute counting numbers for all regions according to (5.11). § 5.

Region-based approximation

Contrary to Bethe approximation that provides a unique joint belief, regionbased approximation is not so restrictive. Constructed from cluster variation method of Kikuchi [R. 51], it offers us numerous choices of the regions set R to factorize b(x). As indicated in previous chapter, variable nodes distributions are not available, that is still true for regions. The joint distribution is accordingly approximated by a joint belief computed with a family of functions {b r (x r )} r,x r that are marginal beliefs on regions:

b(x) = ∏ r∈R b c r r (x r ) (5.12)
Region-based approximation allows us to extend MF approach to longer ranges than Bethe approximation. Any family of functions {b r (x r )} r,x r is then expected to provide better performance than a family {b a (x a )} a,x a , {b i (x i ) i,x i } i,x i .

Example

We present here an evaluation of the accuracy of Mean-Field (MF), Bethe Approximation (BA) and Region-Based Approximation (RBA). Following the experimental protocol of [J.S05], we consider a square lattice of size N = 9, see Fig. 5.1, with soft factors. More precisely, as is commonly regarded in physics of spin glasses, we consider that for each node c a :

f a (x i , x j ) = e J ij e -J ij e -J ij e J ij
where {J ij } <i,j> are centered Gaussian variables of variance σ used as the control parameter of the experiment. We display in Fig. 5.4 the average error:

= ∑ x |b(x) -p(x)| 2 N (5.13)
for every approximation. First of all, we see that the MF error is the largest one. It is consistent with the fact that MF does not consider any graphical link between variable nodes, i.e. it considers them independent although they are correlated4 . We observe then that RBA is globally closer to the true distribution than other approximations. We also observe that the BA reaches a maximum near σ = 1.4 that is not an expected phenomenon, but we do not investigate it here because it does not deal with the current problem.

Region graph

Any factorization stemming from RBA (5.12) is graphically represented by a Bayesian network, e.g. Fig. 5.2 and Fig. 5.3, which is called region-graph, any region r being one of its node. To present construction rules of such a graph, specific subsets have to be defined. Given a regions set R and r ∈ R :

• F r = {q ∈ R|q ⊆ r} is family of r,

• D r = F r \r are descendants of r,

• E r = {q ∈ D r | ∃s ∈ D r , q ⊂ s ⊂ r} are children of r,

• P r = {q ∈ R\F r |r ∈ E q } are parents of r,

• A r = {q ∈ R\F r |r ∈ D q } are ancestors of r, Any region-graph R = {R l } 1≤l≤L is constructed level by level, or generation by generation, a level R l being a set of regions that have the same number of variable nodes and check nodes. First level contains regions called clusters, arbitrary selected according to one rule: all variables nodes and check nodes have to be included in these clusters. We say that clusters, or the clustering, entirely cover the factor graph, see Fig. 5.5.

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
(a) Clustering by adjacent squares

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
(b) Clustering by squares A level R l distinct from the first level is constructed according to the Bayes rule: we search for all intersections between regions of R l-1 . In other words, we exhibit common children between regions of R l-1 :

r = X r ∪ C r ∈ R l ⇐⇒ ∃(p, q) ∈ R 2 l-1 s.t. X r = X p ∩ X q C r = C p ∩ C q
Edges of the region-graph represent links between any region and its parents. The degree of freedom in a region-graph construction is then the Tanner graph clustering. Considering example in Fig. 5.5, we can exhibit at least two region-graphs Fig. 5.2 and Fig. 5.3 constructed from two different clusterings.

The way clusters are picked out widely influences RBA performance. In previous chapter we highlighted that the BP algorithm was suboptimal because of loop-like topology of most LDPC codes. RBA is aimed at absorbing harmful topological structures, e.g. loops and trapping sets, in clusters in such a way that they are quenched: no messages and then no information run on their Tanner graph edges. Thus, effects of topological damageable structures are hoped to be reduced or even annihilated, making region-based approximation a more efficient method than Bethe approximation.

X 1 , X 2 , X 4 , X 5 , X 7 , X 8 X 2 , X 3 , X 5 , X 6 , X 8 , X 9 X 2 , X 5 , X 8 (a) Region-graph by clustering with adjacent squares As we already mentioned, the selection of clusters is not unique, the only conditions that must be fulfilled are:

X 1 , X 2 , X 4 , X 5 X 2 , X 3 , X 5 , X 6 X 4 , X 5 , X 7 , X 8 X 5 , X 6 , X 8 , X 9 X 2 , X 5 X 4 , X 5 X 5 , X 6 X 5 , X 8
• the whole Tanner graph has to be covered i.e. any variable node and any check node has to be included in at least one region,

• any check node cannot be separated from the variable nodes of its neighborhood.

A Tanner graph can then be clustered by a large variety of region-graphs. The difficult task is to single out the region-graph with the most accurate marginal beliefs. It is by no means a trivial task, as already mentioned in [J.S05] "how to optimally choose regions [...] remains at this point more an art than a science".

In our study, we investigate this "art" for particular cases given that, to our knowledge, there is no general and efficient method for any graph. § 5.4 Generalized Belief Propagation

In this section, we make emerge the message-passing equations for any regiongraph following the same method as Bethe approximation.

Variational free energy

As mentioned concerning Bethe approximation, the variational free energy is a cost function which minimization provides beliefs and messages equations. We recall that this function, denoted by F K in the case of region-based approximation, is decomposed as the sum of the variational mean energy U K and the variational entropy S k : 5.14) with:

F K = U K -S K (
U K = ∑ x b(x)E(x) (5.15) S K = -∑ x b(x) log b(x) (5.16)
We recall that the joint belief has the following expression:

b(x) = ∏ r∈R b c r r (x r ) (5.17)
Therefore we obtain for the variational entropy:

S K = -∑ r∈R c r ∑ x r b r (x r ) log b r (x r ) (5.18)
provided that the family is defined such that:

∀r ∈ R, ∀x r , b r (x r ) = ∑ x∪x r b(x) (5.19)
Indeed, the equation ( 5.18) does not equal (5.16) if this condition is not respected. However, this condition is hard to handle, as mentioned in 4.7.1.2.

To be consistent with the construction of Bethe approximation, we consider that the entropy (5. 16) is approximated by a linear decomposition over regions entropies (5.18), that is the region-based entropy. The region-based mean energy is computed using the expression of the energy function from the previous chapter (4.26). It is possible expand it over regions: .20) with:

E(x) = ∑ r∈R c r E r (x r ) ( 5
E r (x r ) = ∑ c a ∈C r E a (x a ) (5.21)
We assign: E r (x r ) =log f r (x r ) (5.22) that amounts to define a potential function of the region:

f r (x r ) = ∏ c a ∈C r f a (x a ) (5.23)
Then the region-based mean energy is:

U K = -∑ r∈R c r ∑ x r b r (x r ) log f r (x r ) (5.24)
We recall that this linear decomposition is exact as it stems from the energy function expression and not from the belief factorization, contrary to the region-based entropy. Finally the region-based free energy becomes:

F K = ∑ r∈R c r ∑ x r b r (x r ) log b r (x r ) f r (x r ) (5.25)

Optimization

As mentioned in previous chapter, we know that minimizing F K under the following constraints provides marginal beliefs equations of regions:

• (C1) ∀r ∈ R, ∑ x r b r (x r ) = 1, • (C2) ∀r ∈ R, ∀s ∈ E r , ∀x s , b s (x s ) = ∑ x r ∪x s b r (x r ).
We include neither inequality constraints on the beliefs, nor global normalizations and marginalizations, see previous chapter for details. We obtain the Lagrange function:

L K = F k + ∑ r∈R γ r ∑ x r b r (x r ) -1 (5.26) + ∑ r∈R ∑ s∈E r ∑ x s λ rs (x s ) ∑ x r ∪x s b r (x r ) -b s (x s ) (5.27)
where {γ r } r and {λ rs (x s )} r,s,x s are Lagrange multipliers associated with con- straints. Anyone can check that RBA satisfies KKT assumptions, as Bethe approximation. Lagrange multipliers are found such that critical points, i.e. beliefs that nullify the gradient of L K , are: ( 5.29) where:

b * r (x r ) = f r (x r )e -1 e -1 cr γ r ∏
l r (x r , y r ) = ∏ X i ∈X r l i (x i , y i ) (5.30)

Rotation of the Lagrange multipliers

Bethe approximation stipulated that the belief {b a (x a )} x a of a check node c a ∈ C was the product of all incoming messages in any of its neighboring variable nodes (4.58). We also saw that these messages were directly related to Lagrange multipliers. Thus, generalizing it for a region would be that the belief {b r (x r )} x r of a region r ∈ R is the product of all incoming messages in any region q of its family F r , messages being equivalently related to the Lagrange multipliers {λ qr (x r } q,x r . However, this assumption is not consistent with equation (5.28). Here we demonstrate that an equivalent formulation of the constraints may make emerge the right generalization of the Bethe approximation. First of all, we define Lagrange multipliers {µ pr (x r )} p,r,x r associated with the alternative constraint:

(C2') ∀r ∈ R, ∀p ∈ P r , ∀x r , c r b r (x r ) + ∑ q∈A r \p∪A p c q ∑

x q ∪x r b q (x q ) = 0

Then we prove that constraint (C2) associated with multipliers {λ pr (x r )} p,r,x r is equivalent to constraint (C2').

From λ to µ

We consider an arbitrary region r ∈ R. We recall that its counting number c r is given by:

c r = 1 -∑ q∈A r c q (5.31)
We rephrase it such that:

c r + ∑ q∈A r c q = 1 (5.32)
We consider a parent region p of r, then: We recall that constraint (C2) associated with {λ qr (x r )} q,x r is assumed to be true. Therefore we obtain:

c p + ∑ q∈A p c q = 1 (5.
c r b r (x r ) + ∑ q∈A r \p∪A p c q ∑ x q ∪x r b q (x q ) = 0 (5.36)
that is constraint (C2').

From µ to λ

This part of the equivalence is induction. We consider a region r ∈ R 2 , and two of its parents p 1 , p 2 . We assume that constraint (C2') is true: 5.37) and:

c r b r (x r ) + ∑ q∈P r \p 1 ∑ x q ∪x r b q (x q ) = 0 (
c r b r (x r ) + ∑ q∈P r \p 2

∑

x q ∪x r b q (x q ) = 0 (5.38) Subtracting ( 5.38) to (5.37) gives:

∑ x p 1 ∪x r b p 1 (x p 1 ) = ∑ x p 2 ∪x r b p 2 (x p 2 ) (5.39)
In other words, marginalization does not depend on the parent of r. The fact that r ∈ R 2 involves that:

• c r = 1 -N p , with N p the number of parents of r,

• ∀p ∈ P r , c p = 1.

Considering a parent p of r, equation ( 5.37) becomes:

(1

-N p )b r (x r ) + (N p -1) ∑ x p ∪x r b p (x p ) = 0 (5.40) that implies: b r (x r ) = ∑ x p ∪x r b p (x p ) (5.41)
provided that N p > 15 . Finally we obtain constraint (C2) associated with λ pr (x r ).

First part of the induction is demonstrated, now we have to work on a region r of an arbitrary generation R l . We can prove by another induction [J.S05],[R.J02] that:

(c r -1)b r (x r ) + ∑ q∈A r c q ∑

x q \x r b q (x q ) = 0 (5.42) Equations ( 5.42) and (C2') imply that for any ancestor t of r:

c r b r (x r ) + ∑ q∈A r c q ∑ x q ∪x r b q (x q ) -∑ x t ∪x r b t (x t ) = 0 (5.43)
For any couple of ancestors t 1 , t 2 of r, subtracting (5.43) for t 1 to (5.43) for t 2 results in::

∑ x t 1 ∪x r b t 1 (x t 1 ) = ∑ x t 2 ∪x r b t 2 (x t 2 ) (5.44)
In words, it means that marginalization does not depend on the selected ancestor, for any region in R. It implies in (C2') that for any ancestor (to some extent for any parent) p: b r (x r ) = ∑

x p ∪x r b p (x p ) (5.45) which is constraint associated with the Lagrange multiplier λ pr (x r ).

As we proved that constraints (C2) and (C2') are equivalent, the Lagrangian function of F K is rephrased:

L K = F k + ∑ r∈R γ r ∑ x r b r (x r ) -1 (5.46) + ∑ r∈R ∑ p∈P r ∑ x r µ pr (x r ) c r b r (x r ) + ∑ q∈A r \p∪A p c q ∑ x q ∪x r b q (x q )
whose stationary points are:

b * r (x r ) = f r (x r )e -1 e -γr cr ∏ p∈P r e -µ pr (x r ) ∏ q∈D r ∏ s∈P q \F r e -µ sq (x q ) (5.47)
which correspond to a generalization of Bethe approximation we searched for.

Considering observations Y r = {Y i |X i ∈ X r } and Bayes rule, we obtain the a posteriori distribution:

b * r (x r , y r ) = l r (x r , y r ) f r (x r )e -1 e -γr cr ∏ p∈P r e -µ pr (x r ) ∏ q∈D r ∏ s∈P q \F r
e -µ sq (x q ) (5.48)

Message-passing

Belief equation ( 5.48) on a region r is, in words, the product of all information concerning it and its family. This information comes from:

• the code itself represented by f r (x r ),

• the channel represented by l r (x r , y r ),

• the rest of the region-graph represented by all the Lagrange multipliers {µ pr (x r )} p∈P r ,x r , {µ sq (x q )} q∈D r ,s∈P q \F r ,x q

In other words we can summarize it in the following formula:

b r (x r , y r ) = l r (x r , y r ) f r (x r )m R→r (x r ) (5.49) where m R→r (x r ) is information coming from the region-graph that can be re-written as:

b r (x r , y r ) = l r (x r , y r ) f r (x r ) ∏ q∈R\F r s∈F r m qs (x s ) (5.50)
where the product implies an implicit condition that regions q is a parent of region s, i.e. (q, s) is an edge of the region-graph. Quantity m qs is a message from q to s, closely related to Lagrange multiplier on the same edge such that:

m qs (x s ) = e -µ qs (x s ) (5.51)
However, values of Lagrange multipliers are not available, this is the dual problem of the optimization. Not solved in the literature, to our knowledge, we shall deal with an iterative solution, as with Bethe approximation.

We make use of the local marginalization constraint (C2) to extract an update equation of a message between two regions (r, s) linked by an edge. We get:

m rs (x s , y s ) = ∑ x r ∪x s l r\s (x r\s , y r\s ) f r\s (x r ) ∏ u∈R\F r v∈F r \F s m uv (x v , y v ) ∏ u∈D r \F s v∈D s m uv (x v , y v ) (5.52)
where:

l r\s (x r\s , y r\s ) = ∏ X i ∈X r \X s l i (x i , y i ) (5.53)
The message-passing that runs on the region-graph is a generalization of the BP algorithm that runs on the factor graph, thus it is called the Generalized Belief Propagation algorithm (GBPA).

Example

We consider the Tanner graph in Fig. 5.76 . We cluster it with squares of four check nodes, that provides us the region-graph in Fig. 5.8. According to equation ( 5.52), message from region r, with C r = {c b , c d , c e , c g } and X r = {X 2 , X 3 , X 5 , X 6 }, to region s, with C s = c g and X s = {X 5 , X 6 }, is:

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9
m 2356→56 (x 56 , y 56 ) = ∑ x 23 c b (x 23 )c d (x 25 )c e (x 36 )l 2 (x 2 , y 2 )l 3 (x 3 , y 3 )m 1245→25 (x 25 , y 25 ) m 25→5 (x 5 , y 5 ) (5.54) c a , c c , c d , c f X 1 , X 2 , X 4 , X 5 c b , c d , c e , c g X 2 , X 3 , X 5 , X 6 c f , c h , c i , c k X 4 , X 5 , X 7 , X 8 c g , c i , c j , c l X 5 , X 6 , X 8 , X 9 c d X 2 , X 5 c f X 4 , X 5 c g X 5 , X 6 c i X 5 , X 8 X 5

Complexity

For a given couple (r, s ∈ E r ), the number of elementary operations required to compute message m rs (x s ), for a given x s , is near 2 |X r |-|X s | . The total number of elementary operations to compute all messages on edge (r → s) is accordingly n eo (r) = 2 |X r | . Therefore, it is necessary to cluster any Tanner graph whose clusters are small enough such that n eo (r) is reasonable for any region r. This is the common trade-off of coding theory:

• too large clusters allow us to absorb and annihilate most loops but they induce a GBP with a very high complexity,

• too small clusters do not help reducing loops effects but they induce a GBP with a competitive complexity.

The GBP algorithm is not currently of a practical interest. To our knowledge, there is no paper in the literature that deals with its complexity. Update equations cannot be as easily reduced as BP ones, because messages are not two-dimensional vectors with binary arguments but multi-dimensional vectors. Nevertheless, the GBP algorithm must be investigated firstly because there are many of its properties that are still concealed behind update equations, then because it reflects the difficulty to solve inference on most factor graphs. Even though its implementation would ask for unreasonable hardware components, it provides us very relevant information about the effects of topology.

Region-graph reduction

Any loopy region-graph cannot perform optimally, only loopfree ones offer optimal approximation as marginal belief on any region r can be derived as:

b r (x r , y r ) = l r (x r , y r ) f r (x r ) ∑ x∪x r ∏ X i ∈X r l i (x i , y i ) ∏ c a ∈C r f a (x a ) (5.55)
that is precisely the a posteriori probability distribution p r (x r , y r ). In [P. 03],

Pakzad proposed a method to modify the region-graph removing few edges according to a specific rule. We call this method the Pakzad reduction of the region-graph. The method is based on the fact that both Lagrange functions (5.26) (5.46) contain few linearly dependent constraints. In optimization theory [D.P01] is stipulated that all constraints must be linearly independent to avoid redundancy. Pakzad then removed corresponding redundant constraints in Lagrange functions. As an example, we consider the small regiongraph in Fig. 5.9.

∀x s , ∑ x p ∪x s b p (x p ) = b s (x s ) (5.56) ∀x s , ∑ x q ∪x s b q (x q ) = b s (x s ) (5.57) ∀x t , ∑ x q ∪x t b q (x q ) = b t (x t ) (5.58) ∀x t , ∑ x r ∪x t b r (x r ) = b t (x t ) (5.59) ∀x u , ∑ x s ∪x u b s (x s ) = b u (x u ) (5.60) ∀x u , ∑ x t ∪x u b t (x t ) = b u (x u ) (5.61)
Constraints (5.57) and ( 5.60) imply that:

∀x u , ∑ x q ∪x u b q (x q ) = b u (x u ) (5.62)
However, constraints (5.58) and ( 5.61) also imply (5.62). As t ∈ P u , we affirm that (5.57), ( 5.60) and ( 5.58) imply (5.61). In other words, constraints (5.61) has to be removed from Lagrange function for it to satisfy the linear independence.

Removing local marginalization constraints between regions t and u means that the resulting GBP algorithm does not compute any message between them. The region-graph is then reduced to the one depicted in Fig. 5.10. We obtain a loopfree region-graph with an exact RBA i.e. with an optimal GBP algorithm. Unfortunately, Pakzad reduction hardly results in loopfree regiongraphs, e.g. we performed extensive experiments on Hamming, square lattice and Tanner codes, and reduction does not provide loopfree region-graphs. Simulations demonstrated that Pakzad reduction did not significantly improve RBA. Eventually, Pakzad reduction is a well-suited method for specific region-graphs, but not for all of them. However, it helps decrease computation time of GBP as it removes a few message computations. As it does not damage the performance but still improves computation time at least, we will use it in the following of the work.

Damping factor

A inherent problem of the GBPA is the lack of convergence it suffers from. To our knowledge, the literature does not provide any explanation concerning this phenomenon, but all experiments tend to prove it. In [J.S05] is introduced a damping on messages to avoid this problem, a well-known solution for many iterative algorithms. We denote by F rs the implicit update function that computes the value of m

(k)
rs for any k such that:

m (k) rs (x s , y s ) = 1 2 F rs + 1 2 m (k-1) rs (x s , y s ) (5.63)
As a matter of fact this equation is a uniform blend of memory m (k-1) rs and update F rs . There is still no visible reason to restrict the modification to a uniform damping, that is why we introduce a general equation:

m (k) rs (x s , y s ) = w k F rs + (1 -w k )m (k-1) rs (x s , y s ) (5.64)
where 0 ≤ w k ≤ 1 is the weight or damping factor. Currently, damping of the GBPA is not much tackled in the literature therefore we do not have many clues to determine properties it should verify. In [P. 03] is assumed that convergence of the GBP decoder is mostly ensured if w k converges to zero as k is getting larger. It contradicts (5.63) but as no general theory has been written, we cannot testify the veracity of any assumption. Investigations may only be made of experiments. This way, first idea is to use the simplest damping laws, that is to say we consider that w k is constant, w k = w 0 for any k. This assumption is against the one proposed in [P. 03] because it does not converge to zero but it is consistent with [J.S05]. After that, second idea is to use decreasing laws of w k , which match with assumption of [P. 03] and not really with (5.63). It is a hard task to extract a global theory of a damped GBP, a very large set of codes should be used to highlight inherent properties of the decoder. In our work that is not dedicated to this problem, we have tested our own different damping laws on Hamming and Tanner codes. We highlight an important assumption that w k does not depend on messages, that is not told in the literature whereas it is not so obvious. As an example, a very interesting work on a reweighted BPA (RBPA) has been presented in [T.G08]. To each edge of the factor graph is assigned a particular weight, according to a particular method based on spanning trees, detailed in the paper. This improves convergence of the decoder. In [H. 12] and [H. 11] is introduced that a uniformly RBPA (URBPA), i.e. a BPA with all edges equally weighted by a value ρ. To extract the weight that offers the lowest BER, the only method is to run decoder for many values of ρ. Displayed in Fig. 5.11 for an LDPC code of length N = 256 with M = 128, results show that the URBPA may significantly surpass the BPA for specific values of ρ. This experimental protocol provides the optimal weight, then we use it in our study to select the constant damping factor.

Constant damping factor

We consider that any message m rs is updated according to the following equation: m takes the advantage over update F rs whereas for high values of w 0 it is the other direction. We display in Fig. 5.12 the BER evolution as a function of w, on Tanner code. Results are averaged over 1000 simulations with 1000 iterations maximum. First of all, we observe that value w = 1.0, i.e. pure GBP algorithm without any blend with memory, provides the worse BER. As a significant example, we see that BER at w = 0.6 for the three SNR values is three times less than BER at w = 1.0. After that, it turns out that the most relevant weight is given by the evolution of the number of iterations K before GBP convergence or appearance of a null syndrome, represented in Fig. 5.13. Indeed for the three SNR values, the shape of K is similar, K reaches a minimum when w ≈ 0.675. We represent in Fig. 5.14 the BER as a function of the SNR for three values of w. We see significant improvement in the BER, for a damping w = 0.675, as much as the number of iterations K. Coming back to the figures Fig. 5.12 and Fig. 5.13, a noteworthy observation is that low values of w do not provide the worst BER, even though corresponding values of K are very high. Furthermore, BER is quite unchanged when w ≥ 0.675. In other words, pure update F rs is less accurate than memory. This leads us to consider a varying damping factor with k, that offers smoother relations between update and memory. We display in Fig. 5.15 and Fig. 5.16 the BER and the number of iterations of GBP, respectively, on Hamming code with a varying damping factor. We observe that the minimum of the BER is reaches for w ≈ 0.85, that also corresponds with the minimum number of iterations K. Contrary to Tanner code, there is no obvious evidence that would validate the assumption for the Tanner code that stipulates that pure update is less accurate than memory. 0 0.2 0.4 0.6 0.8 1 6.62

Decreasing damping factor

Now we consider that the damping factor follows a non-constant evolution that converges to zero. We selected three different profiles depicted in Fig. 5.17:

• the parabolic factor (P) of order n is:

w k = w 0 (1 -k K ) n ,
• the inverse parabolic factor (IP) of order n is:

w k = w 0 1 -( k K ) n ,
• the affine factor (A) is: w k = w 0 -k K w 0 , The affine factor (A) is only a particular case of the parabolic one, or the inverse parabolic one, with n = 1. The key point of w k is the balance between memory m (k-1) rq and update F rq . For functions we introduce, this balance is different from case to case. In the case of a parabolic function, memory is quickly favored at the expense of update, whereas in the case of an inverse parabolic function, memory is almost ignored for a long time. The affine damping coefficient offers a quite balance distribution, in terms of k, between memory and update. On figures Fig. 5.18 are displayed the BER on Tanner code for the different damping laws with an SNR = 2.90dB and order n = 2 of P and IP factors. We observe that all the laws, except the constant one, allows GBP to converge to the same BER, the difference remains in the convergence speed. The fact is that the parabolic coefficient provides the fastest GBP algorithm, as it is significantly visible for high values of w 0 . This property is quite surprising because it means that update functions of GBP has to be ignored as quickly as possible for the benefit of memory. This process seems unnatural because any decoding algorithm is aimed at keeping control on information that pass along edges of the associated graph. Finally, it means that the GBP algorithm with the parabolic damping can be considered as a "starter" decoding algorithm.

On Fig. 5.19 we show the BER for Hamming code with a parabolic damping factor corresponding to n = 2 and w 0 = 0.855 and the BER any without damping. It turns out that the original version is even worse than what the BP algorithm provides. On the contrary, the parabolic damping factor helps the GBP algorithm be more accurate. It should be recommended to carry the same investigation for many other LDPC codes to extract the best damping law. It appears in the current case that the parabolic damping factor with an initial value greater than 0.5 is a very good candidate to improve GBP performance.

Iteration index

Equation (5.52) is not practical because it does not carry any iterative data. According to local marginalization constraints, for any couple of regions r ∈ R, s ∈ E r :

b s (x s ) = ∑ x r ∪x s b r (x r ) (5.66)
that means that left side is computed by means of right side of the equation. More precisely, this equation can be rephrased as a sequence:

b (k+1) s (x s ) = ∑ x r ∪x s b (k) r (x r ) (5.67)
where k is an iterative index. Thus equation (5.52) is practically:

m (k+1) rs (x s , y s ) = ∑ x r ∪x s l r\s (x r\s , y r\s ) f r\s (x r ) ∏ u∈R\F r v∈F r \F s m (k) uv (x v , y v ) ∏ u∈D r \F s v∈D s m (k+1) uv (x v , y v ) (5.68)
It is not a strict update from k to k + 1 given that on the right side of the equation are left messages at iteration k and messages at iteration k + 1. From a practical point of view, it is advantageous to search for a way to parallelize the GBP as it would speed it up. At first sight, we could think of approximating (5.68) forcing the whole right side of the equation to be at iteration k. Doing so would allow us to implement a version where all messages at iteration k are computed in parallel, which would dramatically reduce the computation time. We show in Fig. 5.20 the BER of the GBP algorithm from the parallel version, denoted by GBP p , and the original one, also called the serial one, denoted by GBP s , on Tanner code with the same initialization. Without any ambiguity, GBP p offers poor results compared with GBP s , the BER is significantly damaged. Furthermore, except for SNR = 2.75dB, the GBP p algorithm does not tangibly converge whereas the GBP s algorithm quickly converges towards the final state it will end in. We observe for SNR = 3.20dB, which is a quite high SNR value for this code, that GBP p has not converged yet whereas GBP s has converged in less than 20 iterations.

Finally, experiments show that the parallelization of the GBP algorithm should not be done at once on the whole graph. Actually, equation (5.68) indicates the scheduling of GBP: it requires messages to be computed bottom up in the region-graph. Thus, the best choice is to compute messages generation after generation, with a parallelization inside each generation, beginning by the bottom of R. In this case, parallelization makes sense as we obtain same results as original GBP. § 5.5 Region-graph construction Constructing a loop free region-graph is not a trivial task. Nothing makes sure that a loopfree region-graph exists for a given Tanner graph7 . To our knowledge, there is no general method that provides the region-graph with the optimal GBP algorithm for a Tanner graph because there is no criterion that helps discriminate them.

Systematic construction

In [J.S05] is explained a systematic process that provides a region-graph for any Tanner graph, without using any criterion. Each cluster is made of a single check node accompanied by its neighborhood 8 , someway as in Bethe approximation. Regions of next generations are constructed as explained before, searching for common children. The advantage of this construction is that it is systematic, i.e. we may use it for any LDPC code. The region-graph in Fig. 5.21 is an example of this method applied on Hamming code. In spite of the ease of this construction, it does not offer a very relevant choice for the GBP algorithm according to the goal of the region-based approximation. Previously, we mentioned that it is aimed at absorbing harmful topological structures of the Tanner graph. Therefore, its construction has to be based on the shape of the Tanner graph, which is not the case here.

X 1 X 2 X 4 X 3 X 5 X 6 X 7 c a c b c c (a) Tanner graph c a , X a c b , X b c c , X c X 1 , X 2 X 1 , X 3 X 1 , X 4 X 1 (b) Region-graph

Square lattice

The square lattice Tanner graph in Fig. 5. 22(a) has the great advantage to be topologically regular. It is a pavement made with a constant square basis, then it is very easy to build such a code, even for any other length N. But squares are also the shortest loops, i.e. they are the most harmful topological structures for BP decoding. They accordingly appear to be natural choices for clusters, see Fig. 5.22(b). Absorbing short loops offers a better correction capability, even if the resulting region-graph is not loopfree9 . We display in Fig. 5.23 BER performance, representing pure GBP and uniformly damped GBP compared with the BP algorithm. We observe that even without blending update and memory, the GBP algorithm is more efficient than BP. As an example, at SNR= 5dB, the BER of BP is 2.10 -2 , whereas undamped GBP reaches a BER lower than 2.10 -3 and damped GBP has a BER lower than 2.10 -4 : there is ten times less errors between them, which is a significant result.

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 c a c b c c c d c e c f c g c h c i c j c k c l (a) Tanner graph c a , c c , c d , c f X 1 , X 2 , X 4 , X 5 c b , c d , c e , c g X 2 , X 3 , X 5 , X 6 c f , c h , c i , c k X 4 , X 5 , X 7 , X 8 c g , c i , c j , c l X 5 , X 6 , X 8 , X 9 c d X 2 , X 5 c f X 4 , X 5 c g X 5 , X 6 c i X 5 , X 8
One would suggest to use concatenations of squares as basic clusters but these clusters would be so large that computation complexity would skyrocket, which is not acceptable. Smallest squares succeed in fairly weighting good performance and competitive speed.

Trapping-sets

For most LDPC codes, the Tanner graph topology is not as simple as the square lattice one. It is then close to impossible to construct associated relevant region-graphs 10 . In contrast, in spite of our lack of knowledge on codes topology, significant advances have been made to bring out specific subgraphs, that are trapping sets, responsible for BP failures, see chapter 4 for more details. GBP may be consequently useful provided that region-graphs are constructed taking them into account.

Each code needs a long and deep study to extract a relevant region-graph with a GBP algorithm better than its BP equivalent, the set of trapping sets that compose a given Tanner graph being not always available. In this thesis, we focus on a specific code, the Tanner code of length N = 155 with check nodes and variable nodes degrees d c = 5 and d v = 3, respectively. The reason to choose it is that its Tanner graph is entirely covered by a set of 155 trapping sets TS (5,3), see Fig. 5.24, that is available 11 . Therefore, clusters that build the first generation of the region-graph are naturally selected. 10 notion of relevancy is hard to define, goal being to absorb harmful topological structures 11 the fact that N is also the number of TS (5,3) is not useful in our study

X 1 X 2 X 3 X 4 X 5
The best way to neutralize TS (5,3) is to absorb each trapping set in a cluster. Problem is that it is absolutely not practical. We recall a region-graph construction rule that any check node c a ∈ C, included in an arbitrary region r, has to be accompanied by its neighborhood X a ⊂ X. Counting the number of neighbors of all check nodes inside a TS (5,3), we obtain that such a cluster would contain 9 check nodes and 35 variable nodes. According to the summation in equation ( 5.68), any message from a cluster would then need at least 2 35 elementary operations, which is not acceptable. Eventually, trapping sets need splitting to create clusters of practical size.

Local optimality

As it is very hard to ensure optimality on the whole region-graph, we orient a clustering by considering local optimality. This way, we introduce a basic rule for a region-graph construction.

Local loopfree principle: the local region-graph, made of clusters and their children, created splitting a harmful topological structure, has to be loopfree such that the associated local RBA is exact.

Following this principle, we perform a clustering such that local region-graphs made by splitting trapping sets are loopfree. Unfortunately, when considering all trapping sets, the whole region-graph is not loopfree anymore12 , then RBA remains an approximation. If we do not respect the local loopfree principle, RBA, and so GBP, are likely to provide very bad performance. For a given TS (5,3), we show as examples two different splittings so that a cluster may be:

• either the subgraph made by one of the three shortest loops with two unverified connected check nodes, see Fig. 5.25 where concerned nodes are colored (cluster is called quadruplet):

variable nodes X 1 , X 2 , X 3 , X 4 and check nodes c a , c b , c d , c e , c g , c h , • either the subgraph made by a triplet of variable nodes and one unverified connected check nodes, see Fig. 5.26 where concerned nodes are colored in blue (cluster is called triplet):

-variable nodes X 2 , X 3 , X 4 , X 5 and check nodes c b , c c , c e , c f , c h , c i , -variable nodes X 1 , X 2 , X 4 , X 5 and check nodes c a , c c , c d , c f , c g , c i , X 1 X 2 X 3 X 4 X 5 c a c b c d c c c e c f c g c h c i (a) First cluster X 1 X 2 X 3 X 4 X 5
variable nodes X 1 , X 2 , X 4 and check nodes c a , c d , c g , -variable nodes X 2 , X 3 , X 4 and check nodes c b , c e , c h , -variable nodes X 2 , X 4 , X 5 and check nodes c c , c f , c i . Splitting a trapping set in quadruplets fulfills the local loopfree principle, see Fig. 5.27. However, there is no balance between all the variable nodes as a TS (5,3) contains three quadruplets and not only two, see the third cluster in Fig. 5.28.

X 1 X 2 X 3 X 4 X 5 c a c b c d c c c e c f c g c h c i (a) First cluster X 1 X 2 X 3 X 4 X 5 c a c b c d c c c e c f c g c h c i (b) Second cluster X 1 X 2 X 3 X 4 X 5
c a , c b , c d , c e , c g , c h X 1 , X 2 , X 3 , X 4 c b , c c , c e , c f , c h , c i X 2 , X 3 , X 4 , X 5 c b , c e , c h X 2 , X 3 , X 4
Figure 5.27: Region-graph resulting from the split of a TS (5,3)

in two quadruplets

The whole region-graph is not balanced too, but between all TS (5,3). As a matter of fact, all variables nodes equally participate to several trapping sets 13 . The split in two quadruplets then involves that some variable nodes are included in more clusters than other variable nodes. This makes the GBP algorithm biased in the sense that noise realizations, and then likelihoods, are not fairly introduced to initialize the region-graph. More precisely, variable nodes that appear in a lot of clusters would wield more influence than variable nodes included in less clusters, whereas nothing either in the Tanner code or the channel discriminates between them. Finally, it turns out to be impossible to determine any criterion to select either the first and the second clusters, or the first and the third clusters, or the second and the third clusters. Eventually, we should consider the three quadruplets to build the local regiongraph. However, this would violate the local loopfree principle, as the regiongraph is not loopfree, see Fig. 5.29. Furthermore, quadruplets contain 24 variable nodes that is quite too large to consider a GBP algorithm of reasonable complexity, as a message from a cluster would need 2 24 = 16777216 elementary operations.

X 1 X 2 X 3 X 4 X 5
Embedded inside the whole Tanner graph, any TS (5,3) has many other connections than the edges exhibited in Fig. 5.24. The complete trapping-set in this case is Fig. 5.31. Empty circles in this figure are outer variable nodes, supposed to make the parity-check equations satisfied without assignment to the a inner variable nodes. Construction of the region-graph then results in a loopy Bayesian network, whose first generation is made of the clusters that we have created, the triplets. (5,3) with all variable nodes connected to its check nodes A noteworthy property of the triplets construction is that the two last generations are exactly generations of the Tanner graph:

X 1 X 2 X 3 X 4 X 5 c a c b c c c d c e c f c g c h c i Figure 5.31: A TS
• any of the 93 regions of the penultimate generation R L-1 is only made of one check node and its neighborhood of variable nodes14 , L being the number of generations (3 in the current case),

• any of the 155 regions of the last generation R L is only made of a single variable node.

As a matter of fact, this particular point makes our constructions relevant according to the work made in [M. 04]. In this paper, Welling assumed that a relevant region-graph, i.e. a region-graph whose GBP algorithm performs well, has to be constructed "above" the corresponding Tanner graph 15 : bottom generation R 3 is made of regions each containing a single variable node such that all variable nodes are represented, any region of the upper generation R 2 is made of a check node in such a way that all check nodes are represented, and regions of the upper generation R 1 , called super nodes, are clusters that would allow us to break or absorb harmful topological structures of the two last generations. Our triplets construction satisfies the Welling assumption as we break inner loops of trapping sets.

The Welling assumption is perfectly in agreement with our local loopfree principle provided that super nodes are added to Tanner graphs, and not to region-graphs. As an example, one would note that in the quadruplets construction, in Fig. 5.29, children of quadruplets are triplets, i.e. quadruplets are super nodes of triplets. We recall that a balanced region-graph made with quadruplets contains loops, that means that adding super nodes to the regiongraph in triplets damage its local optimality. Finally, the Welling assumption is relevant for Tanner graphs only.

Experiments

Considering that quadruplets construction implies a very large number of elementary operations, we decide not to hold it for computation time's sake. Experiments we have performed accordingly only concerned triplets construction. All previous results concerning the damped GBP, the iteration index for the scheduling were realized on Tanner code with triplets construction. We now present other comparisons between BP and GBP. Usually, decoders performance are delivered exhibiting the averaged BER as a function of the SNR. Our own experiments, not displayed here, showed that BER of BP and GBP were merged for low and middle values of the SNR on AWGNC 16 . This is an encouraging result as it confirms that the special construction of the regiongraph with triplets does not damage the BER, in average 17 . As the triplets construction is based on trapping sets, and as trapping sets are known to damage BP decoding for weakly noisy channels, performance of GBP will distinguish themselves from BP ones for high SNR values on AWGNC and low p values on BSC.

As a channel noise is lowered, the number of simulations needed to get a trustworthy BER is increased. We previously mentioned that GBP suffers from a large computation complexity compared with BP one, which prevents us from computing any average BER in this case. Instead of considering random noise realizations, we then orient experiments on channel realizations that are deeply linked with trapping sets.

BSC

We deal with low-weight error events 18 on BSC. We only take into account specific low-weight error events that prevent BP from correctly decoding because of TS (5,3). To keep a reasonable amount of results, we select a single error event, EE, that fairly represents all error events that deals with TS (5,3). We roll out an experimental protocol for a given crossover probability p:

• for any variable node X i that is not in EE, we force the channel to copy the input:

y i = x i ,
• for any variable node X i that is in EE, we force the channel to flip the input: y i = xi . 16 and of the crossover probability p on BSC 17 local loopfree principle is an assumption without proof, then experiments are the only means to check it 18 a low-weight error event is a codeword that BSC has slightly disturbed, i.e. very few bits have been flipped

EE is made with 7 errors, a low value compared with the codewords length N = 155. This makes our protocol realistic as getting 7 corrupted bits is much more encountered for low p values than for large ones. First of all, for p = 0.2, GBP BERs do not oscillate as strongly as BP except for GBP-C which local maximum though decreases along k contrary to BP: local maximum values are 0.34, 0.26, 0.25, 0.24, 0.23 for k = 21, 39, 57, 75, 93, respectively. These pseudo-oscillations may be hoped to end to a very low BER, but experiments, not presented here, showed that it is not true, BER of GBP-C always oscillates. GBP-IP seems to behave in the same way, but as for GBP-P, the decreasing law of w k makes the algorithm converge to a single value. This helps local maximum values be dramatically reduced compared with GBP-C: we obtain 0.38, 0.21, 0.14 for k = 21, 40, 60. GBP-P exhibits a small bound around k = 20, as other decoders, but with a lower value 0.1096774. Then during thirty iterations it keeps a constant value unfortunately increased around k = 70. This phenomenon emphasizes that our triplets construction is dedicated to situations where p is significantly low. If p is too high then EE is not realistic, and the experimental protocol too.

As p is getting smaller, behaviors of GBPs merge into a single behavior. All of them converge faster and faster to the same steady state as p decreases20 . For p = 2.10 -2 and p = 2.10 -3 , we observe bounds of no significance for GBP-C and GBP-IP around k = 40, as they disappear for lower values of p. For p = 2.10 -5 , no more than 20 iterations are needed for all GBPs to converge, and BERs reach a value of 0.04, which corresponds to 8 errors out of 155 bits. We may interpret this result telling that GBPs did not succeed in obtaining as many corrected bits as BP could do21 , in addition it has added an error to EE. From another point of view, when comparing both iterative evolutions, we easily see that reaching 0.03 is close to a random process: BP oscillations are not perfectly predicted, then stopping the decoder to an exact local minimum makes the experimenter very lucky. On the contrary, all GBPs reach a unique stable and convergent evolution which allows us to ensure these decoders output, GBP is more reliable and trustworthy than BP.

Finally, we may assume that GBP is able to decrease TS (5,3) effects, which makes it able to decrease LDPC error floor on BSC 4.8.6. We cannot verify this assumption as computation complexity of GBP is too high to simulate it for very low p values, but the experimental protocol is a relevant substitute.

AWGNC

The experimental protocol on AWGNC is similar to the previous one, but the error event notion has to be redefined as it is originally defined on BSC. As a matter of fact, flipping bits does not make sense on AWGNC as noise is a continuous elements independently added on all bits. The concept of error event on this channel is treated using bits corruption: a bit X i of value x i is corrupted if likelihood, p i (y i |x i ), deduced from channel observation, y i , is maximized for xi . A low-weight error event EE is then defined as a codeword that AWGNC has slightly disturbed, i.e. only very few bits have been corrupted. Practically, on BSC, we forced a bit to be erroneous flipping it if channel did not. On AWGNC, we make a bit X i corrupted multiplying the channel output y i by -1: multiplying y i by -1 is equivalent to swapping likelihoods l i (0, y i ) and l i (1, y i ) 22 . As SNR is increased any error event is then more and more emphasized. The drawback of this protocol is that moments, as mean and variance, of random variable Y|X are modified, as depicted in Fig. 5.33. This is due to the fact that random realizations are discriminated to only get specific error events23 . However, as the same protocol is used for BP and GBP, comparing between both decoders still holds, we only have to keep in mind that SNR values on next figures are not SNR f but SNR i , i.e. not estimated ones but original ones. We display iterative evolutions of the BER of BP and GBP-P in Fig. 5.3424 . We specify that GBP-IP and GBP-C are not represented in this figures as they behave in the same way that GBP-P, which means that results stem from the triplets clustering, not from damping. For low values of SNR i , the BP BER pseudo-oscillates up to a steady state of 0.13 at 5dB, that is a substantial value compared with GBP one that quickly reaches 0.06 at 5dB. As SNR i is decreased, the BER of BP stops oscillating by converging in less than 20 iterations to a non-negligible increasing steady value of 0.305 at 11dB, whereas the GBP BER converges to a dramatically lower value of 5.10 -2 at 11dB i.e. about six times less errors than BP.

For SNR i ≥ 13dB, BP does not converge anymore: it either oscillates, see at 14dB and 17dB, or it exhibits an erratic evolution, see at 15dB25 . We see that the range of values of the BER is large enough such that the choice of the last practical iteration widely influences the performance of the algorithm. On the contrary, GBP still convergences in nearly 30 iterations to the same value of 5.10 -2 for any SNR i value. Eventually, for SNR i ≥ 18dB, the BP BER reaches steady and decreasing values. For very high SNR i , it is stuck to the same non-zero value as GBP around 3.10 -2 . We summarize the steady values of the BER of both algorithms in Fig. 5.35. We observe that the triplets clustering provides a robust GBP algorithm against TS(5, 3) compared with BP. For very high values of SNR i , likelihood of any variable node X i of EE is almost:

l i (0, y i ) ≈ 0, l i (1, y i ) ≈ 1 (5.69)
therefore the corresponding bit is too corrupted to be correctly decoded by either GBP or BP. That is why for both decoders converge to a non-zero value of the BER, e.g. at SNR i = 20dB, we obtain a BER of 3.10 -2 for the two algorithms. As EE is of weight 7, such a value of the BER is equivalent to 4 errors, we conclude that BP and GBPs can only correct 3 errors. § 5.6 Conclusion

The Bethe approximation is a particular case of the more general region-based approximation, as the BP algorithm is the GBP algorithm when the regiongraph is exactly made of the Tanner graph 26 . We observed that generalizing the BP message-passing needs a generalization of the node definition, to create region-graphs which edges are support for the GBP algorithm. In addition, new damping rules for GBP help converge, as its original version is unstable, and obtain better results in terms of BER.

We saw that contrary to the Tanner graph, we can find a multitude of region-graphs for any LDPC code, since the clustering is not unique. We then introduced the local loopfree principle and an associate specific construction based on particular clusters, the triplets, stemming from the split of trapping sets. GBP is then aimed at dealing with error events deeply connected with trapping sets. Exhibited simulation results, according to a dedicated protocol, showed that the triplets clustering provided a GBP algorithm more accurate than BP in these cases. These results tend to investigate more on the GBP algorithm, and particularly the construction of the region-graph. This would help bring out other clustering rules and new principles to better understand GBP and RBA. In addition, the study of the damping factor is not closed, we have raised an open problem.

The specific construction in triplets is a technique that can be used as a basis to cluster any LDPC code, provided that known and harmful topological structures entirely cover its Tanner graph. In the literature, clustering the square grid code using squares follows the same idea, and we also observed that GBP outperformed BP in this case. To our knowledge, no other clustering has been proposed on any other code, therefore our construction is likely a new lead in the research field on graphical models.

-Chapter 6 -Dynamical behavior of iterative decoders § 6.1 Introduction Messages m (k) of BP and the GBP are computed at any iteration k, according to update rules (1.12), (1.13), and (5.52). These rules are functions:

• all messages at the previous iteration {m (k) e } e∈E , where E is the set of edges either of a Tanner graph or a region-graph,

• initialization parameters: SNR and likelihoods.

It is accordingly legitimate to consider update rules of these algorithms as iterated maps. This way, it is possible to describe BP and GBP with tools of dynamical systems theory. It allows us to introduce explanations of non-trivial phenomenon, e.g. non monotonous evolution of the BP convergence rate exhibited in 1.7. The goal of this chapter is to extract a few typical behaviors of BP and GBP and to examine them as functions of SNR. More specifically we want to emphasize:

• the variety of behaviors encountered with these decoders,

• how these behaviors wield influence on decoding performance.

First of all, this work needs the definition of a mathematical framework, known as the state space. After that, we will present dynamical quantifiers to describe the decoders. Among the large set of estimators commonly encountered in the literature, we chose a few of them to characterize BP and GBP dynamics. This investigation is aimed at highlighting particular SNR values for which their behaviors are non-trivial, describing the nature of these behaviors and their connections with decoding performance. This would help understand and predict, or even control, decoders behaviors. All experiments in this chapter are performed on Tanner code, to specify dynamics of GBP coupled with the triplets construction. In addition, from now on, the channel used is only the AWGNC.

A Ratio between both variances is around 10 2 , i.e. BER may change ten times more around point B than around point A, which is a considerable difference. One would link this result to the specific couple (X i = 100, X j = 102) but it turns out possible to observe such a behavior for other couples and noise realizations 1 . It then appears that message-passing algorithms are strongly dependent on initialization, i.e. when slight changes are induced.

The likelihoods space provides a BER map for a given couple of decoder and LDPC code. Such a map would give access to the whole description of the error correction capability 2 . Unfortunately, to compute a complete map is needed not only two likelihoods but N likelihoods, N being the number of variable nodes. Even for low values of N 3 , computing the BER map is then impossible, as it would need A N configurations with A the number of values to assign values to a single likelihood. Furthermore, storing this map would ask for unrealistic hardware requirements. Using less informative but tractable estimators finally turns out necessary to understand decoders behaviors. § 6. 3

State space

The study of any iterated map needs state space definition, i.e. a mathematical space where state variables evolve. Usually, state variables are defined as variables handled by the iterated map, e.g. the logistic map transforms x (k) in x (k+1) by function f (x) = ax(1x) therefore x (k) is the state variable. In the same way, update rules of decoding algorithms may be summarized as iterate maps according to the following rules:

BP : ∀(X i , c a )an edge, m (k) ai = F BP ({m (k-1) bj } b,j , {l j } j ), (6.1)
GBP : ∀(r, s)an edge, m

(k) rs = F GBP ({m (k) uv } u,v , {m (k-1) uv } u,v , {l j } j ), (6.2)
1 exhibiting all simulation would not be reasonable to keep a fair paper length 2 whether it is the code or the decoder 3 Tanner code of length N = 155 is considered in coding theory as a small code where update rules F BP , F GBP are provided in previous chapters by equations (1.12),(1.13) and (5.52). F BP maps exactly messages from an iteration to the next one whereas F GBP computes kth iterated messages by former and current messages. Furthermore, BP messages are two-dimensional functions while GBP messages are multi-dimensional. Comparing both algorithms involves that iterated maps and state variables used to describe them are comparable, which is not the case here. To confront both decoders to same situations, we need common quantities.

Beliefs on single variable nodes are common quantities between BP and GBP. We then choose them to stand for state variables. In the following, we consider that any belief of a node X i refers only to b i (0, y i ) 4 . Expanding explicit update rules that compute beliefs as functions of former iteration beliefs, either for BP or GBP, proves to be a hard task therefore we do not present it. Nevertheless, we assume that such functions G BP , G GBP exist, which allows us to define their related dynamical systems:

• the dynamical system related to BP is the state space built from beliefs on the N single variable nodes computed from G BP ,

• the dynamical system related to GBP is the state space built from beliefs on the N single variable nodes computed from G GBP .

We note that in the case of the GBP algorithm, the bottom generation R L of the region-graph is not always made of single variable node regions, sometimes it does not even contain any single variable node region. For such cases, beliefs on missing single variable nodes are computed using the local marginalization on regions of the upper generation R L-1 , or still upper if necessary. Without any loss of generality, we now refer to any decoder and its related dynamical system by the same denomination, in order to lighten writings 5 .

We now have to choose behavior descriptors for decoders. BER is quite poor according to the quantity of information it offers as it firstly maps beliefs, i.e. continuous values in [0; 1], to estimated bits, i.e. discrete values in {0, 1}, and then it associates an average value over these bits. As an example, considering an arbitrary bit X i , an indecisive belief b i (0, y i ) = 0.49, b i (1, y i ) = 0.51 is mapped to xi = 1 in the same way as the belief b i (0, y i ) = 0.01, b i (1, y i ) = 0.99. It is a very rough estimate of the information contained in the belief. Using BER finally deletes a large part of information lying on graph edges, either Tanner graphs or region-graphs. The only advantage of BER is that it delivers any decoder performance with only one value, that makes it easy to store. In this investigation, searching for other descriptors appears absolutely necessary to obtain accurate information on decoders behaviors. § 6.4 Mean square beliefs

Definition

The mean-square-belief (MSB) introduced in [X. 05] is an interesting alternative to BER, as it offers a good compromise: it is a one-dimensional map that reflects decoding performance using directly beliefs. We may compute it at any iteration k for a given channel realization as:

E(k) = 1 N N ∑ i=1 b (k) i (x i ) 2 (6.3)
where x = [x 1 , . . . , x N ] is the emitted codeword, e.g. the null codeword that we always use in our experiments. This function has the following properties:

• E(k) = 1.00: correct decoding,

• E(k) = 0.25: any variable node entropy is maximum, no information is available on the emitted codeword,

• E(k) = 0.00: faulty decoding.

This function implicitly depends on SNR and likelihoods. MSB provides smoother information about a given decoder than BER as it does not include any threshold. Hence, it is a more faithful estimator.

Experimental protocol

MSB reveals performance of a given decoder according to a specific noise realization. To cover a large range of SNR values, we firstly generate N noisy samples n = [n 1 , . . . , n N ], all computed from a Gaussian law N (0, 1). After that, these samples are scaled on the SNR multiplying them by the associated standard deviation 6 . This way, we obtain a sequence of channel realizations:

y 1 = x 1 + σ 1 n (6.4
)

y 2 = x 2 + σ 2 n (6.5
) .7) where y j = [y j,1 , . . . , y j,N ] and x j = [x j,1 , . . . , x j,N ]. Samples y provide us likelihoods that initialize decoding algorithms. In order to obtain a consistent investigation, we use the same noise realizations for BP and GBP. We recall that the goal of this study is to focus on non-trivial decodings, i.e. decodings where the BER does not converge to any steady value7 . To our knowledge, there is no results in the literature that help find such realizations.

y 3 = x 3 + σ 3 n (6.6) . . . ( 6
The method we used in this thesis is to run the desired decoder for many simulations, each one with a different basic noise realization n. Plotting MSB for all of them allows us to discriminate non-trivial error events. In the following of this chapter, we use a specific error event EE 8 to exhibit dynamical properties of the dynamical systems related to BP and GBP. This error event perfectly reflects all typical behaviors that we encountered with other error events in our experiments, not displayed here, i.e. it is a reliable error event to describe non-trivial evolutions of decoders.

Attractor

We now present essential definitions from [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF]. First of all, the trajectory of a decoder is the sequence of points:

T = [[b (1) 1 . . . b (1) N ] T , . . . , [b (k) 1 . . . b (k) N ] T ] (6.8)
of the related dynamical system in its state space as the iteration k is increased.

In addition, an attractor is a set of points, inside the state space, to which the system, related to any decoder, approaches as the number of iterations goes to infinity. Moreover, the basin of attraction of a particular attractor is the set of initial conditions, i.e. likelihoods, that give rise to trajectories approaching the attractor as the number of iterations is increased to infinity.

The difficult task in the study of decoders dynamics is that the state space is so immense that we can list neither all attractors nor their basins of attraction. This is why scanning a large part of the state space is required to extract typical non-trivial error events.

Experiments

BP

We display in Fig. 6 A slight difference remains in transient phases that last longer. Moreover, GBP-IP eventually converges 9 . In contrast, the GBP-P algorithm does not exhibit such a vicinity. We do not observe clear oscillations at SNR = 2.30dB, even though values of E(k) are similar in average. In addition, it reaches an imperfect steady state at 3.09dB. However, as GBP-IP, it is due to converge in the last iterations by means of the damping factor. More generally, we notice that GBP behaviors, either damped or not, can be categorized in the same way that BP behaviors, see 6.1. Therefore, we use the critical SNR values {2.00, 2.30, 2.90, 3.09}dB in the following of the study to describe decoders dynamics. The investigation is now oriented to highlight SNR intervals in which typical behaviors live. § 6.5 Bifurcation diagram

Definition

Analyzing decoders behaviors for any error event and for any SNR value is a substantial work as it results in a tremendous amount of information. Computing the bifurcation diagram is a lead towards a summary of information. This descriptor is a one-dimensional function that provides the evolution of the MSB steady state as a function of the SNR, for any decoder. Strictly speaking, steady states are not systematically found, as MSB often oscillates, see previous figures as examples. We accordingly define the steady state, at an arbitrary SNR, as the MSB value at the last iteration K 10 . The goal is to bring out SNR values for which decoders radically change in their behavior. These changes are called bifurcations [K.T96]. An important practical point is that the bifurcation diagram shape depends on the last iteration K11 . We display in Fig. 6.7 this dependence on K computing the bifurcation diagram, using Tanner code, for four values of K: 30, 100, 1000, 10000. We observe that SNR intervals do not change, but increasing K makes their bounds emerge more significantly. For K = 30, it is hard to split SNR values into five intervals, E(K) does not change abruptly, contrary to the cases in which K ≥ 100. As K is increased, we may distinguish a Hopf bifurcation [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF], at SNR=2.20dB, that brings out two steady states. In other words, increasing K makes BP closer to its steady behaviors, choosing low K values limit BP to its transient phases.

Experiments

We recall that BP is practically run for less than one hundred iterations, i.e. the BP algorithm is not fairly used. As a matter of fact, the lack of convergence may often be solved considering larger K so that BP can reach a steady behavior. Within the state space, the BP dynamical system takes a very long time before falling into an attractor12 . For K = 100004, the attractor is obvious as the fork of the Hopf bifurcation is clearly visible, even if transient points still appear between the two fork branches. In addition, an extra behavior arises, for K ∈ {1000, 10000}, in interval [2.50 dB;2.54 dB], where E(K) stops oscillating but does not begins an erratic evolution, it seems to start up again a smooth increasing as for first interval. In spite of this noteworthy result, we shall not consider such an extra behavior, as selecting a large value of K takes us away from any practical interest. Thus, we use K = 100 to stay close to realistic situations, at the expense of the decoders long-term reliability. Constant damping triggers a behavior similar to the BP one as E(K) pseudooscillates before presenting erratic evolutions. GBP-C firstly differs from BP through the fact that a bifurcation distinctly occurs at SNR = 2.6dB whereas it is less visible for BP 13 . Secondly, GBP-C presents convergence jumps from SNR = 2.87dB, and GBP-IP too, which is lower than for BP as it has to reach larger SNR values than 2.99dB. We note that GBP-P behavior is slightly different: pseudo-oscillations take place earlier than BP and other GBPs, no visible bifurcation occurs around 2.6dB and it erratically behaves only from 2.8dB, i.e. 0.2dB later than other ones.

We note that for SNR values in [2.6 dB; 2.87 dB], BP and damped GBPs with constant and inverse parabolic laws present peculiar evolutions of E(K). The shape of bifurcation diagrams is close to the logistic map one: a cloud of points without particular consistency. It is known [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF] that the logistic map trajectory falls into a chaotic attractor for a ≥ 3.44, see Fig. 6.9. We accordingly assume that the erratic evolution of E(K) is chaos for both message-passing algorithms. Eventually, the bifurcation diagrams of GBPs show the same properties that we brought out for BP on the same error event EE. All of these decoders are similarly influenced by same likelihoods. In other words, the distinction between them is not easy to make, based on their bifurcation diagrams only14 .

The bifurcation diagram underlines critical SNR values and associated intervals in which decoders have particular behaviors. In addition, it helps us introduce assumptions concerning behaviors nature, that we now attempt to describe in depth. § 6.6 Reduced trajectory Iterative BER and MSB are means to draw decoders evolutions according to SNR values. Visualizing bifurcations and peculiar behaviors allows us to conjecture nature of associated attractors that trap decoders. Usually, if state spaces are one, two or three-dimensional, nothing prevents us from observing trajectories evolution, e.g. logistic and Hénon map, Lorenz system lie in humanly observable state spaces of one, two and three dimensions, respectively.

In the current study, trajectories of decoders lie in N-dimensional state spaces, with an often very large N 15 . Visualizing decoders trajectories is then absolutely impossible. BER and MSB are useful to pre-investigate on decoders behaviors but they are too restrictive as they are average of N beliefs. Another method is therefore essential to observe and analyze decoders trajectories.

State space reduction

We introduce an alternative space which number of dimensions is considerably lower than for the original state space. We build it through the State Space Reconstruction (SSR), exposed in chapter 3.5. We recall that SSR helps faithfully recover a state space only with a subset of its state variables. We interpret this property such that a subset of all state variables is enough to fairly describe any dynamical system. In other words, SSR allows us to construct a space of an arbitrary number of dimensions m. This way, we are able to construct three-dimensional reduced state spaces such that dynamical systems related to decoders are now humanly observable.

Decoders state spaces are generated by beliefs assumed to be state variables. We cannot use any beliefs subset to construct a reduced space because we do not have any criterion to discriminate three beliefs among N. One would suggest to take into account only "very variant" beliefs but this would mean that other parameters (SNR, EE, LDPC code) would discriminate three other beliefs, making reduced state spaces hardly comparable.

Therefore, we have to find a way to gather all beliefs without any discrimination in three quantities. As a matter of fact, this is perfectly possible using MSB. We practically compute a sequence, called time series, of K points [E(1) . . . E(K)] that we map to a matrix such that each column is a state variable of the reduced state space:

T =           T 1 T 2 . . . T n . . . T L           =           E(1) E(1 + τ) . . . E(1 + (m -1)τ) E(2) E(2 + τ) . . . E(2 + (m -1)τ) . . . . . . . . . E(n) E(n + τ) . . . E(n + (m -1)τ) . . . . . . . . . E(L) E(L + τ) . . . E(L + (m -1)τ)          
Vector [T 1 . . . T L ] contains L points of a trajectory, called reduced trajectory, inside the reduced space. The value m is the target dimension of the reduced space, called the embedding dimension. Here, m must be three to get a humanly observable space. The term τ is a delay and L is chosen such that:

L + (m -1)τ ≤ K (6.9)

Choice of τ and m are known in the literature to condition the SSR success, research on this topic is still on progress [H. 06]. We set the embedding dimension m = 3, because we need a three-dimensional space. Practically, the number of points K is not very large, around a hundred, to keep realistic decoders whereas a faithful reduced trajectory needs a large number of points L.

We accordingly choose the delay τ as small as possible, i.e. τ = 1. Eventually, we obtain the 3-dimensional trajectory of K -2 points:

T =           T 1 T 2 . . . T n . . . T K-2           =           E(1) E(2) E(3) E(2) E(3) E(4) . . . . . . . . . E(n) E(n + 1) E(n + 2) . . . . . . . . . E(K -2) E(K -1) E(K)          
We display in Fig. 6.10,Fig. 6.11,Fig. 6.12 and Fig. In the case SNR = 2.00dB, we observe the same convergence in spiral to the same fixed-point for any of the decoders. This fixed-point corresponds to what is called a faulty decoding in the study of the MSB. We increase the SNR to 2.30dB meaning that we enter the second SNR interval, where decoders behaviors are distinguishable. We observe that BP and GBP-C clearly oscillates 16 . Associated limit cycles look like ellipses, which axes, either minor or major, are oriented in the same direction for both of them. The only difference is that the GBP-C limit cycle is larger than that of the BP, i.e. oscillation of GBP-C are more important. Considering GBP-P and GBP-IP, we recall the damping helps them converge, making observations of limit cycle impossible, both end in fixed-points. Associated trajectories do not exhibit any obvious convergence as in the previous case SNR= 2.00dB. In addition, we perceive similar continuous lines inside their clouds of points, as for GBP-C and BP. These lines stem from the fact that damping, either parabolic or inverse parabolic, makes GBP eventually leave a "transient" limit cycle after a certain number of iterations k, as we also saw in Fig. 6.4. Limit cycles are oriented in the same direction, as for BP and GBP-C, towards the all-one point, i.e. the perfect decoding state. More precisely, ellipses of all decoders are close to the straight line, denoted by ∆, that goes from the all-zero point (completely faulty decoding) to the all-one point (perfect decoding), which is the unitary linear function in the reduced state space. This means that all decoders do not discriminate between reduced state variables, they are somehow equivalent. In other words, MSB smooths beliefs such that EE do not strongly affect the reduced trajectory.

We now increase the SNR to 2.90dB. Decoders behavior appears chaotic according to the bifurcation diagram conclusions. We do not notice any particular line or known geometric form, clouds of points are themselves chaotic attractors: two close points in such an attractor are mapped by any update function to separate points that belong to this chaotic attractor. However, we observe that GBP has the advantage of stretching the attractor along ∆ line, contrary to BP where the attractor shape is unchanged from any state space point of view. GBP tends to make state variables "more equivalent" than BP, i.e. EE is better smoothed by GBP than by BP, which is a clear advantage of the GBP algorithm. Convergence jumps of GBP-P and GBP-IP confirm that the stretched form of the chaotic attractor helps GBP correctly and quickly converge. When SNR= 3.09dB, we observe that all decoders, whatever the damping factor for GBP, eventually converge to the perfect decoding fixed-point. Furthermore, the number of iterations needed to converge are clearly fewer for GBP than for BP, making GBP consistent with the previous observations concerning the convergence jumps.

The method of the reduced trajectory allows us to visualize and confirm nature of typical decoders behaviors. A numerical method now turns out necessary to associate to any attractor intensity and size, to relate them to decoding performance. § 6.7 Sensitivity to initial conditions Reduced trajectories and bifurcation diagrams bring out that decoders fall into attractors of different kinds, depending on the SNR. In the interval [2.51 dB, 2.99 dB], related dynamical systems oddly behave as attractors, in reduced spaces, do not present any regular geometric shape. This observation is an evidence that decoders do greatly depend on initial conditions, i.e. on likelihoods. In other words, decoders reveal stability: a dynamical system is said stable if and only if it is not strongly affected by slight changes in its initialization, otherwise it is unstable.

Preliminaries

Describing sensitivity to initial conditions is usually done with Lyapunov exponent [A. 85, M.T93, Hil00]. As explained in 3.4, this quantity evaluates the divergence rate λ between two trajectories initially as close as possible 17 . As an example, we consider two very close likelihoods vectors: 6.10) where { i } i are random disturbances between 0 and a maximum value we set to 10 -618 . For each set of N disturbances, likelihoods vectors are injected in a given decoder, either BP or GBP. Beliefs are computed at each iteration k to compute the Euclidean distance d k between trajectories:

L A = [l A,1 (x 1 , y 1 ) . . . l A,N (x N , y N )], L B = [l B,1 (x 1 , y 1 ) . . . l B,N (x N , y N )]. We practically define: ∀X i ∈ X, ∀x i l B,i (x i , y i ) = l A,i (x i , y i ) + i (
d k = N ∑ i=1 b (k) A,i (0, y i ) -b (k) B,i (0, y i ) 2 (6.11)
We observe the iterative evolution of log d k , known as the log-distance. The use of the logarithm function appears legitimate as theoretically, the largest divergence, or separation, found between two trajectories is an exponential function of k. The logarithm allows us to cancel the exponential and to extract the divergence rate. We present in Fig. 6.14 evolutions of the log-distance of BP for SNR values belonging to SNR intervals previously extracted. In figures Fig. 6.14(a) and Fig. 6.14(b), the distance converges to a very small non zero value, meaning that both trajectories have reached stationary and very close points in the state space. We evaluate d K around 7.10 -6 for both SNR, which is around the same magnitude as the distance between the likelihoods. Both trajectories then do not substantially get away one from the other. We note a transient phase for SNR = 2.17dB that lasts about 600 iterations where trajectories seem to pseudo-oscillate one around the other before collapsing. This behavior is an announcement of instability, given that as the SNR is increased, this transient phase lasts more and more.

When increasing the SNR to 2.35dB, see Fig. 6.14(c), the log-distance is getting larger in average, i.e. trajectories move away from each other. Furthermore, oscillations do not stop anymore, contrary to the previous SNR values, meaning that trajectories are also turning around. As the SNR is increased theses oscillations change to strong divergence between trajectories, e.g. the log-distance reaches a value of 0.19 for SNR = 2.90dB, represented in Fig. 6.14(d), meaning that the divergence has been multiplied by 2.10 4 between both, that is a considerable change. In addition, we observe an exponential divergence rate up to 200 iterations, i.e. trajectories have reached the maximum separation, that is a sign of very high sensitivity to the initial conditions. After 200 iterations, d k reaches a stationary state, called plateau, of value d K = 0.19, a substantial value. This plateau only appears because we observe the logdistance, as the log-function compresses values of the true distance. As an example, we display in Fig. 6.15 the distance at SNR = 2.90dB for the BP. Values of d k after the exponential divergence range in [0.1; 0.23], with chaotic oscillation amplitudes. It is then erroneous to conclude that the apparent constant log-distance corresponds to a constant d k . We specify that the maximum value of d k does not exactly reflect the size of the chaotic attractor, we will tackle this point in section 6.8.

Increasing the SNR beyond the next critical value around 3dB makes trajectories getting closer and closer. We still observe, in Fig. 6.14(e), an exponential divergence during first iterations k ≤ 200, but the distance then exponentially decreases, i.e. trajectories that first seemed to move away eventually slowly merge. This phenomenon is confirmed, in Fig. 6.14(f), with SNR= 3.12dB, where trajectories shun one another during 45 iterations and from k = 45, they completely fuse together.

Lyapunov exponent

Lyapunov exponent is a good candidate to accurately reflect these observations. According to 3.4, Lyapunov exponent λ is the divergence rate between two initially close trajectories:

d k = d 0 e λk (6.12)
Sign and magnitude of λ indicate the divergence or convergence nature of a given dynamical systems. Exponentially divergent trajectories are reflected by λ ≥ 0, whereas merging trajectories correspond to λ < 0. We recall 3.4 that the practical method for computing λ is to extract the log-distance slope given by: log d k = log d 0 + λk (6.13)

Usually, computing λ is well-suited for chaotic attractor as they are responsible for exponential divergence. More precisely, the shape of log d k help deduce whether a dynamical system is attracted into a chaotic attractor or not. As a matter of fact, the value at last iteration K strongly conditions the slope value, and then the nature of the corresponding attractor. As an example, for an SNR of 3.05dB, we see, in Fig. 6.14(e), that taking K = 200 provides a higher λ value than taking K = 1000 because the shape of log d k is changing.

Furthermore, for SNR= 3.12dB, selecting K = 45 involves a positive λ whereas K ≥ 200 involves a negative λ. Accordingly, we strongly recommend to define Lyapunov exponent as the linear regression slope not on the transient part of log d k but on the whole values of k. This way, we combine emerging divergence and possible steady state with merging evolution, we do not omit any information. We specify that the last iteration K is restricted to one hundred in order to keep consistency with previous results on bifurcation diagram, reduced trajectory and MSB. The linear regression provides the following equation of Lyapunov exponent:

λ = K ∑ k=1 k log d k - 1 K K ∑ k=1 k K ∑ k=1 log d k K ∑ k=1 k 2 - 1 K K ∑ k=1 k 2 (6.14)
We show in Fig. 6.16 the λ evolution19 as a function of the SNR. We observe that the SNR intervals indicated by the bifurcation diagrams are still valid, especially for BP:

• first interval: no change in distance, d k ≈ d 0 ,

• second interval: distance increases to a stationary value, trajectories turn around inside a limit cycle,

• third interval: distance increases meaning that trajectories follow very different evolutions, trajectories are trapped into a chaotic attractor,

• fourth and fifth intervals: λ reaches infinite negative value, trajectories merge to a fixed-point.

However, non-trivial behaviors of GBPs appear earlier than BP. The values of λ significantly increase near 2.12dB. Furthermore, while Lyapunov exponent of BP is stabilized to λ = 3.10 -3 in the second interval, GBP-P and GBP-IP exhibit a stronger sensitivity to initial conditions as λ ≈ 0.8.

In contrast, GBP-C presents values of λ comparable to those of BP, that indicates that their sensitivity are quite equivalent. We then observe a common decreasing behavior in the third interval for all decoders, we can then conclude that all decoders are subject to a similar chaos intensity. We specify that it does not imply that chaotic attractors are identical, but that their sensitivity to a slight change in the likelihoods are equivalent. By the use of the Lyapunov exponent, we are able to quantify the behaviors nature and their intensity, that helps confront decoders. We now deal with the size of the attractors so as to relate them to decoders performance. § 6.8 Hyperspheres and attractors width

Static hypersphere

The nature of an attractor, e.g. fixed point, limit cycle or chaos, is not enough to describe decoders and to relate with their decoding performance. It appears that widths of chaotic attractors may be very unpredictable and that they largely influence algorithms performance: as the size of a chaotic attractor is increased, beliefs exhibit more and more fluctuations making decoding more and more troublesome.

The size of an attractor is here defined as a measurement of its hypervolume in the state space. As we previously saw, nothing ensures that shapes of chaotic attractors are usual geometric forms (spheres, ellipsoids, hypercubes, etc.), therefore computing the value of the hypervolume of a given chaotic attractor is unfortunately hopeless because it is completely dependent on its shape. To circumvent this problem, we establish a procedure, that we called the hyperspheres method: for a given attractor, it consists in approximating its hypervolume by the hypersphere circumscribed to any of its trajectories. Practically, we only compute the radius R of the hypersphere given that the corresponding hypervolume is proportional to R N . R is then a measure of attractors width. Figure 6.17: Lorenz chaotic attractor: hypersphere from d K in blue, from our own method in red Practically, we compute R first by extracting the mean point of the attractor, then by searching for the attractor point that is the farthest one. The Euclidean distance between them is the radius of the hypersphere. One would point out here that any attractor width should be provided by the maximum of d k , the distance between initially nearby trajectories. This is quite realistic as the radius R resulting from the hyperspheres method is not far from the maximum of d k . As an example, we consider the Lorenz map in its chaotic behavior, see chapter 3.2.2 for details. We recall the profile of the chaotic attractor in Fig. 6.17. The red hypersphere stems from our own computation method whereas the blue one is computed by means of distance d k . We observe that the attractor is enclosed inside the red one without any ambiguity, but not in the blue one, i.e. it is less accurate than our method. To reach same results for both methods, it would be required to compute an average radius over a numerous number of trajectories couples, but it represents a too large computation time, compared with the hyperspheres method.

We display experimental results on BP and GBP in Fig. 6.18, where are exhibited radii of associated circumscribed hyperspheres as functions of the SNR. We observe that BP falls into larger and larger attractors as the SNR is increased, from 2.18dB to 2.80dB, where the width reaches its maximum value R BP = 3.30. As a matter of fact, results show that chaotic attractors of BP, and of GBPs too, are the largest ones. In the last SNR values, R BP decreases slowly up to an abrupt peak at 3.00dB. This may be due to the fact that chaos remains only a transient behavior at this SNR value, the steady state being the right decoding state where R BP ≈ 0.

In other words, convergence jumps, see Fig. 6.13, make inaccurate the method when attractors are not centered around the computed mean point but farther. We also observe this phenomenon with GBPs around the same SNR value. Contrary to GBP-C which static radius follows the same evolution as for BP, we see that GBP-P and GBP-IP present a significant and early increase in the attractors width, at 2.05dB and 2.10dB, respectively. While BP and GBP-C are still enclosed in small-sized attractor, attractors of GBP-P and GBP-IP get substantially larger up to their average steady value R = 3.25. By increasing the SNR, both decoders fall into chaotic attractors which widths are not stable, contrary to BP where fluctuations of R BP are less important.

Finally, the static radius reveals that damped GBPs have their trajectory enclosed in quite constant-sized attractors, their widths are a little larger than the maximum width of BP attractors.

Local hypersphere

The hyperspheres method allows us to find an estimate of any attractor radius R. We observe when SNR reaches value of convergence jumps appearance that R was not very reliable, because transient disturbs the measurement. On one hand, transient phase, especially for chaotic attractors, may last about 200 iterations as we can see in Fig. 6.14(d), which is a non negligible value. On the other hand, we previously mentioned that the number of iterations should not exceed a hundred to preserve a consistency with realistic situations. Therefore, it would be very useful to determine how to modify the hyperspheres method to fit with realistic situations that may not reach steady states.

The idea is to transform static radius R to a dynamical radius R k : rather than considering the whole trajectory, we consider a temporal window, of a finite length W < K, sliding step by step along the trajectory, each step k giving the radius R k of the hypersphere circumscribed to the partial trajectory. The hyperspheres method is then transformed such that it results in a sequence of

local radii {R 1+ W 2 , . . . , R K-W 2 }. Practically, for each point T k of the trajectory, k ∈ {1 + W 2 , . . . , K -W 2 }, we compute the trajectory mean point B k in the temporal window W k = ]k -W 2 , k + W 2 ].
After that, we search for the farthest point to B k inside W k . R k is then computed as the Euclidean distance between both points, see representation in Fig. 6.19.

}

The dynamical radius R k allows us to observe iterative evolutions of attractors using temporal windows. The term attractor has slightly changed as it refers here to the set of points of a given trajectory enclosed in W k , i.e. it is only a partial attractor. The dynamical hyperspheres method then use local average behaviors to observe iterative evolutions of dynamical systems.

For any iteration k, a large value of R k means that the given decoder has significant changes, whereas a small value of R k indicates that the decoder is slightly changing. In other words, local, or dynamical, hyperspheres provide the decoder scattering.

The window length W wields influence on results of the method. On one hand, selecting a too large value of W could make the method find a non varying radius along iterations as it would be too close to the static radius, i.e. the dynamical local hyperspheres would be similar to the static hypersphere circumscribed to the whole trajectory. On the other hand, opting for a too small length W implies that local hyperspheres are completely determined by trajectory points themselves, i.e. no information can be enlightened on the attractor width. In our experiments, we chose a window of length W = 30. The reason for this choice is that this value is around the minimum practical number of iterations used to run any decoder 20 . We present in Fig. 6.20 the R k evolution for BP and GBPs when trapped in attractors previously mentioned 21 . When the SNR equals 2.00dB, decoders are attracted by fixed-points, as depicted in Fig. 6.10. We previously saw that corresponding reduced trajectories drew spirals in reduced spaces, justifying R k as hyperbolic functions of k. We now consider that SNR = 2.30dB. We note that BP radius is quite constant compared with GBPs ones. At this SNR value, BP is trapped into a limit cycle, that explains the relative stability of R k . Concerning the GBP algorithms, we observe that the R k shape of GBP-C behave in the same way that BP one with an increasing dilatation as k increases, meaning that the associated limit cycle is getting larger and larger. A noteworthy point is that the GBP-P radius diminishes for k ≥ 48, i.e. the limit cycle attractor retracts. In contrast, the GBP-IP radius does not present any decreasing values before k = 72, meaning that the damping factor does not help focus on a single point of the state space.

As the SNR is set to 2.90dB, decoders are trapped into chaotic attractors. These attractors are larger than the previous limit cycles, as R k values are increased about a unit. We observe that the BP attractor slowly collapses as k is increased, that could be approximated by an affine function, which makes emerge that inside this chaotic attractor, fluctuations are linearly decreasing. While k ≤ 60, we can distinguish that GBP-C and GBP-IP follow a BP-like evolution, i.e. R k slowly decreases in average too. In contrast, GBP-P evolves pseudo-oscillating in its chaotic attractor with R k values close to GBP-C ones. As k > 60, each GBP changes its behavior, due to damping factors. Radii of GBP-C and GBP-IP grow up to R k > 3 whereas GBP-P radius decreases to R k ≤ 0.5, which is even lower than BP radius. Finally, for a SNR of 3.09dB, all GBPs converge to smaller attractors than BP one. More precisely, GBP-C and GBP-IP both converge to a fixed-point (R k = 0) and GBP-P slowly converges to a smaller and smaller attractor, whereas BP is still stuck to an attractor of radius R k ≈ 2. This piece of information could not be revealed by the static radius, that makes the use of the local radius relevant.

These observations tend to conclude that damping factors wield substantial influence on GBPs attractors, especially when these attractors are nontrivial. We see that to select either a damping rule or another one deeply conditions GBP performance. However, when approaching the major bound of the chaotic interval, the GBP algorithms tend to converge to fixed-points contrary to BP that is still stuck in an attractor of a non-negligible size. In other words, GBP diverge more but shorter than BP.

Thanks to the hyperspheres method, we are able to offer measurements of attractors hypervolumes for all decoders, which particularly helps contextualize chaos in decoding performance. The remaining question now is where these attractors take place in their associated state space, or more simply, where trajectories evolve in the state space. § 6.9 Lengths ratio to determine the relative position in the state space

By means of previous descriptors, we now have information about attractors: we know their nature, their size and the SNR values for which they appear.

The last description we need is about relative positions of these attractors to the perfect decoding state in the state space, that will be offered by the last quantifier, called the complementary length.

Complementary length

The length of any beliefs vector, in the state space, is defined as its Euclidean distance to all-zero point:

L(k) = N ∑ i=1 b (k) i (0)
2 (6.15) which is very close to MSB. To keep connections with concept of length and positions, we prefer to keep using L(k) instead of E(k). We now define the complementary length as the Euclidean distance to all-one point:

L(k) = N ∑ i=1 1 -b (k) i (0) 2 (6.16)
Experimentally, we always consider that the emitted codeword is the zerocodeword, therefore, all-one point in the state space corresponds to perfect decoding 22 . The vector length is then the distance to completely faulty decoding, or wrong decoding. As a matter of fact, only one of these two quantities does not provide decoding performance of a given decoder. We display a toy example in Fig. Coordinates of points A,B and C are given in table 6.3. We include in this table lengths and complementary lengths of vectors -→ OA, -→ OB, -→ OC. Points A, B, C are equidistant to all-zero point O, as L A = L B = L C . However, they are not equidistant to all-one point M as LA < LC < LB , i.e. point B is the closest to point M. This kind of information is not given by the only knowledge of length L, as it does not take into account any other point that M. Using also L helps faithfully represent the relative position of any decoder in its state space.

Lengths ratio

We now apply concept of complementary length to extract positions of decoders inside their own state space and more accurately, their vicinity with the optimal point 23 . Instead of using length and complementary length, we introduce a combination thereof to obtain a unique function that reflects the relative decoders positions to the correct decoding state 24 . This combination is called the lengths ratio and is given by: We link typical values of α k with decoding situations:

α k = L(k) L(k) (6.
• α k = 0: perfect decoding, L(k) = 0,

• α k = +∞: completely faulty decoding, L(k) = 0,

• α k = 1: indecisive decoding, entropy is maximum on the beliefs.

Thus, the smaller α k is, the more relevant the decoding is. In Fig. 6.23 are displayed iterative evolutions of α k for BP and GBP. We observe for SNR= 2.00dB that all decoders are exactly in the same relative position as α k = 0.4. This does not mean that beliefs are common to all of them, but that trajectory points are equidistant to all-zero and all-one points.

Once the SNR reaches 2.30dB, we know that decoders oscillate, especially BP and GBP-C25 . Only small fluctuations are visible around the average value 0.4, meaning that oscillations are not very large compared with the size of the state space26 . We previously mentioned, by means of static radius of hyperspheres, that attractors are getting larger as the SNR is increased. This is why at 2.90dB, significant differences between decoders behaviors are visible. In addition, attractors are chaotic and of non-negligible widths.

Lengths ratio of any decoder is then subject to important variations. GBP-IP lengths ratio has the great advantage of converging to the perfect decoding state, contrary to other ones. However, by increasing again the SNR, e.g. 3.09dB, we observe that GBP-C also presents such a behavior. Surprisingly, GBP-P still does not converge to the right state even though its lengths ratio has been diminished compared with that of the previous SNR. In addition, this lengths ratio is lower than the BP one, indicating that BP is even more remote to the right state than GBP-P. Not displayed experiments have showed that increasing again the SNR makes lengths ratios of all algorithms converge to zero.

The noteworthy point is that the GBP algorithm, especially damped with the inverse parabolic law, presents a better proximity to the perfect decoding state than other decoders. This fact is not particular to the error event EE we selected. Other experiments demonstrated that the lengths ratio of GBP is generally smaller than the BP one. In other words, for any initial condition, GBP approaches the right state more accurately than BP.

Practically, performance of any decoder are measured once either messages have converged, or the estimate word is valid according to the LDPC code, or the last iteration has been reached. We then picked out the value of the steady lengths ratio α K for different SNR values, K being the iteration where the given decoder stopped 27 . We average values of α K over a thousand error events close to the original error event EE, at the Euclidean sense. Results are represented in Fig. 6.24. For SNR values lower than 2.00dB, lengths ratios of all decoders merge in a single line, as well as when SNR values greater than 3.1dB. As previously explained, the α K behavior may fluctuate more and more as the SNR is increased up to 2.30dB. As a summary of previous observations, we observe that from this SNR value, GBP-IP tends to approach the right decoding state significantly earlier than other decoders, that is a manifestation of the convergence jumps. In addition, we see that just after 3.00dB, GBP-P and GBP-IP converge to all-one point, whereas GBP-C and BP has to wait a little bit longer. Furthermore, BP exhibits a rebound once all other algorithms have converged, which is a sign of worse decoding. § 6.10 Conclusion

The study of dynamical behaviors of iterative decoders brought out properties cannot be deduced from the only knowledge of the average BER. We saw that BP and GBP were comparable in terms of BER when the SNR ranges in [2.00; 3.10]dB. But when focusing on non-trivial error events, we observed different evolutions. Bifurcation diagram provided critical SNR values for which message-passing algorithms encountered particular attractors that prevent them from easily converging. We observed that for all decoders, critical SNR values and even SNR intervals are similar, which makes them very interesting challengers. Lyapunov exponents and reduced trajectories indicated that algorithms are trapped by chaotic attractors in the same SNR interval, and it was clear that their sensitivity to initial conditions in this interval was similar, even though it is not the case for lower SNR values, partly due to damping factors. We introduced few original estimators to analyze widths and relative positions of attractors and trajectories in state spaces. We motivated this by trying to relate dynamical behaviors to decoding performance, that is not an easy task. In addition, the usual estimators from the literature on dynamical systems theory cannot reflect the whole influence of dynamical properties on BP and GBP efficiency. By means of the hyperspheres methods, we first extracted attractors size: it appears that damping rules make GBP evolve in attractors of a non negligible but stationary size, contrary to BP that is enclosed in larger and larger attractors as the SNR is increased, up to the critical SNR values where all algorithms become chaotic. Related to reduced trajectories, chaotic attractors of GBP are more stretched in the direction of the right state whereas BP has a quite uniform shape in its state space. By the use of the lengths ratio, we exhibited the vicinity of decoders to the right state, more accurately than what has been done with MSB. We noticed once more that the chaos disappear earlier for GBP than BP, thanks to decreasing damping rules.

Conclusion

We demonstrated in first chapter that fixed-points of BP are stationary point of Bethe approximation. This property allowed us to extend connections between decoding algorithms and statistical physics. Both methods are subject to a significant suboptimality due to the non-convexity of variational free energies, in deep correlation with the loop-like topology of most Tanner graphs.

We demonstrated in second chapter that Bethe approximation could extend to Region-Based Approximation (RBA), which provided the GBP algorithm, a message-passing algorithm running on edges of a region-graph. Goal of the region-graph construction is to absorb harmful topological structures of any Tanner graph, such that resulting RBA is less damaged than Bethe approximation. We made use of this principle to introduce an original region-graph construction based on trapping sets of Tanner code, that are responsible for the degradation of BP performance, especially for high SNR values. By means of an original method, we brought out relevant results of the GBP algorithm for error events based on trapping sets. We observed that this decoder could surpass BP in these situations, making our own construction well-suited.

In third chapter, we introduced a toolbox to examine dynamics of any iterative decoder. This helps understand non-trivial behaviors of algorithms, e.g. the lack of convergence, oscillations or even chaos. By the use of bifurcation diagrams and Lyapunov exponent, we extracted typical attractors that trapped decoders, and SNR values they appear at. We introduced original descriptors of dynamics, that we called hyperspheres and lengths ratio, to discuss size and position of these attractors, and they allowed us to determine their relations with decoding performance of BP and GBP. It appeared that the triplets clustering presented very relevant properties as associated attractors are closer to the right decoding state than BP ones, though Lyapunov exponent demonstrated that BP was a little more stable. In addition, damping rules, that we introduced to overcome the lack of convergence of the pure GBP algorithm, provided that attractors always shrink as the iteration increases, contrary to BP.

Statistically, BP and GBP algorithms are quite comparable. Without any specific construction of the region-graph, results of the GBP algorithm are strictly equal to the results of the BP. This is generally observed as the systematic construction for LDPC codes provides a region-graph which topology is identical to the Tanner graph one. With our dedicated construction, we show that it is possible to make a region-graph such that GBP surpasses BP.

Conclusion

In this thesis, we addressed inference problems on loopy factor graphs. We presented the Belief Propagation algorithm (BP) used as a message-passing algorithm to almost solve inference. We exhibited its deep relations with statistical physics, especially with Bethe approximation, a method that is aimed at minimizing the variational free energy to obtain the log-partition function. By the use of examples and experiments, we showed that loop-like topology of most factor graphs, and more precisely Tanner graphs for LDPC codes, prevents BP from being optimal.

Introduction of the region-based approximation, and the associated Generalized Belief Propagation algorithm (GBP), helped circumvent this problem by mapping the factor graph to the region-graph. Provided that the region-graph is constructed according to the factor graph topology, the GBP algorithm results in more efficient performance than BP. In the case of LDPC codes, we used trapping sets to construct a specific region-graph dedicated to the Tanner code. By means of damping factors and new principles, we developed a method for high SNR values that demonstrated that the GBP algorithm was able to outperform the BP algorithm.

The comparison of both decoders for middle SNR values turned out non conclusive as they exhibit equivalent performance. From the dynamical systems theory point of view, we introduced a relevant investigation that brought out decoders typical behaviors. Bifurcation diagrams offered that the SNR values for which they do not trivially evolve were common. Furthermore, both are subject to chaos for same SNR values. The use of Lyapunov exponent and reduced trajectories revealed that damping factors of the GBP algorithm wield much influence on the decoder stability at a low memory cost. In addition, these functions allowed us to observe that the chaos intensity was similar to both, but that GBP always converged earlier than BP to the right decoding state, especially with a decreasing damping factor. We introduced the static hyperspheres method to approximate the attractors width, especially for chaotic attractors. We observed that the widths of the GBP attractors were quite stationary compared with the BP ones that grew as the SNR is increased. Local hyperspheres helped us conclude that along the iterations, the GBP algorithm is more convergent than BP, due to the damping factor particularly. The introduction of the complementary length and the lengths ratio provided more relevant information than the mean square beliefs as it faithfully reflects the relative position of any decoder to the right decoding state.

We saw that the GBP algorithm presented better length ratios when approaching the major bound of the chaotic interval.

Finally, the specific region-graph construction with a triplets clustering turned out powerful according to the number of errors decoded by the GBP algorithm compared with the BP's. Furthermore, the investigation on dynamical aspects revealed that the GBP behavior provided with damping factors offered a better convergence than the BP algorithm.

Future works

From this work have emerged a certain number of issues that constitute an interesting basis for future works. Points we are about to tell belong to the different research areas that were tackled in this thesis.

• First of all, the clustering in triplets has to be adapted to other LDPC codes of the same family as the Tanner code's, to acquire a larger overview of GBP decoding performance. In addition, it will be necessary to introduce constructions based on other kinds of trapping sets. The investigation here remains on TS (5,3) but nothing prevents from considering other structures, provided that they are small enough to preserve a reasonable computation complexity.

• The computation complexity of the GBP algorithm is its weak point. Therefore, we began to study the decimation of the region-graph. The method that we started to investigate consists in switching off a fixed percentage of edges -to some extent, messages -in the region-graph that we randomly select without any discrimination. After a given number of iterations, we switch them on and we switch off other edges similarly selected. The two parameters of this method are the percentage of switched off edges and the number of iteration before changing in the edges. Preliminary simulations showed that with 30% of switched off edges during about 10 iterations, the GBP algorithm does not present significant difference in its performance in terms of the BER. By increasing the percentage, we induce a too large modification that damage GBP results. Therefore, it would be relevant to consider a more specific rule of edges selection instead of a uniform one.

• Recently was introduced an alternative to the BP algorithm that consists in weighting messages [H. 12],[H. 11],[T.G08]. It reminds us of the damping factor of the GBP algorithm, even though the mathematical details are different and the goal is not to ensure convergence but to improve results. Furthermore, experiments show that regular BP, i.e. BP provided with a unitary weight, is outperformed by non unitary weighted BP. It appeared that this weighting should apply to other iterative decoders, such as GBP, FAID, etc.

• The lack of convergence and the suboptimality of the decoding algorithms are due to the loop-like topology of the Tanner graphs, as we try to minimize non-convex variational free energies. The CCCP proposed in [A.L03] is a powerful tool to treat this problem, even though the double-loops imply a significant increase in the computation time. As a matter of fact, CCCP has been derived for Bethe and region-based variational free energies. Experimenting CCCP on LDPC codes, especially with well-suited region-graphs, e.g. the Tanner code provided with the triplets construction, to compare with the original GBP algorithm would be an interesting work.

• In [T. 03a], authors deal with convexity of the Bethe variational free energy. Studying in the same way the region-based free energy to extract sufficient conditions for its convexity, or even to convexify it, would solve many problems.

• One of the drawback of decoders dynamics is that we do not know basins of attraction of attractors. It would be very useful to know where chaotic attractors are in the state space, that would help to control the BP evolution. Independently from the work presented in this thesis, we have demonstrated that it was possible to write a "reverse BP", i.e. a message-passing algorithm that goes back in time from a given point. It would help associate to any state its initial conditions. "Downdate" equations unfortunately do not result in a unique previous state but in numerous states. In addition, what is true for regular BP is also true for its reverse version, i.e. the sensitivity to initial conditions. Given a point in the state space or a very close neighbor, we obtain by reverse BP very different point when dealing with chaotic attractors. Therefore, we did not go further in this method but it needs a deep investigation to reveal relevant results that would help control BP dynamics, and even GBP dynamics.

List of Figures

  Figure 1.1: Channel

Figure 1 . 2 :

 12 Figure 1.2: Two channel models

Figure 1 . 3 :

 13 Figure 1.3: Bit-error rate of an uncoded transmission

Figure 1 . 4 :

 14 Figure 1.4: Model with an error-correcting code

Figure 1 . 5 :

 15 Figure 1.5: Maximum likelihood decoding on the Hamming code

Figure 1 . 6 :Figure 1 . 7 :

 1617 Figure 1.6: 2-level Tanner graph of the Hamming code

Figure 1 . 9 :

 19 Figure 1.9: Messages in a factor graph

Figure 1 . 10 :

 110 Figure 1.10: Tanner code decoding (N = 155)

Figure 1 . 11 :

 111 Figure 1.11: CCCP decoding of the Tanner code

Figure 1 . 12 :

 112 Figure 1.12: Inner loop of the CCCP decoding of the Tanner code

Figure 1 . 13 :

 113 Figure 1.13: Tanner code (N = 155) decoding by BP on AWGNC

  2.1.

Figure 2 . 1 :

 21 Figure 2.1: Spin glass: the dashed lines are the anti-ferromagnetic interactions and the solid lines are the ferromagnetic interactions, the arrows are the spins (up or down)

Figure 2 . 2 :

 22 Figure 2.2: Spin glass: dashed lines are anti-ferromagnetic interactions and solid lines are ferromagnetic interactions; arrows are spins (up or down). The red line is an interaction modified by an external noise and red spins are spins modified by interactions.

Figure 2 . 3 :

 23 Figure 2.3: Ising chain with N = 3 spins embedded in an external magnetic field

  Figure 2.7: p-range model

  Figure 3.1: Logistic map x k+1 = 4x k (1x k )

Figure 3 . 2 :

 32 Figure 3.2: Logistic map

Figure 3 . 3 :

 33 Figure 3.3: Lorenz attractor As mentioned previously, a detailed study would require many values of the parameters σ, r, b. In this context, we chose values that provide non trivial phenomena: σ = 10, (3.5) b = 8 3 , (3.6) r = 28. (3.7)

Figure 3 . 4 :

 34 Figure 3.4: Lorenz map

Figure 3 . 5 :

 35 Figure 3.5: Logistic map -Period doubling

Figure 3 . 6 :

 36 Figure 3.6: Bifurcation diagram of the logistic map

Figure 3 . 8 :

 38 Figure 3.8: Tinkerbell attractor with a = 0.5

Figure 3 . 9 :

 39 Figure 3.9: Tinkerbell attractor with a = 0.1

Figure 3 . 10 :

 310 Figure 3.10: Tinkerbell attractor with a = 0.9

Figure 3 .

 3 Figure 3.11: Lyapunov exponent estimates

Figure 3 . 12 :

 312 Figure 3.12: Time evolution of the separation

Figure 3 .

 3 Figure 3.13: Fake chaos, K = 300

Figure 3 .

 3 Figure 3.14: Solid line is the alternative Lyapunov exponent, dashed line is the Lyapunov exponent from the differential method

Figure 3 .

 3 Figure 3.15: Hénon map state space

Figure 3 . 16 :

 316 Figure 3.16: Lorenz map state space

Figure 4 . 2 :

 42 Figure 4.2: Tanner tree

Figure 4 . 3 :

 43 Figure 4.3: Square Tanner lattice

Figure 4 . 4 :

 44 Figure 4.4: Square lattice of length N = 9

Figure 4 . 5 :

 45 Figure 4.5: Hamming code

Figure 4 . 6 :

 46 Figure 4.6: Two formulations of the normalization constraint. (a) We only scan X to extract point A, the argument of the minimum. (b) Two points that live on the same line (B, B ), (C, C ) or (D, D ) will have the same mapping by f • ||•|| then the argument of the minimum if not a unique point but a line.

Figure 4 . 7 :

 47 Figure 4.7: Wrong formulation of the normalization constraint. The projection of M on X is different from A.

( 4 . 59 )

 459 Beliefs are solutions to what is called primal problem [D.P01]. At this point, Lagrange multipliers are not valued, their assignments is the dual problem.

Figure 4 . 8 :

 48 Figure 4.8: Tree-like Tanner graph

Figure 4 . 9 :

 49 Figure 4.9: Non tree-like Tanner graph

  It is then impossible to reach canonical equation (4.68), making BP a suboptimal algorithm.

Figure 4 .

 4 Figure 4.10: Local Tanner graph for messages computations

Figure 4 . 11 :

 411 Figure 4.11: Tanner graph of a loop

Figure 4 . 12 :

 412 Figure 4.12: Trapping-set of a = 5 variable nodes and b = 3 unverified check nodes

Figure 4 .

 4 Figure 4.13: BP error floor when decoding the Tanner code covered by TS(5, 3).

Figure 5 . 1 :

 51 Figure 5.1: Square lattice of length N = 9

Figure 5 . 2 :Figure 5 . 3 :

 5253 Figure 5.2: A Bayesian network as a graphical representation of (5.2)

Figure 5 . 4 :

 54 Figure 5.4: Average errors of MF, Bethe and region-based approximations on square lattice N = 9

Figure 5 . 5 :

 55 Figure 5.5: Two clusterings for two factorizations

X 5 (Figure 5 . 6 :

 556 Figure 5.6: Two region-graphs for two factorizations

  observations Y r = {Y i |X i ∈ X r } and Bayes rule, we obtain the a posteriori distribution: b * r (x r , y r ) = f r (x r )l r (x r , y r )e -1 e -1 cr γ r ∏ s∈E r e -1 cr λ rs (x s ) ∏ q∈P r e 1 cr λ qr (x r )

  34)×b r (x r ) gives:c r b r (x r ) + ∑ q∈A r \p∪A p c q b r (x r ) = 0 (5.35)

Figure 5 . 7 :

 57 Figure 5.7: Explicit square lattice of length N = 9

Figure 5 . 8 :

 58 Figure 5.8: Region-graph of the explicit square lattice

Figure 5 . 9 :

 59 Figure 5.9: Small region-graph Part of the Lagrange function dedicated to local marginalization constraints (C2) is made with following equations:

Figure 5 . 10 :

 510 Figure 5.10: Reduced region-graph

Figure 5 . 11 :

 511 Figure 5.11: Uniformly reweighted BP (ρ = 1.00 is BP)

Figure 5 .

 5 Figure 5.12: BER of the GBP with changing w 0

Figure 5 . 13 :

 513 Figure 5.13: Number of iterations of the GBP with changing w 0

Figure 5 . 14 :

 514 Figure 5.14: Performance of the GBP for three weights

Figure 5 .Figure 5 . 16 :

 5516 Figure 5.15: BER of the GBP with changing w 0

Figure 5 . 17 :

 517 Figure 5.17: Damping factor with initial value w 0 = 0.85

Figure 5 .

 5 Figure 5.18: BER of the GBP algorithm on Tanner code with decreasing damping factors

Figure 5 .

 5 Figure 5.19: BER of the GBP algorithm on the Hamming code

Figure 5 . 20 :

 520 Figure 5.20: Parallel version performance

Figure 5 .

 5 Figure 5.21: A systematic region-graph for the Hamming code

Figure 5 . 22 :

 522 Figure 5.22: Suited region-graph for the square lattice code N = 9

Figure 5 .

 5 Figure 5.23: Bit-error rate of the BP and the GBP algorithms on the square lattice code N = 36

Figure 5 . 24 :

 524 Figure 5.24: Trapping set TS(5, 3)

Figure 5 . 25 :

 525 Figure 5.25: Split of the TS(5, 3) -Quadruplets construction

Figure 5 . 26 :

 526 Figure 5.26: Split of the TS(5, 3) -Triplets construction

Figure 5 . 28 :

 528 Figure 5.28: Third quadruplet as a cluster

Figure 5 . 32 :

 532 Figure5.32: BER on the Tanner code for EE completely dependent on TS(5, 3)

Figure 5 . 33 :

 533 Figure 5.33: Modification of moments -Dashed lines are original SNR and mean, solid lines are their practical estimates according to EE

  SNR i = 18.0dB

Figure 5 .

 5 Figure 5.34: BER of the BP and the GBP algorithms on a EE with an AWGNC

Figure 5 .

 5 Figure 5.35: Average BER of the BP and the GBP-P algorithms

Figure 6 . 2 :

 62 Figure 6.2: Likelihoods space at SNR = 2.3dB on the Tanner code for quenched likelihoods except on variable nodes X i = 100, X j = 102 -BER of the BP

Figure 6 . 4 :

 64 Figure 6.4: GBP MSB on EE

Figure 6 . 6 :

 66 Figure 6.6: Bifurcation diagram of the BP on the Mac Kay code

Figure 6 . 7 :Figure 6 . 8 :

 6768 Figure 6.7: Influence of K on the bifurcation diagram of the BP

Figure 6 . 9 :

 69 Figure 6.9: Chaos in the logistic map x k+1 = ax k (1x k ), for a ≥ 3.44

Figure 6 . 11 :

 611 Figure 6.11: Reduced trajectory with SNR = 2.30dB

Figure 6 . 12 :Figure 6 . 13 :

 612613 Figure 6.12: Reduced trajectory with SNR = 2.90dB

Figure 6 . 14 :

 614 Figure 6.14: Evolution of the BP log-distance

Figure 6 . 15 :

 615 Figure 6.15: Evolution of the BP distance at 2.90dB

Figure 6 . 16 :

 616 Figure 6.16: Lyapunov exponent on the Tanner code

Figure 6 . 18 :

 618 Figure 6.18: Radii evolution of BP and GBP

Figure 6 . 19 :

 619 Figure 6.19: Hypersphere of radius R k centered on mean point B k of the partial trajectory {T k-W 2 , . . . , T k-W 2

Figure 6 . 20 :

 620 Figure 6.20: Local radii R k evolution for BP and GBP

Figure 6 . 21 :

 621 Figure 6.21: Length and complementary length

  17

)

  This allows us to judge about the accuracy of BP and GBP. For any couple of beliefs vectors of identical lengths, α k brings out the closest one to the right decoding state which is also the farthest one to the wrong decoding state. In the previous example, lengths ratios associated to points A, B, C are: α A = 1.0279, α B = 1.1692 and α C = 1.0850. As their lengths are equal, we consistently find that α A < α C < α B . Two points of identical lengths ratios are comparable in the state space as their decoding performance are similar. We display in Fig.6.22 the surface covered by α k in the plane generated by all possible values of L(k) and L(k).

Figure 6 . 22 :

 622 Figure 6.22: Lengths ratio

Figure 6 . 23 :

 623 Figure 6.23: Iterative evolution of the lengths ratio

Figure 6 . 24 :

 624 Figure 6.24: Steady lengths ratio along the SNR for BP and GBP

1. 1

 1 Channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Two channel models . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Bit-error rate of an uncoded transmission . . . . . . . . . . . . . 1.4 Model with an error-correcting code . . . . . . . . . . . . . . . . 1.5 Maximum likelihood decoding on the Hamming code . . . . . 1.6 2-level Tanner graph of the Hamming code . . . . . . . . . . . . 1.7 Tanner code (N = 155) hard decoding on a BSC of crossover probability p . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Message m ij in a Bayesian network . . . . . . . . . . . . . . . . . 1.9 Messages in a factor graph . . . . . . . . . . . . . . . . . . . . . . 1.10 Tanner code decoding (N = 155) . . . . . . . . . . . . . . . . . . 1.11 CCCP decoding of the Tanner code . . . . . . . . . . . . . . . . 1.12 Inner loop of the CCCP decoding of the Tanner code . . . . . . 1.13 Tanner code (N = 155) decoding by BP on AWGNC . . . . . . . 2.1 Spin glass: the dashed lines are the anti-ferromagnetic interactions and the solid lines are the ferromagnetic interactions, the arrows are the spins (up or down) . . . . . . . . . . . . . . . . . 2.2 Spin glass: dashed lines are anti-ferromagnetic interactions and solid lines are ferromagnetic interactions; arrows are spins (up or down). The red line is an interaction modified by an external noise and red spins are spins modified by interactions. . . . . . 2.3 Ising chain with N = 3 spins embedded in an external magnetic field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 N = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 N = 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 N = 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 p-range model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3. 1

 1 Logistic map x k+1 = 4x k (1x k ) . . . . . . . . . . . . . . . . . . 3.2 Logistic map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Lorenz attractor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Lorenz map . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Logistic map -Period doubling . . . . . . . . . . . . . . . . . . . 3.6 Bifurcation diagram of the logistic map . . . . . . . . . . . . . . 3.7 Poincaré section P . . . . . . . . . . . . . . . . . . . . . . . . . .

  4) where the notation < ., . > means neighboring sites. The spins variables {s i } i are Ising variables i.e. they equal ±1 (only pairwise interactions are considered but we will see in our research that are often needed high order interactions). If the terms {J ij } ij are random variables, usually Gaussian variables, then the set of spins represent an Edwards-Anderson model [M. 06]. The infinite range model is called the Sherrington-Kirkpatrick model [D. 75].

Chapter 4 - Bethe approximation: an approach to the Belief Propagation

§ 4.1 Introduction BP is a well-known algorithm first introduced by

Pearl [Pea88]

in 1988 to solve inference problem on Bayesian networks in the framework of artificial intelligence. It swiftly turned out that this method could be applied to many research fields with similar issues. As a matter of fact, problems modeled by factor graphs, Markov Random Fields (MRF), or any other probabilistic graphical models, can be solved or at least approximately solved by BP, e.g. neural networks [T. 06b], image processing and video processing [M. 02], spin glasses [M. 06], channel coding

[D. 95]

. For each of these examples, the goal is to extract either an approximate of the most likely state or approximates of marginals. As a very important example developed in this chapter, it was proved in the 90's that BP had deep links with statistical physics where problems were oriented towards a constrained optimization of the free energy.

  1)with E • the floor function. is close to an affine behavior in d min . Therefore, as d min is increased, the corresponding code is able to correct more and more errors. Coding theory is aimed at designing error correcting codes which d min is as large as possible. For very structured codes, defined as codes which construction is easily mapped for different values of N, a constraint is that d min has to evolve the most linearly possible with N, that offers very good codes. However, computing d min is not a trivial task as it is closely related to the rank of the parity-check matrix whose length N can reach very large values. That is partly why LDPC codes design is still an open problem.

  x 4 , x 5 , x 7 , x 8 )p 235689 (x 2 , x 3 , x 5 , x 6 , x 8 , x 9 ) p 258 (x 2 , x 5 , x 8 ) ≈ p 1245 (x 1 , x 2 , x 4 , x 5 )p 2356 (x 2 , x 3 , x 5 , x 6 )p 4578 (x 4 , x 5 , x 7 , x 8 )p 5689 (x 5 , x 6 , x 8 , x 9 )p 25 (x 2 , x 5 )p 45 (x 4 , x 5 )p 56 (x 5 , x 6 )p 58 (x 5 , x 8 )

	(5.2)
	Considering BP assumption 2 , the Markov condition offers us that:
	p(x)

p 5 (x 5 )

Table 6 . 1 :

 61 .3 the iterative evolution of E(k) as a function of the SNR, considering BP decoding on EE. Four typical behaviors may be distinguished, ordered in the table 6.1. BP behaviors in terms of SNR values These four behaviors are not specific to EE, we encountered them for most other non-trivial error events. Difference remains in the SNR values where the behavior changes. We note that continuously increasing SNR values does not make BP behavior change as continuously, but abruptly for specific SNR values. We call these SNR values critical SNR values. To study rough behavior changes, it is relevant to use a bifurcation diagram.We display on Fig.6.4 the iterative evolution of E(k) on EE as a function of the SNR with GBP decoding.

		1
		0.95
	E(k)	0.85 0.9
		0.8			2.00dB	2.30dB
					2.90dB	3.09dB
		0.75	0	100	200	300	400
					k
	Figure 6.3: BP MSB of EE: four SNR values, to each one is found a typical behavior

6.5.2.1 BP Bifurcation

  diagram of BP for the selected error event EE is displayed in Fig.6.5. Five obvious SNR intervals may be emphasized in this figure, each one corresponding to a specific BP behavior. We summarize corresponding SNR intervals in table 6.2.

		1
		0.95
	E(K)	0.9
		0.85
		0.8	2	2.2 2.4 2.6 2.8	3	3.2
				SNR (dB)
	Figure 6.5: Bifurcation diagram of BP
	Behavior of E(K)	SNR
	smooth increasing	[0.00 dB; 2.20 dB]
	pseudo-oscillations	[2.20 dB; 2.51 dB]
	erratic evolutions	[2.51 dB; 2.99 dB]
	convergence jumps	[2.99 dB; 3.04 dB]
	optimal convergence	[3.04 dB; ∞ dB]

Table 6 . 2 :

 62 BP behaviors in terms of SNR intervalsNumerous other experiments, performed with other non-trivial error events, but not displayed here, made reveal that lengths and bounds of these intervals are quite constant. In addition, these behaviors are always ordered in the same way, which indicates a regularity in BP behaviors for Tanner code. When changing LDPC code, BP similarly behaves but with different lengths and bounds of the SNR values, e.g. McKay code of length N = 504 which bifurcation diagram is exhibited in Fig.6.6. Even though critical SNR values and amplitudes of E(K) are different from Tanner code, the five intervals are clearly distinguishable.

		1
		0.98
	E(K)	0.96
		0.94
		0.92	1	1.2	1.4
				SNR (dB)

  6.13 reduced trajectories of BP and GBPs for SNR ∈ {2.00, 2.30, 2.90, 3.09}dB, respectively.

	0.94			0.94
	0.9			0.9
	0.86			0.86
	0.94 0.82			0.94 0.82
	0.9			0.9
	0.82 0.86 0.82 0.86	0.9	0.94	0.82 0.86 0.82 0.86	0.9	0.94
	(a) BP			(b) GBP-C
	0.94			0.94
	0.9			0.9
	0.86			0.86
	0.94 0.82			0.94 0.82
	0.9			0.9
	0.82 0.86 0.82 0.86	0.9	0.94	0.82 0.86 0.82 0.86	0.9	0.94
	(c) GBP-P			(d) GBP-IP

Figure 6.10: Reduced trajectory with SNR = 2.00dB

Table 6 . 3 :

 63 Coordinates, lengths and complementary lengths

a posteriori probability is originally defined as p(x|y). As channels are assumed to be memoryless, likelihoods do not depend on the LDPC code used to encode information, i.e. they only help get an insight on possible values of bits x 1 , . . . , x N . Bayes rule stipulates that p(x, y) = p(x|y)p(y), therefore provided observations y 1 , . . . , y N , searching for distribution p(x, y) is equivalent to searching for distribution p(x|y). That is why we allow us to call a posteriori distribution both distributions.

see documentation on DVB-S2 at http://www.etsi.org/deliver/etsi_en/302300_302399/

302307/01.02.01_60/en_302307v010201p.pdf

denoting H = [h 1 . . . h N ] T and x = [x 1 . . . x N ] T , we obtain that E H (x) = -∑ N i=1 h i x i

channels studied in this thesis (AWGNC, BSC) are memoryless: p(y|x) = ∏ N i=1 p i (y i |x i )

demonstrations of (4.39) and (4.41) implicitly assume that any marginal distribution, either on a potential or on a variable node, is exact to the associated marginalization on the joint distribution, see next section for details

any variable node X i of a unit degree d i = 1 does not appear in(4.42). Therefore, Lagrangian formalism cannot offer any equation of associated belief. Such a belief is obtained according to constraint (C3).

strictly speaking, the trapping sets definition does not depend on BP. As a matter of fact, these are not eventually corrected by any arbitrary iterative decoder, as Gallager-B, Bit-flipping,

BP, etc. In this thesis, we focus on failures of BP, that is why we orient definition towards BP

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 X 9 X 4 , X 5 , X 7 , X 8 X 5 , X 6 , X 8 , X 9 X 1 , X 2 , X 4 , X 5 X 2 , X 3 , X 5 , X 6 X 5 , X 8 X 4 , X 5 X 5 , X 6 X 2 , X 5 X 5

more precisely, couplings {J ij } <i,j> are not ignored but probabilistic correlations are, i.e. any < X i , X j > is approximated by < X i >< X j >

which is true, given that if N p ≤ 1, r cannot exist

which is the explicit representation of the square lattice, i.e. we exhibit the check nodes that belong to the edges

except if we select the region-graph made of a unique region being the Tanner graph itself

short loops are more harmful than large ones as the iterative aspect of BP brings out their effect in very few iterations

trapping sets TS(5, 3) are strongly connected in Tanner code

c a , c b , c d , c e , c g , c h X 1 , X 2 , X 3 , X 4 c b , c c , c e , c f , c h , c i X 2 , X 3 , X 4 , X 5 c a , c c , c d , c f , c g , c i X 1 , X 2 , X 4 , X 5 c b , c e , c h X 2 , X 3 , X 4 c a , c d , c g X 1 , X 2 , X 4 c c , c f , c i X 2 , X 4 , X 5 X 2 , X 4 Figure 5.29: Region-graph resulting from the split of a TS(5, 3) in quadrupletsThe second split of trapping sets in triplets results in local region-graphs, see Fig.5.30, that satisfied the local loopfree principle. This way, RBA is locally exact, i.e. the GBP algorithm provides local optimal performance. In addition, triplets only contain 13 variable nodes, a message would then need 2 13 = 8192 elementary operations, that is considerably less than quadruplets.c a , c d , c g X 1 , X 2 , X 4 c b , c e , c h X 2 , X 3 , X 4 c c , c f , c i X 2 , X 4 , X 5 X 2 , X 4Figure 5.30: Region-graph resulting from the split of a TS(5, 3) in triplets

we recall that a check node has to be accompanied by its neighborhood, as a region is a subgraph of the Tanner graph

any decoder is stopped as soon as it has converged

for BP, a BER of 0.03

on BSC, forcing a bit flipping is also swapping likelihoods

but we obtain a subset of these genuine random realizations, i.e. they are still Gaussian variables

BER is averaged over 10 4 simulations, which is enough to obtain reliable results that reflect behaviors of the decoders with our protocol

we will deal with this peculiar phenomenon in next chapter

for any variable node X i , b i (0, y i ) is enough to describe the distribution of X i as b i (1, y i ) = 1b i (0, y i )

strictly speaking, defining dynamical systems with beliefs is not correct as update rules G BP and G GBP are not bijective, i.e. many different beliefs vectors at iteration k -1 may produce the same beliefs vector at iteration k

if X and Y are two random variables s.t. Y = √ aX, then σ X = aσ Y

a realization, the corresponding error event and the likelihoods are said non-trivial when decoding results in a non-trivial BER

thanks to the decreasing damping factor

E(K) is assumed to be the steady state

for middle SNR values

larger values of K would make the comparison easier but as pointed out earlier, it would not reflect a realistic point of view

the line drawn by oscillations is called a limit cycle

in the true state space, not in the reduced space

averaged over a thousand simulations

for any variable node X i , b i (0) = 1

as GBP-P and GBP-IP present pseudo-oscillations

state space made with beliefs is a hypercube which largest diagonal length is √ N, i.e. around 12.4 for Tanner code

Remerciements

Je remercie les personnes qui ont contribué à mon entrée dans le monde de la recherche. Parmi elles, Sylvain Reynal et David Declercq, qui m'ont pro-

Part I

Overview: methods and tools

Part II

Graphical models and Belief Propagation algorithms § 6.2 Preliminaries BP and GBP are dynamical systems that sometimes behave very differently when initialization parameters are changed. Even though iterated maps that update messages are deterministic, we observe significant dependence on likelihoods, see Fig. 6.1. We represent in Fig. 6.2 the space generated by likelihoods of two variable nodes (X i = 100, X j = 102), randomly selected in Tanner code. For a given channel realization at SNR = 2.3dB, we quench all likelihoods of other variable nodes and we make vary likelihoods l i (0, y i ) and l j (0, y j ) from their minimum value, zero, to their maximum value, one. We then run a BP decoding for each couple of values {l i (0, y i ) ∈ [0; 1], l j (0, y j ) ∈ [0; 1]} and we pick out the BER to plot it in the likelihoods space. We observe that two neighboring couples can reach very different BER, depending on their relative position in the space. As an example, points A and B are defined such that variable nodes X i and X j are:

Computing BER variances σ 2 A , σ 2 B in disks of radius 0.1, centered on A and B, respectively, offers that: σ 2 A = 8.10 -6 , σ 2 B = 3.10 -4 .