Raouia Masmoudi

E Veronica Belmega

Inbar Fijalkow
email: inbar.fijalkow@ensea.fr

Noura Sellami
email: sellaminoura@yahoo.fr

A Closed-Form Solution to the Power Minimization Problem over Two Orthogonal Frequency Bands under QoS and Cognitive Radio Interference Constraints

Keywords: Cognitive Radio channels, power-efficient spectrum allocation problems, green communications

This paper considers a Cognitive Radio (CR) channel composed of a secondary user (SU) and a primary user (PU). An analysis of the power minimization over several orthogonal frequency bands at the SU level under the following constraints is provided: a minimum Quality of Service (QoS) constraint, maximum peak and average interference to the PU constraints. The general solution, when it exists, is a water-filling type of solution which can be computed via iterative algorithms. It turns out that, in the case of two orthogonal bands a closed-form analytical solution can be found and a complete analysis of the feasibility of these opposing constraints is presented in details. Several numerical results that sustain and give inside into the analysis are also discussed.

I. INTRODUCTION

The concept of Cognitive Radio [START_REF] Mitola | Cognitive radio architecture evolution[END_REF] has recently emerged as a promising paradigm for a more efficient use of the available spectrum by allowing the coexistence of licensed (primary) and unlicensed (secondary) users in the same bandwidth. Another major concern that is rapidly gaining momentum is the energy consumption for green communications. Recently, there has been an important switch of focus from a data rate maximization vision to an application-based utility optimization one, that takes into account the cost of power consumption to achieve the target data rates [START_REF] Belmega | An information-theoretic look at MIMO energy-efcient communications[END_REF], [START_REF] Belmega | A survey on energyefficient communications[END_REF].

In this paper, we consider a home-automation scenario where several technologies (i.e., WiFi, PLC, Femto, etc.) are able to operate simultaneously and where different appliances may have different hierarchical priorities. The EconHome-box is meant to allocate frequency bands and powers in a green telecommunications sprint. Notice that the analysis we provide is not restricted to this scenario alone and can be applied to any CR system with several orthogonal sub-channels available This work was sponsored in part by the FUI SYSTEM@TIC Project EconHome.

for communication: OFDM systems, interoperability of several coexisting transmission systems.

To model this home-automation scenario we consider a CR channel composed of a primary user (PU) and a secondary user (SU) that is opportunistically accessing the channel and interfering with the PU. The question under investigation is to know how the SU will optimally allocate the available power over the available bands to minimize its total power consumption under several constraints: target QoS constraint, peak and average maximal interference levels that the SU can create to the PU.

The problem of resource allocation has mainly been studied from a data rate maximization point of view. For exemple, in [START_REF] Cendrillon | Optimal multi-user spectrum managment for digital subscriber lines[END_REF], [START_REF] Scutari | Asynchronous iterative water-filling for gaussian frequency-selective interference channels[END_REF], [START_REF] Bagayoko | Power control of spectrumsharing in fading environment with partial channel state information[END_REF] the authors considered the maximization of the Shannon achievable rate over several sub-channels, subject to transmit powers and mask constraints on each sub-channel. In [START_REF] Cendrillon | Optimal multi-user spectrum managment for digital subscriber lines[END_REF], a centralized approach based on convex optimization (see e.g., [START_REF] Bazaraa | Nonlinear programming-Theory and Algorithms[END_REF], [START_REF] Boyd | Convex Optimization[END_REF]) was proposed to compute, under some technical conditions, the achievable rate region of the system.

To the best of the authors' knowledge, the closest works to ours are [START_REF] Scutari | Competitive Optimization of Cognitive Radio MIMO Systems via Game Theory[END_REF] and [START_REF] Pang | Distributed power allocation with rate constraints in gaussian parallel interference channels[END_REF]. In [START_REF] Pang | Distributed power allocation with rate constraints in gaussian parallel interference channels[END_REF], the authors consider the distributed minimization of transmit powers subject to rate target constraints for the multi-user parallel interference channels. The difference to our work consists in the absence of PU and thus, of the interference constraints that limit the SU transmissions'. In [START_REF] Scutari | Competitive Optimization of Cognitive Radio MIMO Systems via Game Theory[END_REF], the authors study the multiusers multiple-input multiple-output (MIMO) Cognitive Radio channels. In order to design cognitive MIMO SU transceivers, the authors of [START_REF] Scutari | Competitive Optimization of Cognitive Radio MIMO Systems via Game Theory[END_REF] model the distributed maximization of the Shannon achievable rates under interference constraints imposed by the presence of the PU. In order to approach this complex problem, the authors study several particular cases by systematically eliminating one of the interference constraints and solving a simpler problem. In this paper, we consider a simpler channel model (the single-SU case and the multiple orthogonal sub-channels). However, the novelty and interest of our work consists in the fact that we study the dual problem of [START_REF] Scutari | Competitive Optimization of Cognitive Radio MIMO Systems via Game Theory[END_REF] in its generality with respect to the feasibility constraints.

Our contributions can be summarized as follows: i) in the general case of arbitrary number of orthogonal bands, we give the analytical water-filling solution to the power minimization problem under minimum rate and maximum peak and average interference constraints created to the PU, when it exists; ii) in the case of two orthogonal frequency bands, we give an analytical closed-form solution and a complete study of the feasibility of all the constraints is provided; iii) numerical simulations that sustain the analysis are provided for different channel parameters.

The remainder of this paper is organized as follows. In section II, we describe the system model. We provide the water-filling type of solution of the general optimization problem section III. A closed-form solution in the case of two orthogonal bands is given in section IV. In section V, we present some interesting simulation results. We conclude this paper in sectionVI.

II. SYSTEM MODEL

The CR channel under study system is depicted in Fig. 1. The primary/secondary user consists of a Primary/Secondary Transmitter (PT/ST) and a Primary/Secondary Receiver (PR/ST) respectively. Each transmitter or receiver is equipped with only one antenna. The transmission is done over N ≥ 2 orthogonal frequency bands. The transmit power of ST in the frequency band k ∈ {1, ...N } is denoted by p k . We denote the overall power allocation profile by p = (p 1 , p 2 , ..., p N ), p ∈ R N + . The received signal at SR in band k can be written as:

y k = √ p k h k s k + i k,P U + n k , (1)
where h k represents the channel gain, s k the transmit signal of the ST is given by unit-power variable over each orthogonal frequency band. The instantaneous power gain of the ST-SR direct link, the interfering ST-PR link is denoted by are denoted by |h k | 2 , g k respectively. All the channels are assumed to be stationary, ergodic and independent from the noises. The noise n k N (0, σ 2 k ) is a zero-mean circularly symmetric complex Gaussian noise vector and the interfering signal from the PU i k,P U N (0, σ 2 k,P U ). In this scenario, the Shannon achievable rate of the SU is given by

R(p) = N k=1 log 2 (1 + c k p k ), (2)
where c k is related to the interfering ST-PR link by

c k = |h k | 2 σ 2 k + σ 2 k,P U (3)
where σ 2 k,P U is the interfering variance from the PU and σ 2 k is the variance of the thermal noise.

III. GENERAL OPTIMIZATION PROBLEM

We consider the constrained power minimization problem at the SU level. Three constraints are considered in our optimization problem: a) a QoS constraint (4); b) peak [START_REF] Bagayoko | Power control of spectrumsharing in fading environment with partial channel state information[END_REF] and average (5) interference constraints to protect the PU and limit the interference caused by the SU. 1-Target rate constraint:

R(p) = N k=1 log 2 (1 + c k p k ) ≥ R min , (4)
which is imposed to achieve a minimum QoS of the SU transmission. 2-Average interference power shaping constraint:

N k=1 g k p k ≤ P , ( 5
)
where P is the maximum average interference level that can be received at PR. This constraint is an imposed interference limitation at the PR. 3-Peak interference power shaping constraints:

0 ≤ g k p k ≤ P peak k , ∀ k = 1, . . . N, (6)
where P peak k is the maximum peak interference level that can be received at PR in a given band.

The optimization problem under study can be written as:

                           min {p k } N k=1 N k=1 p k s.t. R(p) = N k=1 log 2 (1 + c k p k ) ≥ R min N k=1 g k p k ≤ P 0 ≤ g k p k ≤ P peak k ∀ k ∈ {1, ..N } (7)
Notice that without loss of generality, we can consider only the case where R min > 0. In the case where R min = 0 the solution is trivial and the transmission powers are equal to zero, p * k = 0 ∀k, i.e.,no transmission is performed at the SU.

In the case of applying a flat power allocation policy (i.e., ∀k ∈ {1, ..N } p k = p), all the constraints become more strict and the feasible set is reduced. The constraints in this case will be given by p ≤ min

P peak k g k , P k g k and k log(1 + c k p) ≥ R min .
Therefore, the benefit of applying our optimal allocation policy comparing to this flat power allocation policy, is to have much more flexibility in the feasibility of our constraints and a significant gain in terms of power consumption.

Theorem 1: If the feasible set of the optimization problem [START_REF] Bazaraa | Nonlinear programming-Theory and Algorithms[END_REF] is not void, the optimal power allocation policy is given by the following water-filling type solution 1 :

p * n = ξ n - 1 c n P peak n g n 0 , ∀n ∈ {1, . . . N } ( 8
)
where

ξ n = λ ln(2)(1 + βg n )
, λ and β verify the following conditions: If the feasible set is void then the problem (7) has no solution.

λ > 0 and N k=1 log 2 (1 + c k p * k ) = R min
Proof: We denote the feasible set by S N such that:

S N = p ∈ R N + | 0 ≤ g k p k ≤ P peak k , N k=1 log 2 (1 + c k p k ) ≥ R min , N k=1 g k p k ≤ P . (9)
Two different scenarios are possible:

• S N = ∅: In this case, there is no feasible power allocation policy, i.e., no vector p ∈ R N + that satisfies simultaneously all the constraints in (9).

• S N = ∅: Given the concavity of log(.) function, the linearity of powers constraints and the convexity of the 1 We denote by [x] b a = min{max{x, a}, b} or equivalently

[x] b a =    x, if a ≤ x ≤ b a, if x < a b, if x > b.
objective function

N k=1 p k , it is easy to prove that S N is convex [8], i.e., ∀ p k , pk ∈ S N , ∀ θ ∈ [0, 1] we have θ p k + (1 -θ) pk ∈ S N . (10)
Thus, the optimization problem in ( 7) is a convex problem which has at least an optimal solution. Therefore, we can apply the Lagrangian method to solve [START_REF] Bazaraa | Nonlinear programming-Theory and Algorithms[END_REF]. The Lagrangian is given by:

L(p, λ, β, α, µ) = N k=1 p k -λ N k=1 log 2 (1 + c k p k ) -R min ) + β N k=1 g k p k -P - N k=1 α k g k p k + N k=1 µ k g k p k -P peak k , (11)
where λ, β, α = (α 1 , α 2 , ..., α N ) and µ = (µ 1 , µ 2 , ..., µ N ) are the non-negative Lagrangian multipliers of the corresponding constraints (4)(5)(6), respectively. We denote p * an optimal solution. The Karush-Kuhn-Tucker (KKT) optimality conditions are:

i) Compute the solution that gives null Lagrangian derivative:

∂L ∂p n pn=p * n = 0 ;
ii) Non-negative Lagrangian multiplier of the rate constraint λ > 0 and

N k=1 log 2 (1 + c k p * k ) = R min or λ = 0 and N k=1 log 2 (1 + c k p * k ) > R min ;
The case of null Lagrangian multiplier of the rate constraint cannot occur, because if λ = 0 then the Lagrangian in (11) becomes:

L(p, β, α, µ) = N k=1 p k + β N k=1 g k p k -P - N k=1 α k g k p k + N k=1 µ k g k p k -P peak k ,
and the optimization problem in ( 7) is given by:

               min {p k } N k=1 N k=1 p k s.t. N k=1 g k p k ≤ P 0 ≤ g k p k ≤ P peak k ∀ k ∈ {1, ..N }
In this power optimization problem the optimal power solution p * is null. Then the SU optimal rate is null also (i.e., N k=1 log 2 (1 + c k p * k ) = 0), this is contradictory to the QoS constraint which is strictly positive (R min > 0). Thus, the Lagrangian multiplier of the rate constraint is strictly positive and given by: λ > 0 and As we said before, in order to compute the optimal solution of our problem we derive the Lagrangian in (11),

∂L ∂p n pn=p * n = 0, ⇔ 1 - λc n ln(2)(1 + c n p * n ) + (β -α n + µ n )g n = 0, ⇔ λc n ln(2)(1 + c n p * n ) = 1 + (β -α n + µ n )g n , ⇔ (1 + c n p * n ) = λc n ln(2)(1 + (β -α n + µ n )g n ) , ⇔ p * n = λ ln(2)(1 + (β -α n + µ n )g n ) -1 cn .
Considering the other KKT conditions, the solution in function of the optimal (α n ,µ n ) is as follows: α n = 0 and µ n = 0 (the optimal power solution is strictly 2 between 0 and P peak n /g n )

p * n = λ ln(2)(1 + βg n ) - 1 c n P peak n g n 0 , (12)
α n = 0 and µ n > 0, 2 The notation (x) b a is defined by:

p * n = P peak n g n , (13)
(x) b a = x, if a < x < b 0, otherwise.
α n > 0 and µ n = 0,

p * n = 0. (14)
Finally, the equations ( 12), ( 13) and ( 14) can be re-written in a compact form as follows:

p * n = λ ln(2)(1 + βg n ) - 1 c n P peak n g n 0 , (15)
where λ and β verify the following conditions:

λ > 0 and N k=1 log 2 (1 + c k p * k ) = R min β > 0 and N k=1 g k p * k = P or β = 0 and N k=1 g k p * k < P .
Thus, by making a change of variable

ξ n = λ ln(2)(1 + βg n )
, we obtain the analytical solution of the problem [START_REF] Bazaraa | Nonlinear programming-Theory and Algorithms[END_REF], which is expressed by the well known water-filling:

p * n = ξ n - 1 c n P peak n g n 0 . ( 16
)
When the feasible set is void, one or more constraints need to be relaxed if the system owner wishes to allow a non-trivial secondary transmission.

IV. THE CASE OF TWO ORTHOGONAL FREQUENCY BANDS

In this section, we focus on the specific case where only two orthogonal frequency bands are available, i.e., N = 2. This is motivated by the fact that an analytical study of the feasible set in the general scenario ( 9) is very difficult. As we will show, this difficulty is overcome in this case. We provide a complete and closed-form analysis of the feasibility of the opposing constraints: the QoS constraint on one hand, and the power constraints on the other. Moreover, this two bands case models the systems where two different technologies can coexist, i.e., (WiFi+PLC), or (WiFi+Femto), which is pertinent in current home-automation systems. Because these technologies are inherently different or heterogeneous, it is realistic to assume the orthogonality property of the sub-channels.

In what follows, we focus on the simpler case of N = 2 for the following reasons: a) it can be realistic in several scenarios among which we have mentioned here above two of them; b) it can be solved in an analytical closed-form way; c) it gives insight and intuition on what happens in a general system with N > 2 sub-channels which can be solved only numerically.

In this case, the optimization problem in ( 7) can be reformulated as:

min{p 1 + p 2 } s.t p ∈ S 2 , (17)
where, the feasible set S 2 is given by

S 2 = p ∈ R 2 + | 0 ≤ g 1 p 1 ≤ P peak 1 , 0 ≤ g 2 p 2 ≤ P peak 2 , 2 k=1 log 2 (1 + c k p k ) ≥ R min , 2 k=1 g k p k ≤ P .
Theorem 2: In the case of two orthogonal frequency bands, when a solution to (17) exists, it is unique and its closed-form solution is given here under.

If P < - g 1 c 1 - g 2 c 2 + 2 ( g 1 c 1 g 2 c 2
) exp(R min ) or otherwise if max{ p 1 , p 1 } < 0 then the constraints ( 4) and ( 5) can not be simultaneously satisfied. Also if at least one of the following conditions is met: 6) can not be satisfied. In all the aforementioned cases, a feasible solution does not exist. In all other possible scenarios an optimal solution exists and is given by

P peak 1 g 1 < p 1 or P peak 2 g 1 < min{ p 2 , p2 } or f 1 ( P peak 1 g 1 ) > P peak 2 and f -1 1 ( P peak 2 g 2 ) > P peak 1 g 1 then constraint (
                   p * 1 =   -1 c 1 + exp(R min ) c 1 c 2   y x , p * 2 =   -1 c 2 + exp(R min ) c 1 c 2   z w
, where x, y, w, z are given by

                               x = max max { p 1 , 0} , f -1 1 P peak 2 g 2 , y = min min p 1 , f -1 1 (0) , P peak 1 g 1 , w = max max { p 2 , 0} , f 1 P peak 1 g 1 , z = min min { p 2 , f 1 (0)} , P peak 2 g 2 . ( 18
)
The parameters p, p are the intersection points of the functions f 1 (.) and f 2 (.) defined as: f 1 : R→R (derived from the rate constraint(4)) such that:

f 1 (p 1 ) = 1 c 2 exp(R min ) 1 + c 1 p 1 -1 (19)
and f 2 : R → R (derived from the average interference power constraint( 5)) such that:

f 2 (p 1 ) = P -p 1 g 1 g 2 . ( 20
)
Proof: In order to solve the optimization problem (17), we rewrite the feasible set using the two functions defined in ( 19) and (20) in the following manner:

S 2 = A B A = (p 1 , p 2 ) ∈ R 2 | f 1 (p 1 ) ≤ p 2 , f 2 (p 1 ) ≥ p 2 } B = 0, P peak 1 g 1 × 0, P peak 2 g 2 . ( 21
)
The set B represents simply the peak interference constraints in each band which restricts the feasible set S 2 . Notice that the set B is not void, however the set A can be void.

In order to find the conditions under which the set S 2 is void or not, i.e., a solution exists or not, we have first to discuss the set A. Thus, we analyse the feasible set A in the plane p 1 Op 2 , where O is the origin point of coordinates (0,0) and p 1 and p 2 are the axes of our plane. In order to do that, we investigate the intersection points (p 1 , p 2 ) ∈ A, of the functions f 1 (.) and f 2 (.) by solving the system of equations:

                           p 2 = 1 c 2 exp(R min ) 1 + c 1 p 1 -1 , p 2 = P -g 1 p 1 g 2 , p 1 ∈ 0, P peak 1 g 1 , (22)
We further have:

               P -g 1 p 1 g 2 = 1 c 2 exp(R min ) 1 + c 1 p 1 -1 , p 1 ∈ 0, P peak 1 g 1 , (23)
⇔ P g 2 - g 1 g 2 p 1 = 1 c 2 exp(R min ) (1 + c 1 p 1 ) - 1 c 2 , ⇔ P c 2 -g 1 c 2 p 1 + g 2 (1 + c 1 p 1 ) = g 2 exp(R min ).
We obtain a quadratic equation w.r.t p 1 Thus, we obtain in (24) a quadratic equation where its discriminant is given by

-g 1 c 2 c 1 p 2 1 + P c 1 c 2 -g 1 c 2 + g 2 c 1 p 1 +P c 2 + g 2 (1 -exp(R min )) = 0. (24)
∆ p1 = P c 1 c 2 -g 1 c 2 + g 2 c 1 2 +4g 1 c 2 c 1 (P c 2 + g 2 -g 2 exp(R min )). ( 25
) (a) If ∆ p1 < 0 we have that f 1 (p 1 ) > f 2 (p 1 ) ∀p 1 ∈ R, this
means that the functions f 1 (.) and f 2 (.) do not intersect i.e the two functions in this case are disjoint (Fig. 2) and the set A = ∅. Thus, S 2 = ∅ and we have no solution.

∆ P = 16c 3 1 c 3 2 g 1 g 2 exp(R min ) > 0. (26)
Thus, we obtain two solutions:

P 1,2 = - g 1 c 1 - g 2 c 2 ± 2 g 1 g 2 c 1 c 2 exp(R min ). (27)
Note that one of the two solutions is negative, so it is rejected. Thus, the constraint of average interference powers must satisfy this condition:

P = - g 1 c 1 - g 2 c 2 + 2 g 1 g 2 c 1 c 2 exp(R min ). ( 28
)
Thus the intersection point in this case (∆ p1 = 0) where

f 1 (.) is tangent to f 2 (.), is unique p = p with            p 1 = P 2g 1 p 2 = P 2g 2 (
P > - g 1 c 1 - g 2 c 2 + 2 g 1 g 2 c 1 c 2 exp(R min
), then we have two distinct intersection points (see Fig. 4 and Fig. 5) ( p 1 , p 2 ), ( p1 , p2 ) ∈ A, where

p 1 = (P c 1 c 2 -g 1 c 2 + g 2 c 1 ) -∆ p1 2g 1 c 1 c 2 , p 1 = (P c 1 c 2 -g 1 c 2 + g 2 c 1 ) + ∆ p1 2g 1 c 1 c 2 , ( 30
)
where ∆ p1 is expressed in (25), p 1 < p1 , p 2 = f 1 ( p 1 ) and p2 = f 1 ( p1 ) . • If p 1 or p 1 is negative, (see Fig. 4 (a) when p 1 < 0), then we have an unique solution if A ∩ B = ∅.

In this case, see Fig. 4(a), the optimal solution (p * 1 , p * 2 ) of the problem (17) is not between the intersection points ( p 1 , p 2 )and ( p 1 , p 2 ), but it is computed between ( p 1 , p 2 ) and (0, f 1 (0)). These two last cases (if max{ p 1 , p 1 } < 0 and if p 1 or p 1 is negative) are met, only for small values of R min .

• If both of the two abscissas of the intersection points p 1 and p 1 are positive, (see Fig. 5), then we have an unique solution if A ∩ B = ∅ (see Fig. 5(a)) and in this case, the optimal solution (p * 1 , p * 2 ) exists. However, if A ∩ B = ∅, see Fig. We give a short sketch of this analysis in Summary 1.

Summary 1 Closed-form analytical solution for the N = 2 case

If P < - g 1 c 1 - g 2 c 2 + 2 g 1 c 1 g 2 c 2 exp(R min ), then A = ∅:
No Solution: Average interference power P and QoS constraints are never simultaneously satisfied

Else A = ∅, If max{ p 1 , p 1 } < 0 ⇔ A ⊂ (-∞, 0] × [0, +∞), then A ∩ B = ∅:
No Solution: Average interference power P and QoS constraints are not simultaneously satisfied in B Exists a solution given by: In Fig. 11, we plot the optimal rate R * and SU rate and interference powers caused by SU to PU (g 1 p * 1 , g 2 p * 2 and g 1 p * 1 + g 2 p * 2 ) versus the interference link gain g 1 over two orthogonal frequency channels with the system parameters: In Fig. 11.(a), the optimal rate is equal to the minimum rate constraint R min for values of the interference link gain g 1 below a certain threshold. We see that when the quality of the interfering link is high (i.e. high values of g 1 ) the maximal interference constraint are not met and a feasible solution does not exist.

Else max{ p 1 , p 1 } > 0 ⇔ A ⊂ (-∞, 0] × [0, +∞), If P peak 1 g 1 < p 1 or P peak 2 g 2 < min{ p 2 , p2 } or f 1 P peak 1 g 1 > P peak
                   p * 1 =   -1 c 1 + exp(R min ) c 1 c 2   y x , p * 2 =   -1 c 2 + exp(R min ) c 1 c 2   z w . EndIf EndIf EndIf
In depend on g 1 . The optimal transmit power p * 2 is independent of the interference link gain g 1 except for the same region where p * 2 is proportional to the interference link gain g 1 (see equation (18)). For example, if g 1 1.7, we see that the average interference constraint is satisfied with equality. However, if c 1 < 1.7 the average interference constraint is satisfied with strict inequality i.e. g 1 p * 1 + g 2 p * 2 < P .

In the case where N > 2, the difficulty lies in the fact that the iterative water-filling algorithm converges to the optimal solution provided that the feasibility set of all the constraints is non void. Finding the conditions on the system parameters that ensure that the feasible set is non void, is a difficult problem since the number of constraints and the degrees of freedom becomes larger.

VI. CONCLUSIONS

In this paper, we have analysed the power allocation problem over orthogonal bands at the opportunistic user's level under a target QoS constraint and maximal interference caused to the primary user of the spectrum constraints.

We have seen that, depending on the system parameters, a feasible solution does not always exists. This is mainly caused by the simultaneous requirements of minimum target achievable rate at the SU and maximum allowed interference level caused to the PU. When a solution exist, it takes the form of water-filling and the closed-form solution is provided for the particular case of two available frequency bands.

Future work include the case of multiple secondary and primary users for which efficient distributed algorithms have to be proposed. Also, we will consider the effect of imperfect parameters estimation or knowledge, and investigate the robustness of our power allocation problem to these imperfections.

Fig. 1 .

 1 Fig. 1. Channel model for the communication in band k.

N k=1 log 2 ( 1 +

 21 c k p * k ) = R min iii) Non-negative Lagrangian multiplier of the average interference power constraint β > 0 and N k=1 g k p * k = P or β = 0 and N k=1 g k p * k < P ; iv) Non-negative Lagrangian multiplier of the peak powers constraints for each band: α n > 0 and g n p * n = 0 or α n = 0 and g n p * n > 0 µ n > 0 and p * n = P peak n g n or µ n = 0 and p * n

Fig. 2 .

 2 Fig. 2. No Solution because the set A = ∅ ⇒ S2 = ∅

Fig. 3 .

 3 Fig. 3. ∆p 1 = 0 and the function f 2 (.) is tangent to the function f 1 (.) where (a) : One Solution, (b): No Solution

Fig. 4 .

 4 Fig. 4. Particular cases: either (a): One solution or (b):No solution

Fig. 5 .

 5 Fig. 5. ∆p 1 > 0, A = ∅ either (a):One Solution, or (b):No Solution The discussion in this case is done as a function of the sign of the intersection points ( p 1 , p 2 ) and ( p 1 , p 2 ). • If max{ p 1 , p 1 } < 0 (see Fig.4(b)), then we are back to the case f 1 (p 1 ) > f 2 (p 1 ), ∀ p 1 ∈ R. Thus, S 2 = ∅ and we have no solution.

  5(b), because of the peak interference constraint P peak k , then we have no solution.

1 ,

 1 then A ∩ B = ∅: No Solution: P peak k constraints are not satisfied Else A ∩ B = ∅:

1 Fig. 10 .

 110 Fig. 10. SU optimal rate (a) and interference levels caused by SU to PU (b) vs. the direct link gain c1

c 1 = 2 . 5 , 1 = 5 , P peak 2 = 7 ,

 1251527 c 2 = 1.4, R min = 4, P peak P = 10 and g 2 = 1.3.

P peak 1 g 1 1 Fig. 11 .

 11111 Fig. 11. SU optimal rate (b) and interference levels caused by SU to PU (b) vs. the interference link gain g1

In Fig. 7, we plot the SU optimal rate R * R(p * ) defined in (2)) and the interference powers levels caused by SU to PU (g 1 p * 1 , g 2 p * 2 and g 1 p * 1 + g 2 p * 2 ) versus the minimum rate constraint R min over two orthogonal frequency channels with the system parameters: c 1 = 2.5, c 2 = 1.4, P = 10, P peak In the two sub-figures, we observe that there is a zone where our problem does not admit a solution where the average interference power and QoS constraints are not simultaneously satisfied (for example the case of Fig. 2). However, there is a critical value when R min ≈ 3.866 below which an optimal solution exists. Therefore, if the QoS constraint is too restrictive, i.e., R min is above this threshold, a feasible solution does not exist.

In Fig. 7.(a), we see that the optimal rate is equal to the minimum rate constraint R min which is trivial, because of the 8 KKT condition of λ > 0 in (15). In Fig. 9.(a) the SU rate is equal to the minimum rate constraint R min for P > 6 . In conclusion, if P is below a certain threshold then the average interference constraint is too restrictive and a solution does not exist.

In Fig. 9.(b) we note that the interference power level caused by SU to PU for each orthogonal frequency band is constant and independent from P . Therefore, P does not influence the expressions of the optimal power p * once it exists. For example, if P 6, the average interference constraint is satisfied with equality. However, if P > 6 we obtain g 1 p * 1 + g 2 p * 2 < P , then the average interference constraint is satisfied with strict inequality. In Fig. 10.(a), the optimal rate is equal to the minimum rate constraint R min for all c 1 values (∀c 1 > 0.9). We see that if the quality of one of the direct links (i.e., small values of the direct link gain c 1 ) is very poor, then the QoS constraint cannot be satisfy and thus, a feasible solution does not exist.

In Fig. 10.(b), we note that the optimal power for each orthogonal frequency band decreases exponentially with the direct link gain c 1 . The higher the quality of the direct links, the lower are the powers required to achieve the QoS constraint. For example, if c 1 0.9, we see that the average interference constraint is satisfied with equality. However, if c 1 > 0.9 the average interference constraint is satisfied with


Ben Sarra

Andreas Abbès

Thomas Scheuermann

Meilender

Mathieu D'aquin

Sarra Ben Abbès
email: sarra.benabbes@lipn.univ-paris13.fr

Andreas Scheuermann
email: andreas.scheuermann@uni-hohenheim.de

Thomas Meilender
email: thomas.meilender@loria.fr

' Aquin
email: m.daquin@open.ac.uk

Characterizing Modular Ontologies

Keywords: modularization, patterns, modules, ontology

published or not. The documents may come

Introduction

Difficulties in reusing and maintaining large monolithic ontologies have resulted in an increasing interest in modularizing ontologies. In the past, several approaches and tools (e.g., SWOOP1 , NeOn Toolkit2 ) have been proposed to support the modularization of ontologies. Each of these approaches and tools respectively incorporates its own definition and notion of modular ontologies and criteria underpinning ontology modularization [?]. This proliferation is mainly due to the fact that the area of ontology modularization appears to be still in its infancy. Thus, it lacks the mature, well-defined, well-understood, and commonly agreed upon definitions and concepts as sociated with modularization in the area of software engineering [?]. This lack of knowledge about ontology modularization and, particularly, the lack of knowledge about characteristics of modular ontologies prevents its further development. Being able to characterize modular ontologies would allow for (1) comparing different modularization approaches, (2) assessing the quality of modular ontologies with respect to manually modularized ontologies, (3) customizing ontology modularization by providing modularization criteria, and (4) improving the area of ontology modularization as a whole. The goal of this work is to characterize modularly organized ontologies to contribute to a better understanding of ontology modularization. For this purpose, we (1) extract a number of existing ontologies from the web, (2) select and adopt metrics from software engineering in order to both (3) identify recurring structures (patterns) in modularly organized ontologies and (4) characterize the identified patterns. The contribution is a set of four patterns, which characterize modularly organized ontologies. The rest of this paper is organized as follows: Section 2 provides an introduction to ontology modularization and reviews related work. Section 3 reports on the research design whereas section 4 presents and discusses the results. Section 5 draws a conclusion and points to future avenues of research.

Ontology Modularization

Modular Ontologies

The main idea of modular ontologies originates from the general notion of modular software in the area of software engineering. Correspondingly, ontology modularization can be interpreted as decomposing potentially large and monolithic ontologies into (a set of) smaller and interlinked components (modules). Therefore, an ontology module can be considered as a loosely coupled and self-contained component of an ontology maintaining relationships to other ontology modules. Thereby, ontology modules are themselves ontologies [?]. In general, ontology modularization aims at providing users of ontologies with the knowledge they require, reducing the scope as much as possible to what is strictly necessary. In particular, ontology modules (1) facilitate knowledge reuse across various applications, (2) are easier to build, maintain, and replace, (3) enable distributed engineering of ontology modules over different locations and different areas of expertise, and (4) enable effective management and browsing of modules [?].

Approaches for Ontology Modularization

In recent years, the problem of ontology modularization has attracted more and more attention and, thus, several different approaches for modularizing ontologies appeared. These approaches can be classified in two main categories. The first main category comprises approaches that focus on the composition of existing ontologies by means of integrating and mapping ontologies. On the one hand, approaches addressing integration of existing ontologies are owl:import, partial semantic import, e.g., [?,?], package-based description logics, e.g., [?].

On the other hand, mapping approaches basically aim at (inter-)linking sets of ontology modules. The following approaches can be assigned to these two formalisms: distributed description logics, e.g., [?,?], ε-connections, e.g., [?,?]. Other approaches establish the relationship between various modular ontology formalisms [?] in order to have special syntax in the ontology languages for a modeling perspective. The second main category comprises approaches for modularizing ontologies in terms of ontology partitioning and ontology module extraction. On the one hand, ontology partitioning aims at splitting up an existing ontology into a set of ontology modules. Approaches for partitioning ontologies are proposed by [?,?] whereas [?] proposes a tool. On the other hand, ontology module extraction, which is also called segmentation [?] or traversal view extraction [?], aims at reducing an ontology to its relevant sub-parts. Approaches for ontology module extraction are the subject of [?,?], and the PROMPT tool [?]. More details of this category of approaches are discussed in [?].

Criteria for Ontology Modularization

Criteria for modularizing ontologies generally aim at characterizing modular ontologies. To the best of our knowledge, only [?] explicitly deals with criteria for ontology modularization. Therefore, [?] distinguishes between criteria originating in software engineering, logical criteria, local criteria, structural criteria, quality of modules, and relations between modules. First, criteria from software engineering comprise encapsulation and coherence whereas logical criteria include local correctness and local completeness. Structural criteria, which are also discussed by [?], focus on size and intra-module coherence. It is proposed to determine the quality of modules in terms of module cohesion, richness of the representation, and domain coverage. At least, to assess the relation between modules the criteria of connectedness, redundancy, and inter-module distance can be applied. Against this background, the evaluation of ontology modularization respectively applies a subset of the proposed set of criteria with respect to different scenarios and ontology modularization techniques. Based on best practices in Ontology Engineering, ontology design patterns (ODPs) simplify ontology design by providing a "modelling solution to solve recurrent ontology design problems" [?]. Several types of ODPs has already been identified, e.g., logical patterns that are used to solve problems of expressivity, or naming patterns that are conventions for naming elements. Among these types, architectural ontology design patterns (AODPs) aim at describing the overall shape of the ontology. More precisely, external AODPs describe the modular architecture of an ontology by considering a modular ontology as an ontology network. Involved ontologies are considered as modules and are connected by the import operation. A semantic web portal3 has been proposed as a repository for ODPs. Unfortunately, to our knowledge, no work has been done on proposing and describing external AODPs.

In order to characterize reoccurring structures in modularly organized ontologies, the following approach establishes a methodological basis to guide the research program of this work. This approach comprises six subsequent steps:

1. Search step: the goal of the first step is to gather modularly organized ontologies. We use the Semantic Web gateway Watson4 to search for available modular ontologies from the WWW. The search query focused on importrelationships between ontologies covering the same domain. The result is a set of 77 modularly organized ontologies. 2. Cleaning up step: the second step aims at cleaning up the initial search results in order to establish a thorough basis for further experiments. This is necessary because the set of 77 modular ontologies is afflicted with redundancies and incompleteness. This results in a set of 38 modular organized ontologies constituting a thorough basis for characterizing ontology modularization. 3. Selection of metrics step: the third step selects a set of appropriate metrics to characterize modularly organized ontologies. The modular ontologies could be described by various indicators such as the distribution of classes, the network of links between modules, the number of internal links in modules, etc. In general, the literature provides a plethora of various metrics, which could be applied for characterising modular ontologies. As a starting point, this work focuses on metrics originating in the area of software engineering, due to its maturity. In particular, this work adopts the following metrics from software engineering, which are easier to compute, in order to characterize modular ontologies [?]: (i) size of the module: the number of classes and properties (object and datatype properties), (ii) cohesion of the module: this metric is a value which is between 0 and 1 and is specified as follows: * Hierarchical Class Cohesion (HCC): the number of direct and indirect hierarchical class links. The coupling measure is computed as follows:

Coupling = α * HCD+β * HRD+δ * OP C+γ * AD α+β+δ+γ where, in our case, α = β = δ = γ = 1 4. Metrics implementation step: the fourth step implements the selected metrics. The computation was performed by the OWL API5 and the reasoner HermiT6 . This step sets up the (technical) evaluation framework. 5. Analysis step: the fifth step is the analysis of the basic population of modularly organized ontologies. 6. Result step: the sixth step involves synthesis and discussion of the results from the analysis in order to characterize modular ontologies.

Results and Discussion

A set of four patterns, which characterize M odular Ontologies MO, are proposed using previous metrics (size, cohesion and coupling). This section presents and discusses the results and the characteristics of each kind of pattern.

Pattern type 1: 1 module importing n modules

Pattern 1 contains one module which imports n other modules. For instance (Figure 1), the module WildNET.owl imports several modules such as Animal.owl, AnimalSighting.owl, BirdObservers.owl, Birds.owl, etc. The pattern that we propose conforms to an aggregation. This pattern establishes a relationship between a single module and a set of modules in the same ontology. This link is unidirectional. Applying the size metric (Table 1), the first part of the ontol- and the second part (N modules) is not structured and is respectively bigger in size. Applying the cohesion and coupling metrics (see Table 1), Pattern 1 has a high cohesion compared to the coupling metric. We consider the pattern cohesion metric to be an indicator of the degree to which the elements in the module belong together. The idea of this pattern is that the concepts grouped in an ontology should be conceptually related for a particular domain in order to achieve common goals.

Metrics Size Cohesion Coupling

M O11

WildNET

Pattern type 2: n modules importing 1 module

Pattern 2 contains n modules, which respectively import one module. For instance (Figure 2), there are three independent modules importing one module, which containes general knowledge (biopax-level1.owl ). The pattern that we propose corresponds to inheritance. This pattern establishes a correspondence between a set of modules and a single module in the same ontology. This correspondence is unidirectional. Applying the three metrics (Table 2), the first part of pattern 2 (n modules), biopax-example-ecocyc-glycolysis.owl, biopax-example-Xwnt-b-catenin.owl and Xwnt-b-catenin xref have a high coupling metric with regard to the second part of the pattern (one module) biopax-level1.owl. This means that pattern 2 is characterized by the interconnections between modules.

The degree of coupling depends on how complicated the connections are and on the type of connections between modules. As we can see, the second part of pattern 2 has a high cohesion because it encloses all other modules, which are strongly related. cohesion and coupling metrics (Table 3), the module dublincore.owl has a small size (0 concepts) with regard to other modules dcmitype.owl and terms.owl. All modules have the same degree of relatedness of concepts (cohesion) 20%. The coupling metric of the module dublincore.owl is null. In this case, pattern 3 is transformed to pattern 1 and has the same characteristics.

Metrics Size Cohesion Coupling

M O32

dcmitype.owl 21 0,20 0,02 dublincore.owl 0 0,20 0 terms.owl 47 0,20 0,25

Table 3. coupling metric with regard to the cohesion one. Two modular ontologies isometadata and iso-19115, have the same size, cohesion and coupling but they do not have a relationship like import.

Occurrence of Pattern Types

Having introduced and defined four types of patterns in order to characterize modularly organized ontologies, we consider how often these types of patterns interesting to observe that pattern type 1 accounts for about 37%. The reason for this may be the fact that this type of pattern appears to be very intuitive. It could therefore be concluded that it implicitly constitutes the rationale underlying a large part of (semi-)automatic or manual approaches for modularizing ontologies. Similarly, pattern type 4 also accounts for about 37% of the basic population, i.e. 14 modularly organized ontologies. Pattern type 4 combines Pattern 1, Pattern 2, and Pattern 3. On the one hand, it is obvious that not all modularly organized ontologies have a rather straightforward structure, which could be easily characterized. This is especially true when assuming (semi-)automatic or manual modularizing approaches, which do not use clear and precisely defined criteria. And even when these criteria are clearly and precisely defined, the modularly organized ontologies could also have such a structure depending on the overall purpose of modularization. On the other hand, it is interesting to see that even more complex structures can (almost completely) be characterized by more simple and straightforward structural forms. Moreover, pattern type 2 and pattern type 3 equally account for about 13%, i.e. 5 modularly organized ontologies. This is particularly interesting because pattern type 2 is reasonable. This is due to the fact that it appears obvious that there exists an ontology that is of significance for several other ontologies. On the contrary, pattern type 3 is much less reasonable than pattern type 2. It is really hard to understand why ontology modules respectively import each other.

In this context, it can be observed that domain ontologies combine a clear structure and organization. This means that modularization of domain ontologies tend to rely on pattern type 1 or pattern type 2. In contrast, it appears that top-level ontologies (which represent relevant knowledge to a particular domain such as medical domain) have a less straightforward structure and organization particularly when compared to domain ontologies (which represent upper (generic) ontologies, covered the knowledge of many domain types such as Biomedical ontology, Dolce). An example for this is dublincore.owl, terms.owl and dcmitype.owl, which can be characterized by pattern type 3 (Figure 3).

Conclusion and Future Work

This work aims at characterizing modularly organized ontologies to contribute to a better understanding of ontology modularization. We introduced the notion of modular ontologies, reported on approaches for ontology modularization, and reviewed existing efforts to characterize modular ontologies. To characterize modular ontologies, we followed an approach comprising six consecutive steps. This approach mainly includes the extraction and selection of modular ontologies, the selection of a set of metrics from software engineering to analyse modular ontologies, and the evaluation of the analysis results. The evaluation results in a set of four patterns, which allow for characterizing the modular organization of ontologies. These patterns show amongst other things that modularly organized domain ontologies have a clear structure whereas top-level ontologies tend to have a rather confusing modular organisation.

In the future work, we aim at using firstly further Semantic Web gateways such as Falcons or Swoogle to identify and extract additional ontologies to gain a larger basic population. Second, extending the set of metrics and applying them to the ontologies should provide further insights to modular ontologies. Third, it would be interesting to create a comparison framework to conduct experiments with different modularization approaches, comparing them to each other or to manually modularized ontologies.

  HCC = 2 * (N dHC+N idHC) N C 2 -N C where: N dHC: Number of direct Hierarchical relationships between Classes, N idHC: Number of indirect Hierarchical relationships between Classes, and N C: Number of Classes. * Role Cohesion (RC): the number of direct and indirect hierarchical role links. RC = 2 * (N dR+N idR) N Roles 2 -N Roles where: N dR: Number of direct roles between Classes, N idR: Number of indirect roles between Classes, and N Roles: Number of Roles. * Object Property Cohesion (OPC): the number of classes which have been associated through the particular object property (domain and range). ri) * N rC(ri) N Roles * (N c 2 -N C) where: N dC(r i ): Number of ontology Classes in the domain of the role r i , N rC(r i ): Number of ontology Classes in the range of the role r i , N C: Number of Classes, and N Roles: Number of Roles. The cohesion measure is computed as follows: Cohesion = α * HCC+β * RC+δ * OP C α+β+δ where: α, β and δ specify the impact of each type of hierarchical class, role or object property cohesion. In our case, we choose α = β = δ = 1. (iii) coupling of the module: it takes an estimation of the inter-dependency of different modules and is specified as follows: * Hierarchical class dependency (HCD): the number of all direct and indirect hierarchical class relationships to foreign ontologies. HCD = 1 2 * ( N edHC N dHC + N eidHC N idHC ) where: N edHC: Number of direct Hierarchical class dependencies between local classes and external classes, and N eidHC: Number of indirect Hierarchical class dependencies between local classes and external classes. * Hierarchical role dependency (HRD): the number of all direct and indirect hierarchical role relationships to foreign ontologies. HRD = 1 2 * ( N dHR N dHR + N eidHR N idHR ) where N dHR: Number of direct roles dependencies between local classes and external classes, and N eidHR: Number of indirect roles dependencies between local classes and external classes. * Object property dependency (OPD): the number of roles that associate external classes to local ones. OP D = N eRoles N Roles where: N eRoles: Number of all roles that have an external class in their domain or range, N Roles: Number of all existing roles in the ontology. * Axiom dependency (AD) : a role or a class is associated to an external ontological element through an inclusion axiom. LS(axm): the size of the left sides of the axiom axm, RS(axm): the size of the right sides of the axiom axm, LSE(axm): the number of external elements in the left sides of the axiom axm, RSE(axm): the number of external elements in the right sides of the axiom axm, and externalAssociationN umber(axm): the number of all external ontological elements that have been associated through the axiom axm to internal elements. externalAssociationN umber(axm) = LSE(axm i ) * RS(axm i ) + LS(axm i ) * RSE(axm i ) -LSE(axm i ) * RSE(axm i ).

Fig. 1 .

 1 Fig. 1. 1 module importing n modules

Fig

  Fig. 2. n modules importing 1 module

  Results of Pattern n to n-1 4.4 Pattern type 4: Pattern mix Pattern 4 combines all previous patterns (Patterns 1, 2 and 3). For instance (Figure 4), we find partterns 1 (5 * Pattern 1) and 2 (3 * Pattern 2). The proposed pattern is pattern mix. The correspondence can be undirectional and bidirectional. The major characteristic of this type of pattern is the highest

Fig. 4 .

 4 Fig. 4. Combination of patterns 1, 2, and 3

Fig. 5 .

 5 Fig. 5. Occurrence of Pattern Types

Table 1 .

 1 Results of Pattern 1 to n

	.owl	0	0,25	0
	SnakeSightings.owl 18331	1	0,25
	Snakes.owl	95	0,5	0,25
	Sites.owl	8	1	0,25
	Observers.owl	7	0,53	0,25
	Geography.owl	16	0,41	0,25
	Combined.owl	3497	0,67	0,12
	ClimateSensors.owl 255	0,61	0,6
	ClimateReadings.owl 63	0,5	0,25
	Climate.owl	24	0,61	0,15
	BirdSites.owl	437	1	0,25
	BirdSightings.owl 10401	1	0,5
	Birds.owl	1745	0,5	0,25
	BirdObservers.owl 303	0,5	0,25
	AnimalSighting.owl 26	0,66	0,25
	Animal.owl	9	0,82	0,11

Table 4 .

 4 Results of Pattern mix respectively occur. Figure 5 provides an overview of the occurrence of the different types of patterns in a population of 38 modularly organized ontologies. It is

		Metrics	Size Cohesion Coupling
		iso-metadata 2214 0,02	0,15
	M O41	iso-19108 159 iso-19103 224	0,08 0,15	0,16 0,16
		iso-19115 2214 0,02	0,15

http://www.mindswap.org/2004/

http://neon-toolkit.org/wiki/Main_Page

http://ontologydesignpatterns.org

http://watson.open.ac.uk/

http://owlapi.sourceforge.net/

http://hermit-reasoner.com/

Acknowledgment

The authors would like to thank the organizers of Summer School on ontology engineering and the Semantic Web 2001 (SSSW'2011). This work was partly realized as part of the Quaero Program, funded by OSEO, French State agency for innovation.


J C Sibel

S Reynal
email: reynal@ensea.fr

D Declercq
email: declercq@ensea.fr

Evidence of chaos in the Belief Propagation for LDPC codes

Keywords: LDPC, iterative map, chaos, Lyapunov exponent, bifurcation diagram

In this paper, we investigate the behaviors of the Belief Propagation algorithm considered as a dynamic system. In the context of LDPC (Low Density Parity-Check) codes, we use the noise power of the transmission channel as a potentiometer to evaluate the different motions that the BP can follow. The computations of dynamic quantifiers as the bifurcation diagram, the Lyapunov exponent and the reconstructed trajectory enable to bring out four main behaviors.

In addition, we propose a novel measure that is the hyperspheres method, which provides the knowledge of the time evolution of the attractor size. The information collected from these different quantifiers helps to better understand the BP evolution and to focus on the noise power values for which the BP suffers from chaos.

Introduction

The channel coding is a research field whose purpose is to protect an information to transmit from environmental disturbances. The first step is the encoding of the information, a procedure in which the information, modeled as a sequence of k bits u 1 , . . . , u k , is mapped to a larger sequence of N bits x 1 , . . . , x N . The map consists in artificial correlations called constraints or parity-check equations. In [START_REF] Gallager | Low-Density Parity-Check Codes[END_REF] are introduced the Low-Density Parity-Check (LDPC) codes which are a widespread technique to encode the information. Such a code can be represented by a Tanner graph [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF], a graphical representation which turns out to be very useful in the second step, the decoding. In this part, the bits transmitted though a random noisy channel are iteratively handled by a decoding algorithm to create an associated output sequence of N bits that verify the whole set of parity-check equations and that must be as close as possible to the input sequence. One of the most famous decoding algorithm is the Belief Propagation (BP) [START_REF] Pearl | Probablistic Reasoning in Intelligent Systems: Networks of Plausible Inference[END_REF], also used in statistical physics [START_REF] Bethe | Statistical Theory of Superlattices[END_REF], extensively studied in [START_REF] Richardson | Error floors of LDPC codes[END_REF][START_REF] Chung | Analysis of Sum-Product Decoding of Low-Density Parity-Check Codes Using a Gaussian Approximation[END_REF], which is deemed to be the optimal message-passing algorithm in the case the Tanner graph of the LDPC code is loopfree. However, in most cases the Tanner graph is not loopfree [START_REF] Murphy | Loopy belef propagation for approximate inference: an empirical study[END_REF] that involves that the BP becomes suboptimal. Moreover, the BP presents some complex behaviors in terms of the noise power of the transmitted channel, as periodic and chaotic motions [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF]. Along the whole paper, we present some measures to bring out these different behaviors. The paper is organized as follows: in the second section are presented preliminaries about the LDPC codes and the BP, in the third section we present the dynamic environment of the BP, the measures to identify the behaviors and the associated results.

LDPC codes and Belief Propagation

We consider a set of N binary random variables X = {X 1 , . . . , X N } whose global state is denoted by x = [x 1 , . . . , x N ]. An LDPC code is built by a set of M constraints, or parity-check equations, C = {C 1 , . . . , C M } such that for each check C j its state is c j = Xi∈Nj x i where N j is a subset of X called the neighborhood of C j depending on the LDPC code, and the sum is computed over the Galois field GF [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF]. We consider that C j and X i such that X i ∈ N j , form an edge e ij between two nodes inside an undirected bipartite graph called the Tanner graph G = (X ∪ C, {e ij }). An example of a Tanner graph is displayed on the figure 1. The purpose of the BP is to estimate the joint probability distribution p(X = x), written simply p(x), by a distribution called the belief b(x) and to extract the most likely state x = arg max x b(x) considered as the estimate of the input sequence. Each variable X i has its own marginal probability distribution p i (x i ) and its own belief {b i (x i )} 1≤i≤N . The BP is an iterative algorithm that consists in passing messages between the variables and the constraints on the edges of a given Tanner graph at each iteration k. For each edge e ij :

X0 X1 X2 X3 X4 X5 X6 C0 C1 C2
• the messages n

(k) ji (x i ) from C j to X i are: n (k) ji (x i ) = f ji ({m (k-1) xy } (x,y) )
• the messages m

(k) ij (x i ) from C i to C j are: m (k) ij (x i ) = g ij ({n (k) yx } (x,y) , l i (x i ))
where the functions f ij and g ji are detailed in [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF], and {l i (x i )} i are the likelihoods L(x) = {l 1 (x 1 ) = p(y 1 |x 1 ), . . . , l N (x N ) = p(y N |x N )} computed from the channel obervations {y 1 , . . . , y N }. The noisy transmission channel that we use in all of our simulations is an additive white Gaussian noise channel whose power is σ 2 . We summarize the message-passing by a unic implicit equation ∀e ij , n

(k) ji = G ij ({n (k-1) nm } (m,n) , {l i } i ). The convergence of the BP is: ∀e ij , n (k) ji = n (k-1) ji
. These messages provide the computation of the beliefs, whose equations are defined in [START_REF] Kschischang | Factor Graphs and the Sum-Product Algorithm[END_REF]. All the simulations that follow are done on a commonly used example which is the Tanner code [START_REF] Tanner | A Class of Group-Structured LDPC Codes[END_REF].

Dynamics

State space definition and properties

We consider {G ij } (i,j) as a set of iterated maps on the state variables {n

(k) ji } (i,j) and T (k) = {n (k)
ji } (i,j) as the points of the trajectories of the BP at iteration k in the associated state space E.

Parameters and scaling

In [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF] the value of σ 2 , or the corresponding Signal to Noise Ratio (SNR), is used as a parameter such that different values imply different motions of the BP. However, most of their simulations are done for particular noise realizations and scaled on the SNR, that prevents from evaluating a statistical behavior. A reason is that the noise realizations that lead the BP not to converge or to converge to a wrong estimate are rare events, essentially because the LDPC codes and the iterative algorithms are created to this end. A way to have some statistical evaluations of the behavior of the BP is:

1. finding some of these noise realizations, 2. storing the corresponding initializations on the state variables, 3. averaging the quantities to measure for a sufficient set of initializations that are close in the state space in the sense of the Euclidean distance.

By this way, we target the critical values, i.e. the SNR that correspond to a blatant change in the behavior of the algorithm, which are the bifurcations.

Bifurcation diagram

A relevant method to extract the critical values of the SNR is the use of the bifurcation diagram. It consists in evaluating the value of a function E that is computed from the state variables at their steady state for J different values of the SNR. We get a sequence [E σ1 , . . . , E σ J ] that represents the behavior of the dynamic system in terms of the SNR. We consider the following function exposed in [START_REF] Zheng | Study of bifurcation behavior of LDPC decoders[END_REF] called the mean square beliefs:

∀σ, E σ = 1 N N -1 i=0 b 2 i (0) (1)
where the input sequence in the channel is the null sequence and the beliefs are computed at the last iteration K of the BP algorithm. Obviously, there is no reason that the associated dynamic system has reached any steady state at this iteration but we need to suppose it for computation time's sake. This function presents three important values:

• E σ = 1: all the beliefs indicate that the ouput sequence is the null sequence which is a good decoding,

• E σ = 1
4 : all the beliefs are uniform thus there is no information about the state of the transmitted bits, which is a missed decoding, • E σ = 0: all the beliefs indicate that the output sequence is the complementary of the sent sequence, which is a completely wrong decoding. The display of [E σ1 , . . . , E σ J ] enables to know two properties of the used algorithm: the amplitudes provides information about the decoding performance, and the variation between successive values gives us the critical values of the SNR. We display on the figure 2 the mean bifurcation diagrams of the BP. We observe that for SNRs lower than 2.19 dB, the BP follows a regular increasing steady motion. When the SNR is greater than this critical value, the algorithm follows a periodic motion. However, for values in [2.5 dB, 2.98 dB] it does not appear any known evolution which can be an indication of chaos. For values greater than 2.99 dB, the BP converges to the good decoding state.

Reduced trajectory

Another use of the mean square beliefs function is the representation of the trajectory in a 3-dimensional state space. To this end, we use the phase space reconstruction [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF]. The method is first to compute this function at each iteration to get the following sequence E σ = [E σ (k)] 0≤k≤K . Afterthat we share this one dimensional sequence in a three dimensional sequence Ẽσ . On the figures 3 and 5 are displayed some reduced trajectories of the BP for typical values of the SNR deduced from the previous bifurcation diagram. It appears four typical behaviors that match with the four intervals exposed in the previous paragraph. We obtain a very small sized attractor for SNR = 2.10 dB that can be considered as a fixed point, whereas the reduced trajectory transforms to a limit cycle when the SNR is between 2.19 dB and 2.49 dB. A crucial point is that the thickness of the trajectory along this limit cycle increases as the SNR is getting greater up to 2.50 dB. At the same time this limit cycle interleaves with other limit cycles, that can be understood as a sequence of period doubling bifurcations in terms of dynamic system, as is displayed on the figure 4 with two interleft cycles. Such a phenomenon is a typical route to chaos, that is observable from 2.51 dB. A chaotic motion means that there is not any periodic motion or fixed point convergence anymore, as it is displayed for 2.70 dB. When the SNR reaches 2.99 dB the trajectory collapses to a single point that is a true fixed point. We have to be cautious because Ẽσ is not a true trajectory, it does not respect the Cauchy-Lipschitz condition [START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF] due to the non bijection between the messages and the beliefs. Thus, this sequence only has the role of giving clues about the true behavior as the possible shape of the actual trajectory in E that are: convergence to a fixed point, convergence to a limit cycle, convergence to a chaotic attractor. To distinguish these shapes, we need a criterion that reflects the behavior by its own value. A common measure is the Lyapunov exponent λ, that consists in evaluating at each iteration k ≤ K the log-ratio of the Euclidean distance between two initially close trajectories. As detailed in [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF][START_REF] Hilborn | Chaos and Nonlinear Dynamics: An Introduction for Scientists and Engineers[END_REF] the sign of λ reveals the behavior of the system around the corresponding initialization of the trajectories: λ ≥ 0 means the trajectories have moved away one from the other, which is an evidence of a chaotic behavior,λ ≤ 0 means the trajectories have got closer, which is an evidence of a convergent behavior to a small sized volume of the state space. This volume is reduced to a fixed point if and only if λ → ∞. When λ crosses the x-axis the system suffers from a bifurcation meaning that the algorithm has changed of motion. The corresponding SNR are the critical values. We display on the figure 6 the Lyapunov exponents of the BP on the Tanner code. As we have observed on the bifurcation diagrams, the evolution is really different as soon as the SNR is greater than 2.09 dB. The BP curves is perfectly consistent with the associated bifurcation diagram in the sense that the critical values we extract are the same and the behaviors we could imagine by the bifurcation diagram are also revealed by the Lyapunov exponent. A relevant analyze we need to effect is the comparison with the reduced trajectory we exposed previously so as to associate accurately with each reduced motion a particular evolution of λ. For SNR ∈ [0dB; 2.19db]: the reduced trajectory converges to a very small sized volume of E BP that we can consider as a fixed point whereas λ is close to the null value, for SNR ∈ [2.20db; 2.49db]: the reduced trajectory is trapped into a limit cycle whereas λ has gone over a stair, for SNR ∈ [2.50db; 2.98db]: the reduced trajectory does not converge to any fixed point, limit cycle or quasi-limit cycle but to a chaotic attractor whereas λ soars to high values, for SNR ∈ [2.99db; +∞db]: the reduced trajectory converges to a fixed point corresponding to a good decoding.

Lyapunov exponents

Hyperspheres method

We propose here a novel method to evaluate the unstability of the BP, based on its own trajectory in E. This method is complementary to the Lyapunov exponent measure because it reveals the size of the attractor that the trajectory falls into and some other properties about the limit cycles. This method consists in computing the rays R k of the hyperspheres circumscribed to the trajectory inside a given temporal window centered around each point T (k) . On the figure 7 are displayed the evolutions of two rays that correspond to two initially close trajectories in the Euclidean sense. The motion we observe for SNR = 2.30 dB is consistent with the limit cycle we observed on the reduced trajectory. The curve of the ray enables to estimate the period of the trajectory around 23 iterations. Moreover we can assert that this limit cycle is stable because the two rays cannot be distinguished. For 2.70 dB the rays moved away one from the other as it was predicted by the Lyapunov exponent observations. More accurately we can see that the rays has different oscillations step. This is due to the period doubling bifurcations explained previously. During 92 iterations in average for k ≤ 500, the trajectory is trapped in a given limit cycle and for the next 92 iterations the trajectory falls into another limit cycle of different ray. For both it is possible to measure the period or pseudo-period that is the same as the period of the first limit cycle, that is 23 iterations. For k ≥ 500 we cannot distinguish these different phases of evolution, the period doubling has led to chaos. Such an observation makes our method relevant to bring out crucial information by a one dimensional function. Another important aspect of the hyperspheres method is the raising of the behaviors difference between two initially close trajectories: we easily observe that the evolution of the rays cannot be distinguished while k ≤ 200 but as k is getting greater, the evo-lution of the rays move away one from the other but they follow the same kind of motion. For both of them the hypervolumes of the state space in which they are locking in are quite of the same size. When the SNR reaches the last critical value we observe that one of the rays decreases to the null value because the BP has converged to a fixed point. The other ray has not collapsed yet because the SNR is just at the critical value, if it was increased a little we would see the two rays going to zero.

Conclusion

In this paper, we raised the dynamics issue of the BP by the use of known and new measures. The most important result is that the BP follows a systematic scheme when the decoding is not trivial: convergence to a smallsized attractor, locking in a limit cycle, chaos and convergence to a fixedpoint. Such a property is really useful because it helps to bring out the critical values of the noise power of the channel for which the BP could present complex behaviors. Another advantage of our study is that all the measures we used can be applied on any other decoding algorithms. Therefore, we have started to create a toolbox for the dynamics study of the LDPC decoders.

Fig. 1 .

 1 Fig. 1. Tanner graph of the Hamming code

Fig. 2 .

 2 Fig. 2. Bifurcation diagrams of the BP

Fig. 4 .

 4 Fig. 4. Reduced trajectory for the BP on the Tanner code with SNR = 2.40 dB

Fig. 5 .

 5 Fig. 5. Reduced trajectory for the BP on the Tanner code with SNR = 2.70 dB and 3.00 dB

Fig. 7 .

 7 Fig. 7. Evolution of two hyperspheres rays corresponding to two initially close trajectories of the BP on the Tanner code at SNR = 2.10 dB, SNR = 2.30 dB

Fig. 8 .

 8 Fig. 8. Evolution of two hyperspheres rays corresponding to two initially close trajectories of the BP on the Tanner code at SNR = 2.70 dB, SNR = 3.00 dB