Nof Abuzainab
email: nof.abuzainab@inria.fr

Sai Rakshit Vinnakota
email: sai-rakshit.vinnakota@inria.fr

Corinne Touati
email: corinne.touati@inria.fr

Coalition Formation Game for Cooperative Cognitive Radio Using Gibbs Sampling

This paper considers a cognitive radio network in which each secondary user selects a primary user to assist in order to get a chance of accessing the primary user channel. Thus, each group of secondary users assisting the same primary user forms a coaltion. Within each coalition, sequential relaying is employed, and a relay ordering algorithm is used to make use of the relays in an efficient manner. It is required then to find the optimal sets of secondary users assisting each primary user such that the sum of their rates is maximized. The problem is formulated as a coalition formation game, and a Gibbs Sampling based algorithm is used to find the optimal coalition structure.

I. INTRODUCTION

Towards a more efficient utilization of the wireless spectrum, there has been recent interest in employing cooperation among secondary nodes in cognitive radio networks, and modeling it using game theoretic tools [START_REF] Saad | Coalitional Games for Distributed Collaborative Spectrum Sensing in Cognitive Radio Networks[END_REF]- [START_REF] Wang | Pricing-based spectrum leasing in cognitive radio networks[END_REF]. Cooperation in cognitive radio was initially considered from a collaborative sensing and access perspective in [START_REF] Saad | Coalitional Games for Distributed Collaborative Spectrum Sensing in Cognitive Radio Networks[END_REF]- [START_REF] Wang | Evolutionary Game Framework for Behavior Dynamics in Cooperative Spectrum Sensing[END_REF], where secondary users were considered transparent to the primary users. The secondary users actively listen to the primary users' channels and opportunistically transmit when the primary users' channels are idle. More specifically in [START_REF] Saad | Coalitional Games for Distributed Collaborative Spectrum Sensing in Cognitive Radio Networks[END_REF], the secondary users collaborate to improve their sensing information. In [START_REF] Saad | Coalitional Games in Partition Form for Joint Spectrum Sensing and Access in Cognitive Radio Networks[END_REF], the tradeoff between channel sensing and channel access time is captured, and the secondary users form coalitions within which they share their channel sensing information in order to improve their channel sensing and access time. Both problems in [START_REF] Saad | Coalitional Games for Distributed Collaborative Spectrum Sensing in Cognitive Radio Networks[END_REF] and [START_REF] Saad | Coalitional Games in Partition Form for Joint Spectrum Sensing and Access in Cognitive Radio Networks[END_REF] are formulated as coalition formation games, and a distributed algorithm for coalition formation that adapts based on topology and environment changes is presented. In [START_REF] Wang | Evolutionary Game Framework for Behavior Dynamics in Cooperative Spectrum Sensing[END_REF], the problem of collaborative spectrum sensing is formulated as an evolutionary game, where each secondary user can choose whether to disclose its sensing information to others or not, and the objective is to study the behaviour of selfish secondary users who make use of other secondary users sensing information in order to maximize their channel access time.

More recent problems in cooperative cognitive radio assume that the primary users are aware of secondary users, and that secondary users cooperate with primary users so as to reduce the primary user transmission time and increase the channel idle time in order to get a better chance in accessing their channels. A Stackelberg formulation is utilized in [START_REF] Afghah | A reputationbased Stackelberg game approach for spectrum sharing with cognitive cooperation[END_REF]- [START_REF] Hao | A Stackelberg game for cooperative transmission and random access in cognitive radio networks[END_REF] where the primary users are considered as leaders and the secondary users are considered as followers. The secondary users assist the primary user in delivering its information in order to get a chance of accessing the primary user's channel. Hence, the transmission time duration is divided to the primary user transmission phase, the cooperation phase, and the secondary user transmission phase, and the objective is to find the optimal time duration for each phase. In [START_REF] Afghah | A reputationbased Stackelberg game approach for spectrum sharing with cognitive cooperation[END_REF], the case of a single primary user and a single secondary user is considered, and a reputation based model is used in which every round of the game, the different time durations are updated based on the behaviour of users from the previous rounds, thus encouraging cooperation. In [START_REF] Wang | Pricing-based spectrum leasing in cognitive radio networks[END_REF] and [START_REF] Hao | A Stackelberg game for cooperative transmission and random access in cognitive radio networks[END_REF], the case of multiple secondary users assisting one primary user is considered. In [START_REF] Wang | Pricing-based spectrum leasing in cognitive radio networks[END_REF], it is assumed that all the secondary users are willing to help the primary user, and that they access the channel using time division for their own transmissions. Hence, each secondary user tries to find its optimal channel access time to maximize its utility. Also, an algorithm that selects which secondary users to assist the primary user is provided. In [START_REF] Hao | A Stackelberg game for cooperative transmission and random access in cognitive radio networks[END_REF], it is assumed that the primary user selects only one of the secondary users as a relay, and random access (Slotted Aloha) is used as a multiple access scheme among the secondary users for their own transmissions.

In this problem, we consider a similar cooperative cognitive radio model as in [START_REF] Afghah | A reputationbased Stackelberg game approach for spectrum sharing with cognitive cooperation[END_REF]- [START_REF] Hao | A Stackelberg game for cooperative transmission and random access in cognitive radio networks[END_REF]. However, our contributions are summarized as follows: 1. We extend the models in [START_REF] Afghah | A reputationbased Stackelberg game approach for spectrum sharing with cognitive cooperation[END_REF]- [START_REF] Hao | A Stackelberg game for cooperative transmission and random access in cognitive radio networks[END_REF] to the case of multiple primary users and multiple secondary users. Thus, each secondary user has to choose which primary user it should assist in order to maximize its utility, where the utility is defined as a function of the achieved throughput and the energy spent during the cooperation phase. Thus, secondary users serving the same primary user collaborate to maximize their utilities. 2. In [START_REF] Wang | Pricing-based spectrum leasing in cognitive radio networks[END_REF], it is assumed that all active relays should decode the information before transmitting in the cooperation phase, and thus the achieved rate will be severely degraded by the relays with low channel quality. Thus in this problem, we assume that the secondary users employ sequential relaying to improve the received rate by the primary user, and hence increase their channel access duration. 3. We formulate the problem as a coalition formation game, and we draw a connection between the coalition formation game and potential games [START_REF] Monderer | Potential Games[END_REF]. To our knowledge, this has not been previously addressed.

We propose a distributed algorithm based on Annealed

Gibbs Sampling to find the optimal coalition structure.

II. SYSTEM MODEL

We consider the downlink scenario in a Cognitive Radio (CR) network composed of a set P of primary users (PU) and a set S of secondary users (SU). It is assumed that each primary user p (p ∈ P) has a fixed rate requirement of R p bits per seconds over the time interval [0, T ] where T is in seconds, and that the rate requirement is less than or equal to the link capacity between the base station and the primary user. Also, it is assumed that the primary users are assigned orthogonal frequency channels. The channels between any pair of nodes is modeled as slow Rayleigh fading i.e. the value of the fading coefficient does not change over the interval [0, T ], and that all channels are independent. Additive white Gaussian noise of variance N 0 is assumed to be present at each of the users.

In order to give the secondary users the chance to access the primary users' channels, the base station allows the secondary users to assist it in satisfying the demands of the primary users. Hence, cooperation is employed, and during the cooperation phase, each secondary user s (s ∈ S) chooses which primary user it will assist to satisfy its demand. Hence, for each primary user p, we define coalition P made of the set C P of secondary users serving p. Also, the time interval [0, T ] is divided into two main phases:

1) The cooperation phase: During the first fraction α P of T , the secondary users in set C P assist primary user p in satisfying its rate demand.

2) The SUs' transmission phase: During the time fraction 1-α P of T , the secondary users in set C P will share the channel of primary user p for their own transmissions.

A. The cooperation phase

To serve each primary user p during the cooperation phase, we assume that sequential relaying is used, and we employ the same transmission scheme as in [START_REF] Beres | Optimal Relay-Subset Selection and Time-Allocation in Decode-and-Forward Cooperative Networks[END_REF]. In this scheme, it is assumed that the secondary users in set C P have a certain order in the coalition such that σ P (k) is the k th secondary user in coalition set C P . Hence, the cooperation phase is divided into several time fractions t P k (0 ≤ k ≤ C P = |C P |), and transmission occurs as follows:

• During the first time fraction t P 0 , the base station transmits to the first secondary user σ P (1). Then, σ P (1) decodes while the remaining SUs and the PU store the received information.

• During any subsequent t P k (1 ≤ k ≤ C P -1), secondary user σ P (k) transmits to secondary user σ P (k + 1) which decodes while the remaining SUs σ P (k + 2),...,σ P (C P ) and primary user p store the newly received information.

• During the last fraction t P C P , secondary user σ P (C P ) transmits to PU p, which ultimately decodes the message. A coalition of SU P is therefore a triplet (C P , σ P , t P ), where t P denotes the vector of time fractions t P = (t P 0 , . . . , t P C P ) and is such that t P > 0 and t P .1 = 1, with 1 refering to the vector of size C P + 1 whose elements are all equal to 1 and t P .1 refering to the inner product of vectors. Notation t P > 0 means that the inequality holds componentwise.

In the remainder of this section, we focus on a single coalition P and therefore drop the corresponding index for clarity of notations.

Based on this scheme, for a given coalition P, the mutual information at each secondary user σ(k), and at primary user p can be written as [START_REF] Beres | Optimal Relay-Subset Selection and Time-Allocation in Decode-and-Forward Cooperative Networks[END_REF]:

   R σ(1) (t) = t 0 L B,σ(1) R σ(k) (t) = t 0 L B,σ(k) + ... + t k-1 L σ(k-1),σ(k) R p (t) = t 0 L B,p + ... + t C L C,p , (1)
where L ,m is the channel capacity of the link between transmitter and mobile m, and B refers to the base station, and is given by shannon's capacity formula

L ,m = log 2 1 + |h ,m | 2 P R ,m N 0 , (2)
where h ,m is the value of the fading coefficient of the link of transmitter and mobile m. In both phases, it is assumed full channel state information and that the users transmit with fixed power value P . The received power P R ,m by each node follows the pathloss model. Hence, it is given by P R ,m = |d ,m | -a P where d ,m is the distance between the transmitting and receiving nodes, and a is the pathloss exponent.

Thus, the transmitted rate from the base station B to the primary user p is

α.R with R = min k {R σ(1) (t), ..., R σ(k) (t), ..., R p (t)}. (3)
The time fraction α is chosen so as to satisfy the primary user constraint:

α = R p /R. (4)
For ease of notations, we define σ(0) = B the base station and introduce the matrix of channel capacities L = (L σ(a-1),σ(b) ) a,b :

L=         L B,σ(1) L B,σ(2) L B,σ(3) • • • L B,p L σ(1),σ(2) L σ(1),σ(3) • • • L σ(1),p L σ(2),σ(3) • • • L σ(2),p 0 . . . . . . . . . . . . L C,p         . (5)
Then, equation (1) becomes:

∀k ∈ C ∪ {p}, R k (t) = (t .L) k .
And finally, R can be rewriten

R = min 1≤k≤C+1 (t .L) k . (6)

B. The SU transmission phase

Time division multiplexing is assumed, and the transmission time is divided based on the contribution of the secondary users in C P . Hence, the time allocated for secondary user k is given by (1 -α P )t P k . Thus, its reward is made proportional to the amount of energy spent by the secondary user to assist the primary user in the cooperation phase. The corresponding achieved throughput of each secondary user is:

u k (P) = (1 -α P (P))t P k L B,k . (7)
This is the utility of secondary user p when in coalition P.

III. COALITION STRUCTURING

In this section, we study the performance and optimal structure of a coalition P for a given PU p, given the set of SUs composing it C P .

Recall that a coalition is a triplet (C P , σ P , t P ). For a given C P , the maximum achievable rate of the coalition is given by α P .R * (P) with R * (P) = max

C ⊂C P max σ P ,t P k >0, t P .1=1 min 0≤k≤C P (t P L P ) k ,
with σ P being a permutation σ P :

C P → C P . (8)
Again, as we focus on the optimal relay organization for a given PU p, we omit the p index in what follows. Also, the proofs of the following propositions and theorems of this paper have been omitted due to space limitation. Please check [START_REF] Abuzainab | Coalition Formation Game for Cooperative Cognitve Radio Using Gibbs Sampling[END_REF] for the detailed proofs.

A. Optimal Time Fractions

Consider first the optimization problem given by equation [START_REF] Beres | Optimal Relay-Subset Selection and Time-Allocation in Decode-and-Forward Cooperative Networks[END_REF]. We first consider the sub-problem: σ R(P) = max

t P k ≥0, t P .1=1 min 0≤k≤C P (t P L P ) k , (9)
In other words, we fix the set of relays and the relay order σ, and discuss about σ t the maximizing time fraction vector, and σ R the corresponding solution.

Note that the relaxed optimization problem corresponds to the 2 player zero-sum game with matrix game L (which depends on σ), when the coalition relays play the role of the row player which plays against nature. Then, σ R is the value of the game and σ t is the corresponding optimal strategy. We define the set PK = {q ∈ C, σ t > 0}. We have two important properties:

Proposition 1. Assume that PK = C P . Then, σ t = L -1 1 1 L -1 1 and σ R = 1 1 L -1 1 . ( 10
)
Proposition 1 illustrates the fact that if all relays are active, the time allocated for the transmission of each relay k should be such that it passes on to the next active relay k + 1 the information that k + 1 hasn't yet received from the previous transmissions. Then, the corresponding time fractions and throughput can be computed by matrix inversions, as already noted in [START_REF] Beres | Optimal Relay-Subset Selection and Time-Allocation in Decode-and-Forward Cooperative Networks[END_REF]. Also, as the L matrix is lower triangular, the inversion process can be obtained by forward substitution.

In the subsequent propositions, we will first state the conditions for a relay to be active and then formulate the problem of finding the optimal time allocation for the general case when some relays can be inactive. Proposition 2. ∀q ∈ PK, and ∀k ∈ C, (

σ t L) k ≥ ( σ t L) q .
A corrollary is that ∀k, q ∈ PK, (

σ t L) k = ( σ t L) q .
Simply put, Proposition 2 shows that, at the optimal time vector σ t, any used relay forwards the same amount of information to the next one: there is neither loss of information nor wasted time. The fact that all unused relays k would have a greater value of ( σ t L) k means that the σ(k + 1) have already received all information without the help of relay σ(k) which is therefore unused. Theorem 3. For a given matrix L, the optimal σ t and σ R are the solution of the linear program:

max t,R R s.t.                ∀i, t i ≥ 0, C i=0 t i = 1, ∀k, R ≤ C j=0 L k,j t j . (11)

B. Optimal Relay Ordering

For a given C P and σ P , Theorem 3 gives the optimal corresponding t P . In this subsection, we focus on the joint optimization of σ P and t P for a given coalition set.

The optimal choice is crucial in terms of performance, as shown in the next proposition.

Proposition 4. For a given coalition set, the degradation of throughput due to the choice of order σ may be unbounded. In other words,

∃C P , ∀Z > 0, ∃σ 1 , σ 2 , σ2 R/ σ1 R > Z.
Since relay order significantly impact the performance of a coalition, we propose an algorithm for constructing a relay order, which is built on Proposition 5.

Proposition 5. For a certain order of relays, the optimal time fractions of relays are (strictly) positive if and only if each relay at the k th position (k = 1, 2, ..., C) satisfy:

           R σ(k) (t) > t 0 L B,r + k-1 m=1 t m L σ(m),σ(r) , ∀σ(r) > σ(k), R σ(k) (t) > t 0 L B,p + k-1 m=1 t m L σ(m),p .
Corollary 6. It follows from Proposition 5 that each optimal time fractions is positive if and only if each relay k is such that ,k .

σ(k) = argmin k t k with t k = R σ(k-1) (t) - k-2 i=1 L σ(i),k t i L σ(k-1)
Remark 1. In this problem, since we are considering Rayleigh distribution of the channel coefficients, the channel capacities have real values, and thus we assume that two relays having the same t k occurs with probability 0. Hence, there exists only one relay that minimizes t k .

The proposed relay ordering amounts in taking the next relay as the one for which the needed time to transmit the information will be minimized. Algorithm 1 returns such order σ as well as the corresponding optimal time fractions σ t. Further, the ordering will be the one where all relays have positive optimal time fractions (i.e. where all relays participate to the transmission), if such order exists (which then is unique, following Remark 1). If not, the algorithm returns a relay order and the set of unused relays.

Algorithm 1: Relay Coalition Ordering Algorithm Input: The set C and matrix L. Output: Permutation σ, time fractions t and set of unused relays V .

1 Set J = C to be the current set of relays; 2 Set the array index k = 1 to be the current relay position; 3 Initialize the variable T j = 1 for all j in J; 4 Initialize the variable T p = 1 for primary user p; 5 Initialize the variables sum = 0; mini = 1/L 1,1 6 while J is not empty or V is empty do 7 foreach secondary user j in J do 8

If (k > 0) Set T j = T j -t[k -1]L k,j ; 9 Set t j = Tj L k,j ; 10 If (t j < mini) Set mini = t j ; 11 end 12 Select secondary user s = mini; 13 Set t p = Tp L k,p ; T p = T p -t[k -1]L k-1,p ; 14 if t s < t p then 15 Increment k; 16 Set σ[k] = s; t[k] = t s ; J = J\{s}; 17 Set sum = sum + t[k]; mini = 1/L 1,1 ;

C. Optimal Coalition Set

We finally look at the problem of jointly optimizing the coalition set C P . Braess-like situations can occur: adding a relay may decrease the throughput achieved by the coalition. Indeed, a user with low channel capacities will require a higher transmission time and therefore be detrimental to the coalition.

For example, consider a system with symmetric capacity links and

L B,1 = 6, L B,2 = 8, L B,3 = 2, L 1,2 = 1, L 1,3 = 2, L 2,1 = 5, L 2,3 = 2.
The achieved throughput is 2. Then add secondary user 0 with L B,0 = 10, L 0,1 = 4, L 0,2 = 1, L 0,3 = 1. Then the total resulting achieved rate falls to 20/11 (see [START_REF] Abuzainab | Coalition Formation Game for Cooperative Cognitve Radio Using Gibbs Sampling[END_REF] for the details on the corresponding relay orders and time fractions). Hence, it is necessary to find the optimal set of relays to assist each primary user. This is why we propose in the following section an algorithm based on Gibbs Sampling.

IV. COALITION FORMATION GAME

As each secondary user is devoted to assisting a single primary user, the secondary users are partitioned into disjoint sets. Hence, each group assisting a given primary user is considered as a coalition, and the whole set of secondary users is mapped into a coalition partition. The value V of a coalition P is the sum of the utilities of each of its member for their own tranmission, as given by Equation [START_REF] Monderer | Potential Games[END_REF].

V (P) = k∈P u k (P) = (1 -α P (P)) k∈P t P k L B,k . (12)
A. Allocation Coalition Game Congestion games are games where the set of players share a set of resources and where each player takes an action by selecting which of the resources to use. The payoff of each player depends on the number of players using the same resources. Allocation games are more general games, where the payoff of each user depends on the set of players using the same resources.

Congestion games have interesting optimization properties. Indeed, the class of congestion games is known to be the class of exact potential games, which are games where there exists a function F such that the change in the utility of any player (due to a change of his strategy) can be computed as the change in the value of F due to the change of that strategy. It has been shown in [START_REF] Monderer | Potential Games[END_REF] that the local maximizers of the potential function F are the Nash equilibria of the potential game.

Allocation games are in general not potential games, but it has been shown in [START_REF] Coucheney | Selection of Efficient Pure Strategies in Allocation Games[END_REF] that a simple tranformation on the utilities can turn an allocation game into a potential game.

Note that these classes of games do not yet have, to the best of our knowledge, their counterparts in the coalition game theory framework. Note also that coalition games differ from these games in that the players cooperate within each coalition and compete with the other coalitions.

Inspired from the definition of allocation games, we introduce the coalitional allocation games as follows: Definition 8. A coalition game satisfying the 3 following properties is said to be an coalitional allocation game:

1) The number of coalitions A is given (although a coalition may be empty) 2) Coalitions are indexed by parameter a, 1 ≤ a ≤ A.

3) The value of each coalition is a function on parameter a as well as the set of members of the associated coalition C a , but does not depend upon the coalitions formed by other members.

Note that allocation coalition games are not in characteristic form in that the value of a coalition depends on the coalition index p. Proposition 9. The coalition games formed by secondary users that cooperate to assist the primary users with value function of Equation ( 12) is an coalitional allocation game.

Our game can also be thought an allocation game where the players are the secondary users, the resources are the channels of the primary users, and the payoffs are the secondary users utilities as given by Equation [START_REF] Monderer | Potential Games[END_REF].

The value of a coalition V (P) not only depends on the set of secondary users that it consists in (i.e. C P ) but also of which primary user p the set is assisting. Hence, we define the coalition structure CS to be an |P| dimensional vector in S |P| , where each entry P, (1 ≤ p ≤ |P|) is the set of SUs assisting PU p. Also, the entries should satisfy p C P = ∅.

Pursuing the analogy with allocation games, we have the following result: Proposition 10. Suppose that the advertised utility for player k when inside coalition a is the repercussion utility:

r k (C a ) = u k (C a ) - j∈ ,j =k u j (C a \{k}) -u j (C a ) . (13)
Then, the set of stable coalition partitions CS * are the maximizers of the social welfare, i.e. the sum of valuations V of the different coalitions.

W (CS * ) = max CS a V (C a ). (14)
The remaining of this section presents our solution method.

B. Finding the Optimal Coalition Partition

In order to find the optimal coalition structure, we propose a randomized algorithm based on Annealed Gibbs Sampling [START_REF] Bremaud | Markov Chains, Gibbs Fields, Monte Carlo Simulation, and Queues[END_REF]. For our algorithm, we will allow one user at a time to move to a new coalition thus forming a new coalition structure. The Annealed Gibbs Sampling based algorithm for our problem is defined in Algorithm 2. Starting from an initial coalition structure CS 0 , the algorithm first picks a secondary user j at random.

Then at each time step t and for all coalitions C P , the algorithm computes the repercussion utility r j (C P ) of secondary user j when moved to coalition C P . Finally, the algorithm computes the probability of moving j to coalition C P . The expression of the probability used is the one used in Gibbs Sampling and it is often known as the Gibbs measure: when T is a parameter, commonly known as the temperature, and is often used to control the randomness in jumping to suboptimal solutions. Then, the algorithm moves secondary user j to a coalition C P according to the computed probability distribution. This process is repeated until convergence.

In what follows, we study the convergence of the algorithm.

Proposition 11. When T = 1 log(t) , Algorithm 2 converges to a global optimal solution.

Although Gibbs Sampling can converge to the global optimal solution, there are no guarantees on its convergence time. In our problem, the algorithm is supposed to run in real time. Hence, the algorithm is not allowed to exceed a certain time duration. Thus, we set a maximum number of iterations for the algorithm. As for the location of the nodes, the coordinates of the the primary and secondary users are chosen randomly according to a uniform distribution on a 10×10 rectangular grid, whereas the base station is placed at the center of the grid. Without loss of generality, we neglect the effect of fading in the simulations by setting the value of all fading coefficients to be one. The rate demand of each primary user is set to be the channel capacity of the link between the base station and the primary user. All transmission power values are set to be 0.5 Watts.

For the Gibbs Sampling coalition formation algorithm described in Algorithm 2, we consider two cases for the value of the temperature T . The first one when the temperature is set to be T = 1/log(t) where t is the current iteration value, while the second case when the temperature is very low and set to be T = 0.001. The second case corresponds to the greedy algorithm that chooses the next best solution at each iteration. The maximum number of iterations is set to be 1500.

Then, we compute the total utility obtained at each iteration of Algorithm 2 for both cases when T = 1/log(t) and when T = 0.001. Also, we compute the maximum utility obtained using brute force optimization in order to determine the value of the global optimal solution. Figure 1a shows the the evolution of the total utility as Algorithm 2 elapses.

Based on Figure 1a, we first notice for the case when T = 1/log(t) considerable fluctuation in the value of the total utility, and that sometimes the total utility value drops in the subsequent iteration, and this follows from the random nature of Gibbs Sampling that allows jumping to less optimal solution for the purpose of escaping from local optimum solutions and eventually reaching the global optimal solution. Also, we see that the global optimum (obtained from brute force optimization and which is found to be 0.5914) is attained for the first time at iteration number 277 which is relatively fast. Due to the stochastic nature of the algorithm, the total utility keeps fluctuating but keeps closer to the global optimal solution. For the case of T = 0.001 (i.e. the greedy choice), the total utility converges fast to a suboptimal solution where the total utility value is found to be 0.1579, which is considerably lower than the global optimal solution. Figure 1b shows the optimal coalition structure and the relay order for each coalition. The secondary users and the base station are represented by black circles, where the base station is the circle at the center of the grid. Primary user 1, 2, and 3 are represented by the blue, green and red circles respectively. Each of the blue, green and red lines connects the secondary users assisting primary user 1, 2 and 3 respectively. The secondary users are connected from the base station to the primary user based on their order obtained from the relay ordering algorithm. Since fading is not considered in this case, it is clear to observe from Figure 1b that the secondary users are connected based on the relative proximity to each others, to the base station and to the primary users while the secondary users that are far away from the primary users and the base station do not assist any of the primary users. Hence, this shows the effectiveness of our Gibbs Sampling algorithm in selecting the secondary users that are mostly beneficial to each primary user thus reaching the optimal solution.

VI. CONCLUSION

We have formulated the problem of cooperation among primary users and secondary users in a cognitive radio network as a coalition formation game, and proposed a Gibbs Sampling based algorithm in order to find the optimal coalition structure. The results show that our Gibbs Sampling based algorithm can reach the global optimum value within an acceptable time duration, unlike the case of greedy algorithms which they are more likely to converge at a local optimum value. The results also show the dependence of the coalition structure on the system parameters such as the distance between the nodes and the rate demand of the primary users.

  j = 1, ..., k do 23 Set t[j] = t[j] sum ; 24 end Proposition 7 (Relay Coalition Ordering Algorithm). The algorithmic complexity of Algorithm 1 is of O(C 2 ).

Algorithm 2 :1 8 end 9

 289 Primary User Selection Using Annealed Gibbs Sampling Set u j = 0 for all secondary user j 2 Define the set C P = ∅ for all primary user p 3 Initialize temperature T 4 forall the time epoch t do 5 Pick randomly a secondary user j; 6 forall the coalition set C P do 7 Compute the value of the repercussion utility r j (C P ) of secondary user j when moved to C P Randomly pick up coalition C according to the Gibbs measure (Gibbs-Measure) 10 Move secondary user j to coalition C 11 Update temperature T (e.g. according to T = 1 log(t) ) 12 end V. NUMERICAL RESULTS In this section, we study the performance of the coalition formation algorithm. For evaluation, the following values are used: |P| = 3, |S| = 10, N 0 = -40.87 dBm, a = 3.4.

  Optimal coalition structure and relay order


Nof Abuzainab
email: nof.abuzainab@inria.fr

Corinne Touati
email: corinne.touati@inria.fr

Multihop Relaying in Millimeter Wave Networks: A Proportionally Fair Cooperative Network Formation Game

Millimeter wave channels suffer from considerable degradation in the channel quality when the signal is Non Line of Sight (NLOS) between the source and the destination. Multihop relaying is thus anticipated to improve the communication between a source and its destination. This is achieved by transmitting the signal to a sequence of relays in which a Line of Sight (LOS) signal exists between two nodes along the path, or more generally when the signal is better than the transmitted signal directly from the source to the destination. In this paper, we consider a millimeter wave network composed of multiple sourcedestination pairs and a set of deployed relays. We formulate the problem of multihop relaying as a cooperative network formation game in which each relay chooses which source-destination pair to assist in order to improve the end-to-end performance, that is, the multihop delay between the source and the destination. Further, we present an algorithm based on the Nash Bargaining Solution to ensure fairness among the different source-destination pairs and assess its efficiency on numerical simulations.

I. INTRODUCTION

The interest in Millimeter Wave Communications has been tremendously increased as a viable technology for fifth generation wireless cellular systems. This is due to the fact that millimeter wave communications support the very high data rates necessary for broadband and multimedia communications thanks to the availability of large bandwidth at the high frequencies. However, communications at these high frequencies suffer from two main drawbacks. The first is that the millimeter wave signal suffers from severe pathloss. To overcome this, there is an active research going on designing beamforming techniques in order to extend the signal range and enable communication between the targeted transmitter and receiver [START_REF] Roh | Millimeter-wave beamforming as an enabling technology for 5G cellular communications: theoretical feasibility and prototype results[END_REF]- [START_REF] Hur | Millimeter Wave Beamforming for Wireless Backhaul and Access in Small Cell Networks[END_REF]. The second drawback is that the millimeter wave signal gets severely attenuated in the case of Non Line of Sight (NLOS) [START_REF] Rangan | Millimeter Wave Cellular Networks: Potentials and Challenges[END_REF], [START_REF] Pietraski | Millimeter wave and terahertz communications: Feasibility and challenges[END_REF]. To improve communications in case of NLOS, the use of intermediate relays that have LOS (or in general a better) signal with the source, the destination or among each others is suggested. Hence constructing a path between the source and destination using those relays improves the sourcedestination communication. In this paper, we focus on the second challenge and attempt to design a multihop relaying technique for a millimeter wave networks.

Relaying in general is a well studied topic in wireless communications, and there is a vast literature covering multihop relaying. However, multihop relaying in the context of millimeter wave communications is yet a new topic, and there has been still few works that deal with this issue. The work in [START_REF] Qiao | Efficient concurrent transmission scheduling for cooperative millimeter wave systems[END_REF] selects for a given source-destination pair employing millimeter wave RF the best relay within the beamwidth of the source to assist in the transmission in case of Non Line of Sight (NLOS) between the source and the destination. Also based on this relay selection mechanism, a scheduling algorithm is presented for the case when multiple sourcedestination pairs are present. In [START_REF] Kim | Quality-Aware Millimeter-Wave Device-to-Device Multihop Routing for 5G Cellular Networks[END_REF] a centralized algorithm for multihop relaying routing that takes into account the characteristics of the millimeter wave transmissions is presented. In particular, the presence of multiple source-destination pairs, where each source is interested in video streaming to its destination is assumed. Further, the performance is measured in terms of a differentiated quality function of each flow. The algorithm then finds a feasible route of relays for each sourcedestination pair, and the objective is to maximize the sum of differentiated quality functions for all flows.

Our work considers the multihop routing problem for multiple source-destination pairs employing millimeter wave RF as in [START_REF] Kim | Quality-Aware Millimeter-Wave Device-to-Device Multihop Routing for 5G Cellular Networks[END_REF]. However, our approach is different because we formulate the multihop relaying problem as a network formation game. Network formation games have been recently used for multihop relaying in wireless networks (see [START_REF] Saad | Network Formation Games Among Relay Stations in Next Generation Wireless Networks[END_REF], [START_REF] Li | Cooperative multihop relaying via network formation games in cognitive radio networks[END_REF])

but not yet in particular for millimeter wave networks. In [START_REF] Saad | Network Formation Games Among Relay Stations in Next Generation Wireless Networks[END_REF], an algorithm based on network formation game is presented that constructs a uplink multicast tree of relays, to which the mobiles can connect to in order to communicate with the base station. In this game, the relays are the players, and their objective is to connect to the tree in such a way that maximizes their utilities, where the utility is measured in terms of per hop delay and bit error rate. In [START_REF] Li | Cooperative multihop relaying via network formation games in cognitive radio networks[END_REF], a network formation game is formulated for multihop relaying in Cognitive Radio Networks. The game has also a Stackelberg approach in which the primary source-destination pairs are considered as the leaders, and the secondary users are considered as the followers, and the objective is to construct a path of secondary users (that act as relays) between each source-destination pair so as to improve its transmission and to eventually give the secondary users chance of channel access. In both problems, the network formation games are non cooperative i.e. the player moves are based on maximizing their individual utilities.

In contrast, our approach is based on cooperative network formation games and in particular a coalitional graph game in which the path between each source destination pair is constructed in a distributed fashion i.e. each relay decides on joining the path of a certain source-destination pair. Each group of relays along the same path forms a coalition. But as opposed to all previously mentioned works, our algorithm achieves proportional fairness in order to maintain an acceptable quality for every source-destination pair and to ensure fairness in relay assignment among the different source destination pairs.

II. SYSTEM MODEL

We consider a set M of source-destination pairs {(s i , d i )}, (i = 1, 2, ..., M , M = |M|) where each source s i has a file of B i bits to deliver to destination d i , and a set N of deployed relays. It is assumed that the nodes employ millimeter wave RF. Then, the received power P R is given by

P R = AM T M R d -α P T ,
where M T and M R are the antenna gains at the transmitting and receiving nodes respectively, A and α are the pathloss coefficient and exponent, P T is the transmitted power, and d is the distance between the transmitting and receiving nodes.

Millimeter wave signals get severely attenuated with distance. Hence, we assume that all nodes employ directional beamforming, and that each pair of communicating nodes engage in a beamstearing algorithm in order to achieve the maximum directivity gains. It is further assumed that the beamwidth is very small as transmitting at very high frequencies permits very narrow beamforming. (Some current products such as [START_REF]About Millimeter Wave[END_REF] demonstrate that the beamwidth can be as small as 2 degrees and that interference can be eliminated even with nodes along the same path.) This makes it very unlikely for two pairs of nodes to interfere with each other and therefore interference is neglected. Additive white Gaussian noise with zero mean and power spectral density N 0 is assumed to be present at each node. Hence, the received Signal to Noise Ratio (SNR) is given by SNR = P r /N 0 and we assume that the achieved rate R is related to the SNR through Shanon's capacity formula i.e.

R = W log(1 + SNR) = W log 1 + AM T M R d -α P T N 0
where W is the available bandwidth. Further, millimeter wave signals can get severely attenuated with blockage, and thus the signal can get considerably degraded in the case of non line of sight (NLOS). Hence, the channel quality between any pair of nodes is dependent whether a line of sight (LOS) signal exists or not. In particular, measurements (such as in [START_REF] Rangan | Millimeter Wave Cellular Networks: Potentials and Challenges[END_REF]) have shown that different pathloss models exist for the LOS and the NLOS cases. We define A N and α N to be the pathloss coefficient and exponent respectively for the NLOS case and A L and α L to be the pathloss coefficient and exponent for the LOS case.

Since the direct channel between each source-destination may be NLOS, the achieved rate using direct transmission may be low and incur significant delay to deliver file from the source and the destination. The objective is then to devise multihop relaying i.e. to try to find a path between each sourcedestination pair using relays (as some relays might have LOS signal with the source and destination and among each other) so as to improve the communication between each sourcedestination. The performance is assessed by computing the multihop delay i.e. the time spent through the path to deliver the file from the source to its destination. In order to find a path for each source-destination pair, we design an algorithm using a cooperative network formation game or more specifically a coalition graph game in which each relay chooses to connect to one source-destination pair. Our algorithm also ensures proportional fairness among the source-destination pairs.

III. COALITION GRAPH GAME FORMULATION

A. Problem Formulation

We formulate our problem as a cooperative network formation game or more specifically a coalition graph game, where the players are the relays. Each relay chooses to assist one source-destination pair by connecting itself along the path between the chosen source and destination in a way that achieves the best performance possible. Thus, a group of relays assisting the same source-destination pair is considered as a coalition. This leads us to the following definitions: Definition 1 (Path). A path P i between source s i and destination d i is a sequence σ 0 , σ 1 ,...,σ Ni ,σ Ni+1 , where N i is the number of relays in the path, σ 0 = s i and σ Ni+1 = d i and σ 1 ,...,σ Ni are the relays along P i . In other words, it is the set of edges given by

P i = {< σ j , σ j + 1 >, 0 ≤ j ≤ N i }.
Definition 2 (Action Set). Each relay r (which is presently either unused or assisting a source-destination pair (s i , d i )) can decide to perform action a k and assist source-destination (s k , d k ) by inserting itself between two consecutive nodes along the path of (s k , d k ) in a way that achieves the minimum possible multihop delay. In other words, if P k is the current path between s k and d k , P k (r, j) is the path formed by inserting relay R between nodes j and j + 1 along the path:

P k (R, j)=(P k \{< σ j , σ j+1 >}) ∪ {< σ j , r >, < r, σ j+1 >}.
Relay r inserts itself between nodes σ j * and σ j * +1 such that j * = argmax j D k (P k (r, j)) where D k (P k (r, j)) is the multihop delay along path P k (r, j). We denote by P * k (r) the resulting path. Also, we define action a 0 where the relay decides not to assist any pair. The action set for each relay is then

A = {a k , 0 ≤ k ≤ M }.

B. Proportional Fairness Maximization

We are interested in allocating the relays to the sourcedestination pairs in a fair way so as to avoid situations in which all relays would be allocated to a few source-destination pairs (that have better channels with the relays) enjoying very enhanced performance, while other source-destination pairs would not be adequately assisted by the relays and get poor performance. One approach in cooperative game theory is the Nash Bargaining Solution (NBS) [START_REF] Nash | The Bargaining Problem[END_REF]. In the NBS, the objective is to choose the strategies of the players that maximize the following objective function, commonly known as the Nash product:

max s∈S N i=1 (U i (s) -a i ),
where N is the number of players and U i (s) the utility for player i when all players take actions represented by vector s and vector a = (a 1 , a 2 , ..., a N ) is known as the disagreement point, where each a i corresponds to the utility value of player i when no agreement is reached. More often, the value of the d i s are assumed to be zero. Using this assumption and by taking the logarithm of the objective function we get our objective function: max s∈S N i=1 log(U i (s)) which corresponds to the proportional fairness [START_REF] Kelly | Rate control for communication networks: shadow prices, proportional fairness and stability[END_REF].

Definition 3 (Coalition Value). The value (or utility) of the coalition of each source-destination pair (s i , d i ) is expressed in terms of the multihop delay, i.e. the time required to deliver the file from s i to d i , which we denote D i (P i ), with P i the path corresponding to the relays in the coalition C i . Hence, in order to minimize the proportional fair sum of D i s, we define the value of coalition C i to be V (C i ) = -log(D i (P i )). Then, maximizing the sum of the values of the coalitions amounts to maximizing the proportional fairness of the utilities of the source-destination users with their utilities being proportional to their transfer rates U i (P) = 1/D i (P i ).

We consider the multihop delay as being the sum of delays of all edges from s i to d i . Hence, it is given by the following expression

D i (P i ) = B i Ni j=0 1 R σj ,σj+1 , (1)
where N i is the number of relays along the path, B i is the size of the file to transfer, R σj ,σj+1 is the rate achieved between node σ j and node σ j+1 . Again, we assume that node σ 0 is the source s i and node σ Ni+1 is the destination d i .

Any other node σ j (1 ≤ j ≤ N i ) is the j th relay along the path. Note that in the delay expression of Equation 1, it is assumed that each relay decodes the whole file before transmitting it to the next relay along the path. The multihop delay can be improved in the case where the file is divided into packets and transmissions occur packet by packet. However, the expression is more complicated to handle and relies on assumptions on the packet based system (see, e.g. [START_REF] Saad | Network Formation Games Among Relay Stations in Next Generation Wireless Networks[END_REF] and [START_REF] Li | Cooperative multihop relaying via network formation games in cognitive radio networks[END_REF], in which packets arrive at each source at a certain rate, and the average delay at each hop is computed based on modeling the packet service system as M/G/1 queue). Hence, our expression constitutes a simple upper bound on the delay for packet based transmissions.

The following theorem shows how the optimal actions of the relays maximize the proportional fairness sum.

Theorem 1. If each relay r chooses to connect to the sourcedestination pair (s i , d i ) that maximizes its marginal contribution (i.e. that has the maximum log(D i (P i ))log(D i (P i (r))) if it is positive and to chooses not to assist any pair -i.e. to remain unused -otherwise), the corresponding equilibria are the maximizers of the proportional fair sum.

Proof. We consider the game of transferable utility in which the utility (welfare) of each relay in a coalition is proportional to the collective contribution of all relays in the coalition i.e. the coalition value. Hence, when relay r is assisting sourcedestination (s i , d i ), the coalition value is divided equally among the N i relays in the coalition:

u r (C i ) = V (Ci)
Ni . We also set the utility of any unused relay to be zero (as if they were belonging to a dummy path with null coalition value). Now, we define the repercussion utility of relay r in coalition C i as

r r (C i ∪ r) = u r (C i ∪ r) - Ni k=1,k =r (u k (C i ) -u k (C i ∪ r)).
By substituting the values of the utilities u r (C i ∪r), u k (C i ∪ r), and u k (C i ) into r r (C i ∪ r), we get:

r r (C i ) = u r (C i ∪ r) - Ni k=1,k =r (u k (C i ) -u k (C i ∪ r)) = V (Ci∪r) Ni - Ni k=1,k =r V (Ci) Ni-1 -V (Ci∪r) Ni = V (Ci∪r) Ni -(N i -1) V (Ci) Ni-1 + V (Ci∪r) Ni = V (C i ∪ r) -V (C i ) = log(D i (P i )) -log(D i (P i (r))).
Recall that we assume that the utility of each unused relay is zero. The importance of this assumption is to prevent the relay to join a path that it would harm. This happens in the case when all its repercussion utilities for all paths are negative.

Due to our assumption that the network is interference free, the value of each coalition is not dependent on the other coalitions. It has been proven in [START_REF] Abuzainab | Coalition Formation Game for Cooperative Cognitive Radio Using Gibbs Sampling[END_REF] that a coalition game that satisfies this property and where repercussion utilities are used is an exact potential game with the sum of the original utilities as the potential function. Hence, our game is a potential game where the potential function is the negative of the proportional fair sum of delays of all source-destination pairs. The result in Theorem 1 follows since an exact potential game has the property that (at least) one pure Nash equilibrium exists and that the Nash equilibria are the local maximizers of the potential function.

C. Algorithm

We present in Algorithm 1 a distributed algorithm for our cooperative network formation game in which the relays select to connect to a particular source-destination (s i , d i ) based on proportional fairness. We assume that the relays have full knowledge of the network topology, and that they store the current value of the multihop delay as well as the current path of each source-destination pair. We assume the size of the broadcasted messages is small (i.e. do not require high bit rates) and thus omnidirectional transmission is used during this phase and that a round-robin algorithm is chosen to select each relay r periodically.

Note that we introduce some randomness as we allow each relay to take some non-optimal decision. Indeed, each relay joins the path that yields the maximum repercussion utility with probability 1ε, where ε is commonly known as the mutation probability [START_REF] Kandori | Learning, mutation, and long run equilibria in games[END_REF]. Otherwise, the relay will randomly join the path of any other source-destination pair.

D. Convergence

Due to the mutation probability, the evolution of paths of all source-destination pairs forms a Markov chain which is irreducible and aperiodic. Hence, it has a unique stationary distribution. It is shown in [START_REF] Kandori | Learning, mutation, and long run equilibria in games[END_REF] that as ε tends to zero, the process converges to a unique limiting distribution. Also, since our game is a finite exact potential game, it admits one or several pure Nash equilibria that are the local maximizers of the Nash product. Hence, as ε tends to zero, the algorithm converges to a deterministic Nash equilibrium. In order to reach the global maximum of the Nash product, it is useful to incorporate Gibbs Sampling techniques (such as the algorithm in [START_REF] Abuzainab | Coalition Formation Game for Cooperative Cognitive Radio Using Gibbs Sampling[END_REF]). The drawback of Gibbs Sampling is that the convergence time might be unacceptably large for some scenarios.

IV. NUMERICAL SIMULATIONS

In order to simulate our multihop relaying algorithm, we consider M = 3 source-destination pairs and N = 10 deployed relays. The coordinates of all nodes are generated randomly according to a uniform distribution on a 1000×1000 meters rectangular grid. We set all direct channels between source-destination pairs to be NLOS. All other links are chosen LOS or NLOS randomly. The probability of LOS, p L is taken according to the statistical blockage model of [START_REF] Bai | Coverage and Rate Analysis for Millimeter Wave Cellular Networks[END_REF]: p L = e -βd , where d is the distance between the nodes and β is a parameter known as the average LOS range of the network and is related to the density and average blockage sizes and set to 1 β = 141.4 meters. As for the pathloss models for both cases of LOS and NLOS, we use the the values obtained from the measurements in [START_REF] Rappaport | Broadband millimeter-wave propagation measurements and models using adaptive-beam antennas for outdoor urban cellular communications[END_REF], i.e. α N = 3.88 and α L = 2.20 for the pathloss exponents, A N = A L = 1 for the pathloss coefficients and M T = M R = 4 for all antennas gains. The transmission power for all nodes is set to be P T = 1 Watts and the AWGN variance is set to N 0 = -40.87 dBm. The bandwidth is set to W = 1 GHz and the size of all files is B 1 = B 2 = B 3 = 1 Gb. For the algorithm, we choose the mutation probability ε = 10 -4 .

In order to investigate the potential benefits of proportional fairness, we compare its results to a modified version of the algorithm in which each relay joins the path of the sourcedestination pair that has the minimum delay. This modified version can be interpreted as a greedy approach whose convergence points are the Nash equilibria of the system. Figures 1 and2 show the paths formed between each source destination pair by using our multihop relaying algorithm and the modified minimum delay algorithm respectively. The red, green, and blue circles represent source nodes 1,2, and 3 respectively. The red, green, and blue squares represent destination nodes 1,2, and 3 respectively, and the black circles represent the relay nodes. In this run, most of the formed edges are LOS. Table I shows the delay values computed for the cases of direct transmission, the proportional fair multihop relaying algorithm, and the modified minimum delay algorithm respectively. First, it is easy to see that the delay values have significantly decreased when multihop relaying (for both cases of proportional fairness and minimum delay) is employed, which confirms the benefits of the multihop relaying algorithm in improving the transmission of the different source destination pairs. Also by comparing the delay values obtained from the two algorithms, we find that when using the minimum delay algorithm, the delay of (s 2 , d 2 ) has slightly dropped from 0.2505 sec (for the proportional fairness algorithm case) to 0.1058 sec while the delay of (s 1 , d 1 ) has increased from 0.1428 sec (for the proportional fairness case) to 1.2451 sec. Further, we compute the variances of the delays of all paths for both the proportional fairness algorithm and the minimum delay algorithm. We find out the the variance is 0.282 for the minimum delay algorithm while it is 0.0028 for the proportional fairness algorithm. This shows that proportional fairness can provide a better distribution of the relays among the source destination pairs. In each simulation, we randomly generate the coordinates of the nodes and choose which links are LOS. After each run, we record the sum of delays of all paths for both algorithms. Then, we compute the average sum of delays for each algorithm out of the 1000 runs. We find that the average is 6.8 sec for our algorithm while it is 9.2 sec for the minimum delay algorithm and 86 sec when no multihoping is used (i.e. through the direct path). These values demonstrate the power of proportional fairness compared to the minimum delay approach: while the minimum delay approach can slightly benefit to some users, it does so at the cost of a decreased global performance. Both techniques exhibit excellent global performance compared to the direct transmission. Further, although it converges to a local optima of the Nash product, we observe very good performance metrics. Hence, while Gibbs sampling techniques could ensure convergence to global optima, it is anticipated that the cost of the convergence time will not be compensated by significant performance improvement.

V. CONCLUSION

We have considered a cooperative network formation algorithm in order to construct a multihop path through relays to improve the transmission of source-destination pairs employing millimeter wave RF. Also, we have considered proportional fairness in our algorithm. Due to the assumption of negligible interference in millimeter wave networks, we could show that our network formation game can be turned into a potential game, whose Nash equilibria maximize the proportional fair sum of the transmission rates. We further proposed a distributed algorithm and assessed its performance by numerical simulations. The results show the considerable performance improvement brought by multihop relaying especially in the case where the sources and destinations are NLOS. They also confirm the benefits of proportional fairness in achieving balanced allocations among the source-destination pairs while maintaining a good overall performance of the system.

Algorithm 1 : 2 foreach 4 r 7 with probability 1

 12471 Multihop Relaying Algorithm 1 repeat Relay r in N do 3 foreach source-destination (s k , d k ) in M do computes the repercussion utility of connecting to (s k , d k ): r k (r) = log(D k )log(D k (r)) 5 Let k = argmax log(D k )log(D k (r)) 6 end ε:

8 { 9 Update path P k 10 Relay: 12 { 13 Update path P j 14 Relay 17 r

 891012131417 Relay r connects to (s k , d k ) r broadcasts the updated path P k and the new multihop delay D k (r) to all other relays} 11 otherwiser connects randomly to (s j , d j ) (j = k) r broadcasts the updated path P j and the new multihop delay D j (r) to all other relays } 15 if r was previously connected to different (s i , d i ) then 16Update path P i by removing r broadcasts the updated path P i and the new multihop delay for (s i , d i )

Fig. 1 :Fig. 2 :

 12 Fig. 1: Paths formed by Algorithm 1

TABLE I :

 I Delay Values (in seconds) Further, we run both both algorithms for a thousand times.

This research was supported in part by the French National Research Agency project NETLEARN (ANR13INFR004)


Nicolas Duchateau

Maxime Sermesant

Prediction of infarct localization from myocardial deformation

We propose a novel framework to predict the location of a myocardial infarct from local wall deformation data. Non-linear dimensionality reduction is used to estimate the Euclidean space of coordinates encoding deformation patterns. The infarct location of a new subject is inferred by two consecutive interpolations, formulated as multiscale kernel regressions. They consist in (i) finding the low-dimensional coordinates associated to the measured deformation pattern, and (ii) estimating the possible infarct location associated to these coordinates. These concepts were tested on a database of 500 synthetic cases generated from a realistic electromechanical model of the two ventricles. The database consisted of infarcts of random extent, shape, and location overlapping the whole left-anterior-descending coronary territory. We demonstrate that our method is accurate and significantly overcomes the limitations of the clinically-used thresholding of the deformation patterns (average area under the ROC curve of 0.992±0.011 vs. 0.812±0.124, p<0.001).

Introduction

In clinical routine, imaging of the heart often aims at evaluating the local cardiac tissue viability. A decrease in this viability directly a↵ects the electrical propagation and the muscle contraction, and therefore hampers the resulting cardiac function. However, the transfer function linking the local deformation to the tissue viability is complex, due to advanced physiological interactions between the muscle and the blood, the fibers arrangement, or the influence of the neighboring segments or the opposite wall. Late-enhancement imaging is generally accepted as ground truth to localize the regions with scarred tissue [START_REF] Ortiz-Perez | Correspondence between the 17segment model and coronary arterial anatomy using contrast-enhanced cardiac magnetic resonance imaging[END_REF][START_REF] Sjøli | Diagnostic capability and reproducibility of strain by Doppler and by speckle tracking in patients with acute myocardial infarction[END_REF]. However, this modality is costly and requires the injection of a contrast agent. Moreover its post-processing is still challenging due to the limited contrast and number of acquired slices. On the other hand, 3D echocardiography is non-ionizing and cheaper, and allows the quantification of local wall deformation (myocardial strain). Cardiologists use it in a daily practice to assess the local tissue viability. Nonetheless, the localization of infarct by thresholding deformation patterns [START_REF] Sjøli | Diagnostic capability and reproducibility of strain by Doppler and by speckle tracking in patients with acute myocardial infarction[END_REF] is inaccurate, as the optimal threshold depends on the infarct position and grade, and the relationship between these two parameters is not straightforward.

In this paper, we address this issue by proposing a framework to predict the infarct location from local wall deformation data. While the literature abounds of This a pre-print version. The final document is available at http://www.springerlink.com works that threshold the deformation data, our work is novel in designing a more advanced and accurate prediction strategy. Several works investigated variables inference for di↵erent cardiac applications, such as the reconstruction of fibers architecture from cardiac shape [START_REF] Lekadir | Statistical personalization of ventricular fiber orientation using shape predictors[END_REF] or to relate cardiac shape remodeling and myocardial infarction [START_REF] Zhang | Atlas-based quantification of cardiac remodeling due to myocardial infarction[END_REF]. However, these works used linear or logistic regression techniques, which are not necessarily suited to compute statistics on cardiac motion and deformation patterns. For our concrete application, we preferred a formulation based on spectral embedding [START_REF] Yan | Graph embedding and extensions: a general framework for dimensionality reduction[END_REF] and kernel regression [START_REF] Duchateau | Adaptation of multiscale function extension to inexact matching. Application to the mapping of individuals to a learnt manifold[END_REF][START_REF] Gerber | Manifold modeling for brain population analysis[END_REF] against a Bayesian formulation, in order to minimize the amount of a-priori knowledge in our method.

A notable asset of the proposed work consists in the large database of synthetic cases created to test the methods against a large variety of infarct configurations. Indeed, in such a framework, the transfer function between the model parameters of healthy/damaged tissue and the local deformation is fully controlled, a requisite to test the methods accuracy on ground truth data. Our experiments show that the prediction of infarct location with our method is significantly more accurate than the clinically-used techniques.

Methods

In the following, we denote {(d k , i k )} k2 [1,K] the pairs of local deformation and infarct position data for the set of K training samples. This information is available at each point of the volumetric mesh for each subject, and each of these two parameters is treated as a column vector. In the present work, the infarct position is defined as a binary value 0/1 at each vertex. The relative change in the cell volume with respect to the beginning of the cycle is used as a surrogate for local deformation. A single scalar value is therefore associated to the deformation at each point of the mesh. This was preferred over more advanced measures of local deformation such as the strain tensor, which would require more advanced statistics. However, using scalar deformation data does not limit the concepts demonstrated in the present work: the observed deformation patterns still reflect local infarct (Figs. 1 and4), and focus is kept on the core of the inverse problem.

Non-linear dimensionality reduction (training set)

The (high-dimensional) local deformation data {d k 2 D} k2 [1,K] are first mapped to a Euclidean space of (low-dimensional) coordinates {c k 2 C} k2 [1,K] by means of standard non-linear dimensionality reduction (Isomap [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]). This assumes that there exists a lower-dimensional manifold that can explain the main variations in such data.

In brief, a nearest-neighbors graph is built for the samples {d k } k2[1,K] , using the Euclidean distance as metric. Then, the geodesic distance between each pair of samples is approximated as the shortest path connecting them along the graph, and put into an a nity matrix. The set of coordinates {c k } k2 [1,K] is finally obtained by the diagonalization of a centered version of this geodesic distance matrix.

No restriction is made on the manifold learning technique used. We preferred a non-linear one, as linear operations directly on the deformation data may generate unphysiological patterns [START_REF] Duchateau | Constrained manifold learning for the characterization of pathological deviations from normality[END_REF][START_REF] Gerber | Manifold modeling for brain population analysis[END_REF] (e.g. the linear average of two deformation patterns corresponding to two disjoint infarcts would provide a mixed widespread pattern reflecting a larger but less accentuated infarct, instead of the pattern of an infarct of similar grade at an intermediate location). The Isomap algorithm was preferred over other spectral embedding techniques [START_REF] Yan | Graph embedding and extensions: a general framework for dimensionality reduction[END_REF] as we assumed that the samples distribution follows a uniform random distribution, in contrast with kernel-based methods more relevant for more clustered distributions.

From deformation patterns to infarct prediction (testing set)

Given a new case for which only the local deformation d is known, our method provides a prediction of the infarct position B( î). Its estimation is obtained by two consecutive interpolations formulated as kernel regression, and a multiscale strategy to prevent from artifacts due to non-uniformities in the density of the samples. An overview of the processing pipeline is given in Fig. 1.

From deformation patterns to low-dimensional coordinates Given a new deformation pattern d 2 D, a single-scale formulation would compute its corresponding coordinates ĉ 2 C as:

ĉ = K X k=1 k D (d, d k ) • a k . ( 1
)
This a pre-print version. The final document is available at http://www.springerlink.com

Here a k is the k-th column of the matrix (K D + 1 D I) 1 C, where I is the identity matrix, C = (c 1 , . . . , c K )

T , D is a scalar weight balancing the adherence to the data and the smoothness of the interpolation, and

K D = (k D (d i , d j )) (i,j
) is a kernel-based a nity matrix between the input samples. The kernel function is defined as k D = exp( kd i d j k 2 / 2 D ), D being its bandwidth. Detailed explanations of this single-scale process can be found in [START_REF] Duchateau | Constrained manifold learning for the characterization of pathological deviations from normality[END_REF]. The expression in Eq.1 corresponds to the analytical solution of the following inexact matching problem:

argmin f 2F 1 2 kf k 2 F + D 2 K X k=1 kf (d k ) c k k 2 ! , (2)
where ĉ = f (d), and k.k F stands for the norm on the reproducible kernel Hilbert space F of functions D ! C.

In the current application, this process is iterated from large to small scales by dividing the bandwidth D by a factor 2 at each iteration s, and looking for the remainder function f F (s 1) , where F (s) stands for the s-th scale approximation of the original function f [START_REF] Duchateau | Adaptation of multiscale function extension to inexact matching. Application to the mapping of individuals to a learnt manifold[END_REF]. In practice, scales across iterations range from the overall spread of the samples until getting lower than the average density of the samples [START_REF] Duchateau | Adaptation of multiscale function extension to inexact matching. Application to the mapping of individuals to a learnt manifold[END_REF].

From low-dimensional coordinates to infarct prediction An infarct map î 2 I is estimated from the coordinates ĉ using a multiscale regression process similar to the one described in the previous subsection, where î, ĉ and c k now respectively stand for ĉ, d and d k in Eq.1. The interpolating function is now denoted g. The infarct map therefore corresponds to î = g(ĉ) = g f (d). Data adherence and smoothness are balanced by a weight C , and the kernel function is k C , of bandwidth C . In a single-scale formulation, this corresponds to:

î = K X k=1 k C (ĉ, c k ) • b k , (3)
where b k is the k-th column of the matrix (K C + 1 C I) 1 I, I = (i 1 , . . . , i K ) T , and

K C = (k C (c i , c j )) (i,j
) . This process is also made multiscale, and follows similar rules to the ones used for the estimation of f . Due to this regression process, the values of î at each vertex lie within the continuous [0, 1] interval. Thus, the prediction of the infarct position B( î) is finally obtained by applying a relevant threshold to î (in our case, previously determined by a ROC analysis, as described in Sec.3.2).

Experiments and results

Dataset

Infarct generation A database of 500 synthetic cases was generated to evaluate the methods. Starting from a volumetric tetrahedral mesh of the two ven-This a pre-print version. The final document is available at http://www.springerlink.com tricles (46876 cells and 9673 points for the left ventricle [LV], corresponding to a myocardial volume/mass of 173 mL / 182 g), a fully-connected region corresponding to an infarct of random extent, shape, and location was constructed for each case. The algorithm determined the diseased region iteratively, as illustrated in Fig. 2. It first randomly selected a starting point within the leftanterior-descending coronary territory. This specific territory was retained due to its higher prevalence and agreement in its delineation [START_REF] Ortiz-Perez | Correspondence between the 17segment model and coronary arterial anatomy using contrast-enhanced cardiac magnetic resonance imaging[END_REF]. Note that this only concerns the location of the starting point for the infarct generation algorithm, and that infarcts can spread out of this territory (Fig. 4). Further testing should be extended to other coronary territories and di↵erent geometries to better evaluate the performance of the method. A spherical neighborhood of random radius between 2 and 12 mm was marked as diseased, and a new starting point was randomly selected within this new region. The process was iterated a random number of times (values from 1 to 16). Infarct extents were of 5.2±2.6 mL (3.0±1.5 % of the LV myocardium). A total of 400 cases were used as training set and the other 100 served as testing data.

Electromechanical simulation A realistic electromechanical model [START_REF] Marchesseau | Personalization of a cardiac electromechanical model using reduced order unscented Kalman filtering from regional volumes[END_REF] was then used to simulate the cardiac function along a full cycle of duration 1 s.

This model was previously evaluated on invasive clinical data and has a realistic

This a pre-print version. The final document is available at http://www.springerlink.com behavior. Simulations here use a real anatomical mesh with fibers architecture from an atlas. The contractility and sti↵ness parameters were altered in the random zone defined by our algorithm to model a local infarct [START_REF] Alessandrini | A pipeline for the generation of realistic 3D synthetic echocardiographic sequences: methodology and open-access database[END_REF]. They partially reflect changes in active force and tissue elasticity. These parameters were retained among many others to limit the training data to manageable amounts, and corresponded to the ones with major influence on the deformation [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF]. Border zone was not set to only evaluate the algorithm on binary prediction. Neighbor locations were therefore only passively influenced by the infarct.

Then, local deformation was computed for each tetrahedral cell of the mesh. Figure 3 depicts the deformation values of two infarct configurations at opposite AHA segments, against a healthy case. In this illustrative example, the whole AHA segment was infarcted. Curves represent volume change values averaged over each of the 17 AHA segments of the LV. Notably, deformation patterns are less a↵ected by a mid-anterolateral infarct than by a mid-inferoseptal one, due to the influence of the right ventricle and the surrounding septal regions. Testing our methods on a wide territory of infarct configurations will therefore allow evaluating the sensitivity of the algorithm to marked or moderated alterations of local deformation.

For the sake of simplicity in the infarct prediction process, we limited the input to our algorithm to the deformation data at end-systole (the {d k 2 D} k2 [1,K] in Sec.2), spatially smoothed by a Gaussian filter of bandwidth 1 cm to prevent from inconsistencies due to point-cell correspondences and the non-homogeneity of cell sizes and orientation across the whole volumetric mesh. Note that this might lower the accuracy of the infarct prediction near the border of the ground truth location. The use of the deformation data along the whole cycle may add robustness to the results in more complex configurations, where local post-systolic abnormalities may be more marked or in the presence of asynchronous hearts.

Parameters setting

The number of nearest neighbors in the Isomap algorithm was set to 5. The number of dimensions retained for the estimated coordinates was set to 30. This value corresponded to the limit from which eigenvalues weigh less than 5% of the first eigenvalue. The scalar weights used in the multiscale regression were determined by a leave-one-out procedure, as the value that minimizes the generalization ability (the reconstruction error for samples lying within the range of noise of the training set). Such values were of D = C = 1.

Results

Representative examples of the outputs of our method are shown in Fig. 4, to compare with the thresholding of the deformation data from the same cases, in Fig. 5 (animated version available as Supplementary Material 1 ). The latter notably failed on the transmurality of the infarct location. Qualitatively, our method correctly predicted the infarct location, even for infarcts of small size, location internal to the myocardium, and reduced e↵ect on the deformation curves (e.g. infarcts closer to the septum).

ROC analysis On each case, a ROC analysis was used to determine the optimal threshold leading to the infarct prediction B( î) from the infarct map î 2 I (Fig. 6a). A similar process was applied for the direct thresholding of the deformation data d 2 D (Fig. 6b). These thresholds were defined as the average of the optimal thresholds for each individual case. The performance of directly thresholding the deformation was rather poor in terms of sensitivity and specificity (average area under the curve: 0.812±0.124). Some cases even led to a ROC curve worse than randomly selecting points of the mesh (the diagonal line). Our method significantly outperformed this technique in all cases (average area under the curve: 0.992±0.011, p<0.001).

Conclusion

We presented a method to predict the location of a myocardial infarct from local deformation patterns. This approach is novel and contrasts with the simpler and This a pre-print version. The final document is available at http://www.springerlink.com clinically-used thresholding of the deformation patterns. Notably, our method significantly outperformed this technique. A notable asset of the proposed work also resides in the large database used to test the methods, made of synthetic cases with infarcts of random extent, shape, and location. We are currently collecting a database of 3D echocardiographic sequences and late-enhancement images to extend the evaluation of our method to real data. This will allow evaluating the scalability of the method towards the use of di↵erent heart geometries, more complex deformation patterns, and possibly less contrasted local changes in the deformation data to the acquisition and post-processing of ultrasound sequences.

Fig. 1 .

 1 Fig. 1. Proposed method to predict infarct location from local myocardial deformation.

Fig. 2 .

 2 Fig. 2. Iterative generation of an infarcted region (red ) with random extent, shape, and location, initiated within the mid-anterolateral segment (black arrow).

Fig. 3 .

 3 Fig. 3. Average myocardial deformation for each AHA segment in a healthy case (a), and two infarcts in opposite AHA segments (b and c).

1Fig. 4 .

 4 Fig. 4. Examples of myocardial deformation pattern, ground truth infarct location, and estimated infarct location.

Fig. 5 .Fig. 6 .

 56 Fig. 5. Ground truth infarct location (red ) against the thresholding of the deformation patterns (blue) for the cases shown in Fig.4. Animated version available as Supplementary Material 1 .

This a pre-print version. The final document is available at http://www.springerlink.com

Acknowledgements. The authors acknowledge the European Union 7th Framework

Programme (VP2HF: FP7-2013-611823) and the European Research Council (Me-dYMA: ERC-AdG-2011-291080). They also thank their colleagues R. Mollero and S. Gi↵ard-Roisin for their support on practical aspects of the SOFA simulations, and their This a pre-print version. The final document is available at http://www.springerlink.com

collaborators M. De Craene (Philips Suresnes, France) and E. Saloux (CHU Caen, France) for discussions on these concepts.