Sara Berri

Vineeth Varma

Samson Lasaulce

Mohammed Said Radjef

Correlated equilibria in wireless power control games

come

Introduction

In this work, the notion of correlated equilibrium, which is a generalization of the Nash equilibrium, is applied in the context of power control in wireless network to determine efficient cooperative strategies. Power control in wireless networks has been studied using game theory in literature [START_REF] Goodman | A Power Control for Wireless Data[END_REF][START_REF] Meshkati | A Game Theoretic Approach to Energy-Efficient Power Control in Multicarrier CDMA Systems[END_REF] by characterizing the Nash equilibrium. However, as the Nash equilibrium is often inefficient, introducing the concept of correlated equilibrium will improve the players utility.

Several works apply the concept of correlated equilibrium to the wireless communications paradigm, we will now present some of the relevant papers. The papers [START_REF] Rozario | Optimization of Non-Cooperative P2P Network from the Game Theory Point of View[END_REF][START_REF] Wang | Peer-to-peer file sharing game using correlated equilibrium[END_REF] look at peer-to-peer (P2P) networks, where the behavior of greedy users can degrade the network performance. Here, the authors introduce correlated equilibrium to improve the player utilities. The paper [START_REF] Wu | Cooperative Strategies for Energy aware Ad Hoc Networks: A Correlated Equilibrium Game Theoretic Approach[END_REF] studied the energy efficiency in ad hoc networks, and proposed a cooperative behavior control scheme to determine cooperative strategies, and help noncooperative players to coordinate their strategies using the correlated equilibrium. An efficient broadcasting strategy in wireless ad hoc networks is proposed in [START_REF] Hakami | Efficient broadcastinginslowfadingwirelessad-hocnetworks[END_REF], modeling the interaction among nodes as a game, the action set comprises two actions, to forward or drop the received message from the source. To achieve the correlated equilibrium, linear programming, and a distributed learning algorithm based on the regret matching procedure [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF] are used.

The main objective of this paper is to propose another equilibrium concept (i.e., correlated equilibrium) in the context of wireless power control games, which allows players to obtain a larger equilibrium set and more efficient points in the presence of a correlation device. Assuming that we can add a correlation device to the game, is it possible to create a mechanism such that the equilibrium payoff set of the obtained game includes payoffs that are not in the initial set (game without correlation mechanism)? We provide answers to this question in this paper.

The key contributions and novelty of our paper are as follows:

1. Introduce the concepts of correlated and communication equilibrium to power control in wireless networks. 2. Provide an algorithm to achieve a correlated equilibrium via regret matching. 3. Provide an algorithm to obtain the pareto-optimal correlated equilibrium via Linear programming. 4. An extensive numerical study comparing the efficiency of the proposed correlated equilibrium with the standard Nash equilibrium.

System model

We consider a system comprised K ≥ 2 pairs of interfering transmitters and receivers, where the transmitters want to communicate with their respective receivers. The channel gain of the link between Transmitter i ∈ {1, ..., K} and Receiver j ∈ {1, ..., K} is denoted by g i j ∈ G (here g i j are the real valued channel gains), where G = {g 1 , . . . , g N } represents the alphabet of the possible channel gains. The transmitter i transmits at discrete power level p i ∈ P i = {p 1 i , . . . , p M i }, with p 1 i = P min i and p M i = P max i . Note that G , P i ⊂ R ≥0 ∀i. We denote by ϕ a communication efficiency function which measures the packet success rate as a function of the signal to interference and noise ratio (SINR). It is an increasing function and lie in [0, 1], and is identical for all the users. Let A i denote a finite discrete set of actions that can be taken by player i. In a power control game, this action corresponds to the wireless signal power used by the transmitter i. The SINR at receiver i ∈ K writes as:

SINR i = a i g ii σ 2 + ∑ j =i a j g ji . (1)
where a i ∈ A is the power of the transmitter i and σ 2 is the noise power. Using these notations, the power control game denoted by G , is defined in its normal form as follows.

Definition 2.1 A power control game is a triplet:

G = {K , {A i } i∈K , {u i } i∈K }; (2)
where:

1. K = {1, . . . , K} is the set of players. 2. A i = {a 1 i , . . . , a M i }, is the corresponding power level set of the player i ∈ {1, . . . , K}; with a m sorted such that a m < a m+1 m, a 1 i = P min i and a M i = P max i , are respectively the minimum and maximum transmitting power of player i, M ≥ 1.

3. u 1 , . . . , u K are the function utilities of the K players for a combination of choices a = (a 1 , . . . , a K ) = (a i , a -i ),

where a -i = (a 1 , . . . , a i-1 , a i+1 , . . . , a K ) denotes the power levels of all other players except player i. For each player i, the utility function u i depends on the success of its transmission, which is a function of all players' actions through a, and on the energy spent in transmission a i . Mathematically, the players' utilities are defined by the following formula (3):

u i (a 1 , . . . , a K ) = ϕ(SINR i ) -αa i , (3)
The parameter α > 0, is introduced to weigh the energy cost.

Note that this kind of payoff, different from the traditional energy efficiency (the ratio of the data rate to the power) has been studied in literature [START_REF] Alpcan | CDMA Uplink Power Control as a Noncooperative Game[END_REF], and is relevant when the payoff corresponds to the profit (in terms of money) for each step. In the following section, we introduce the problem of correlated equilibrium from the power control game, where an observer is added to help the players to correlate their actions.

Problem formulation

In the proposed game (2), the users are modeled as rational players, which means they are expected to choose actions from the possible choices to maximize their utilities. An important concept to characterize the outcome of the game is the Nash equilibrium, which states that every player will select an action which maximizes its utility given the actions of every other player. It corresponds to an action profile from which no player has interest in changing unilaterally its action. However, it is well known that the Nash equilibrium does not always lead to the best performance for players [START_REF] Goodman | A Power Control for Wireless Data[END_REF][START_REF] Meshkati | A Game Theoretic Approach to Energy-Efficient Power Control in Multicarrier CDMA Systems[END_REF].

Therefore, other concepts to reach such a more efficient equilibrium needs to be investigated. In this paper, we will study the concepts of correlated equilibrium and communication equilibrium, introduced by Aumann [START_REF] Aumann | Subjectivity and Correlation In Randomized Strategies[END_REF], and Forges [START_REF] Forges | An Approach to Communication Equilibria[END_REF], respectively. These concepts allow players to use an external mediator which provides each of them information about the action to be played in the game. Such a coordination scheme between players may sustain some equilibrium payoffs that are not achievable by an equilibrium without it (Nash equilibrium). The same conclusions hold for the sub-game perfect correlated equilibrium in repeated games, [START_REF] Lehrer | Correlated Equilibria in Two-Player Repeated Games with Nonobservable Actions[END_REF]. Moreover, the correlated equilibrium is simpler to compute than the Nash equilibrium [START_REF] Gilboa | Nash and Correlated Equilibria: Some Complexity Considerations[END_REF]. In the following sub-section, we define the extended game, including an outside observer, and the correlated equilibrium concept.

Power control games using a correlated device

The concept of correlated equilibrium is developed by considering an extended game that includes an outside observer (mediator), which provides each user with a private recommendation regarding which action to perform. The recommendations are chosen according to a probability distribution over the set of action profiles. This probability distribution is called a correlated equilibrium, if the action profile in which all players follow the observer's recommendations is a Nash equilibrium of the extended game. Definition 3.1 Define p as a joint probability distribution on the action profile set A = ∏ i∈K A i . The distribution p is a correlated equilibrium if and only if for every player i ∈ K :

∑ a -i ∈A -i p(a i , a -i )u i (a i , a -i ) ≥ ∑ a -i ∈A -i p(a i , a -i )u i (a i , a -i ), ∀a i , a i ∈ A i . (4)
Inequality ( 4) means that for each user i, choosing power level a i while it received recommendation to choose the power a i , does not provide a higher expected utility. Thus, it is in the best interest for the users to follow the recommended action. The set of correlated equilibria is nonempty, closed and convex in every finite game [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF], and it can include distributions that are not in the convex hull of the Nash equilibria distributions. Indeed, Nash equilibrium is a special case of correlated equilibria, where p(a i , a -i ) corresponds to the product of each individual probability. Thus, the set of Nash equilibrium points are a subset of the set of the correlated equilibrium points.

The correlated equilibrium considers the ability of users to coordinate actions, and computes the optimality by the joint distribution, so it provides a better solution compared with the non-cooperative Nash equilibrium, where each user acts in isolation. Further, the correlation can provide important insights, when we face the problem of equilibrium selection in games admitting multiple equilibria, which could be the case of the proposed game in some setups. Now, if the problem is to add a mediator to the game and enable players to coordinate their actions, but also assume that the mediator receives information from the players before gives them recommendations, we could construct another kind of mechanism, that makes the information coming from players as inputs, and uses them to find the suitable outputs that correspond to the recommendations. In this scenario, the mechanism is also created such that the players have no interest to deviate from the recommendations, which is the communication equilibrium concept. This information exchange could bring performance improvement with respect to the case where players coordinate their actions without reporting any information to the mediator. In the following Section 3.2, we present how could we apply this concept to the power control game (2).

Communication equilibrium in power control games

Here, we assume that the mediator collects information from the players before making them recommendations. In the studied power control game, we assume the players have type sets represented by T i = {(g 1 ii , . . . , g 1 Ki ), . . . , (g N ii , . . . , g N Ki )}. In the case of a wireless channel, the channel gains are randomly chosen following a certain (known) probability distribution. Therefore, the 'type' for each player is chosen randomly according to the given probability distribution over each player's type set. Each player i sends information about his type, i.e., t i ∈ T i to the mediator (the player might lie if it brings some gain). Thus, a communication device consists of a system p of probability distributions p = p(.|t) t∈T =∏ i∈K T i over A = ∏ i∈K A i . The interpretation is that every player i ∈ K reports its type t i = (g ii , g ji ) ∀ j =i to the mediator, which privately recommend a i according to p = p(.|t). The system p defines a communication equilibrium if none of the players can gain by unilaterally lying on its type or by deviating from the recommended action. Definition 3.2 The system p defines a communication equilibrium if:

∑ t -i ∈T -i q(t -i |t i ) ∑ a∈A p(a|t)u i (t, a) ≥ ∑ t -i ∈T -i q(t -i |t i ) ∑ a∈A p(a|t i ,t -i )u i (t, a i , a -i )∀i ∈ K , ∀t i ,t i ∈ T i , ∀a i ∈ A i . (5)
where q(t -i |t i ) = q(t) ∑ t i ∈T i q(t i ,t -i ) , is the subjective probability assigned to the event t -i , that is the actual profile of the other players' types, if t i is the type of i. The definition 3.2 implies that player i does not get a higher expected utility if it lies about its true type t i and reports t i , or if it plays another action a i instead the recommended action a i . In the following section, we propose some techniques to achieve correlated and communication equilibria.

Implementation of communication and correlated equilibria

The sets of correlated equilibria and communication equilibria are the subsets defined by the intersection of the half-spaces given by the inequalities ( 4) and ( 5), respectively. In this section, we investigate methods to obtain correlated and communication equilibria.

Linear programming method

In this paper, among the multiple correlated equilibria, we consider the one that provides the highest social welfare. Thus, the problem reduces to computing a correlated equilibrium that maximizes the sum of the players' expected utilities. In order to characterize that optimal equilibrium, we propose a linear programming method in which the optimization problem of the power control game can be formulated as follows:

max c x A i x > 0 ∀i ∈ K ∑ M K j=1 x j = 1; 0 ≤ x j ≤ 1 ∀ j = {1, . . . , M K } (6)
Where x = (p(a 1 1 , . . . , a 1 K ), . . . , p(a M 1 , . . . , a M K )); c = (∑ i∈K u i (a 1 1 , . . . , a 1 K ), . . . , ∑ i∈K u i (a M 1 , . . . , a M K )).

A i =                                  (u i (a 1 i , a 1 -i ) -u i (a 2 i , a 1 -i )) . . . (u i (a 1 i , a M -i ) -u i (a 2 i , a M -i )) 0 . . . 0
i , a 1 -i ) -u i (a M i , a 1 -i )) . . . (u i (a 1 i , a M -i ) -u i (a M i , a M -i )) 0 . . . 0 0 . . . 0 0 . . . 0 (u i (a 2 i , a 1 -i ) -u i (a 1 i , a 1 -i )) . . . (u i (a 2 i , a M -i ) -u i (a 1 i , a M -i )) 0 . . . 0
i , a 1 -1 ) -u i (a M i , a 1 -i )) . . . (u i (a 2 i , a M -1 ) -u i (a M i , a M -i )) 0 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 0 . . . 0 (u i (a M i , a 1 -i ) -u i (a M-1 i , a 1 -i )) . . . (u i (a M i , a M -i ) -u i (a M-1 i , a M -i )) . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 0 . . . 0 (u i (a M i , a 1 -i ) -u i (a M-1 i , a 1 -i )) . . . (u i (a M i , a M -i ) -u i (a M i , a M-1 -i ))                                  . ( 7
)
p(a i , a -i ) is the probability that the action profile (a i , a -i ) is chosen. Thus, a distribution x * is said to be correlated optimal if it is solution of the linear program [START_REF] Bacci | Understanding Game Theory via Wireless Power Control[END_REF].

In the same manner, we can characterize the optimal communication equilibrium. The solution could be obtained by solving the following optimization problem: Algorithm 1 Algorithm leading to the optimal correlated equilibrium in power control game 1: In the beginning, the mediator chooses a power profile (a i , a -i ) ∈ A according to p * , that is obtained solving the described linear program, using an appropriate method. 2: The mediator informs each user i of the power to choose a i . max ∑ a∈A t∈T q(t)p(a|t) ∑ i∈K u i (t, a).

(

) 8
subject to constraint ( 5)

∑ t -i ∈T -i q(t -i |t i ) ∑ a∈A p(a|t)u i (t, a) ≥ ∑ t -i ∈T -i q(t -i |t i ) ∑ a∈A p(a|t i ,t -i )u i (t, a i , a -i )∀i ∈ N , ∀t i ,t i ∈ T i , ∀a i ∈ A i . (9)
and for all t ∈ T .

∑ a∈A p(a|t) = 1. ( 10)

p(a|t) ≥ 0. ( 11
)
In the following we summarize the different steps to reach the optimal communication equilibrium.

Algorithm 2 Algorithm leading to the optimal communication equilibrium in power control game 1: The mediator simulates the sequence of reports that could receive from the users, that correspond to channel gain profiles, and the power profiles that could be received by the users given a channel gain profile. 2: Using a method to solve the linear program constituted by the objective function [START_REF] Goodman | A Power Control for Wireless Data[END_REF] and the constraints ( 9)- [START_REF] Forges | Correlated Equilibria and Communication in Games[END_REF], to find the optimal probability distributions p(.|t) for all type profile t. 3: For each player i, the Nature randomly chooses type, that corresponds to the channel gain profile (g ii , . . . , g Ki ), according to a given probability distribution over the type set T i . 4: Each player i reports its type, t i , to the mediator. 5: The mediator performs lotteries according to the received type profile, and sends private recommendations to the players, that corresponds to the transmitting power.

However, with the linear programming method, the computation complexity grows exponentially with the number of users and actions since an increase in the number of users and actions, results an increase in the number of constraints. There exists a distributed learning approach, i.e., regret matching [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF] to achieve a correlated equilibrium. However, regret matching does not ensure the pareto optimality of the given correlated equilibrium as it converges to an arbitrary correlated equilibrium, whereas the linear program as defined in [START_REF] Bacci | Understanding Game Theory via Wireless Power Control[END_REF] gives a correlated equilibrium maximizes the social welfare. In the following, we present the distributed learning approach.

Regret matching procedure

A game procedure is proposed in [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF], called 'regret-matching'. In which the players measure the regret for not choosing other actions in the past, and change their current action with probabilities that are proportional to these measures. Thus, the game is played with probability distribution over the action set. The details of the regret-matching algorithm [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF] is shown in 3.

Algorithm 3 Regret matching algorithm for power control game 1: For any two distinct power levels a i = a i ∈ A i calculate the average regret of user i at time t for not choosing a i as:

R t i (a i , a i ) = max{ 1 t ∑ τ≤t [u i (a i , a τ -i ) -u i (a τ i )], 0}. (12)
2: Let a i ∈ A i the last power chosen by user i, a t i = a i . Then the probability distribution over the power levels for the next period, is defined as

p t+1 i (a i ) = 1 µ R t i (a i , a i ), ∀a i = a i ; p t+1 i (a i ) = 1 -∑ a i ∈A i a i =a i p t+1 i (a i )
, where µ is an enough large constant. 3: Let For every t, the empirical distribution of the power profile a is:

p t (a) = 1 t |{τ ≤ t : a τ = a}|. ( 13
)
where |{τ ≤ t : a τ = a}| is the number of times the power profile a has been chosen in the periods before t.

It is shown in [START_REF] Hart | A Simple Adaptive Procedure Leading to Correlated Equilibrium[END_REF] that if the players follow the proposed algorithm 3, the empirical distribution (13) converges to an arbitrary correlated equilibrium of the game if it is not unique. The obtained correlated equilibrium by applying this procedure is not always Pareto-optimal (but the procedure can achieve the PO equilibrium). However, the one provided by the linear programming method is Pareto-optima.

Numerical Results and Analysis

In this section, we present numerical results. The simulation setup is as follows. The number of Transmitter and Receiver pairs is equal to 2, K = {1, 2}. Set of possible powers, that is the action set A i , ∀i ∈ K : M = 25, P min i (dB) = -20, P max i (dB) = +20. The presented results correspond to the expected values over different values of (g 11 , g 12 , g 22 , g 21 ) that lie in a discrete set G i = {g 1 i , . . . , g N i }, with g 1 i = g min i , g N i = g max i , ∀ ∈ iK : N = 10, g min i = 0.01, g max i = 3, the channel gain increment equals 3-0.01 10 . The means of the channel gains are given by: ( ḡ11 , ḡ12 , ḡ22 , ḡ21 ) = (1, 1, 1, 1). The communication efficiency function is: ϕ(x) = (1e -x ) L , L being the number of symbols per packet (see e.g., [8][18][15]). In most of the simulations provided we take L = 100. The aforementioned parameters are assumed, otherwise they are explicitly mentioned in figures.

Fig 1 compares the payoffs attained by the regular Nash equilibrium and the correlated equilibrium with the region of all possible payoffs (in a centralized setting). We observe that the correlated equilibrium can improve the utility of both the players when compared to the Nash. We are limited to the two action case, but this is a very relevant case as seen from literature [START_REF] Larrousse | Coded Power Control: Performance Analysis[END_REF]. Fig 2a plots the utility of the correlated equilibrium and the communication equilibrium for higher action sets, but we plot to a maximum of four actions due to complexity issues in evaluating the communication equilibrium.

Conclusion

We have introduced the concepts of correlated and communication equilibria to the paradigm of wireless power control games and provided algorithms based on linear programming and regret matching to achieve these equilibria. Our numerical results indicate that the resulting equilibria can be pareto-optimal and outperform the classical Nash equilibrium. While using linear programming to get the communication equilibrium results in a pareto-efficient solution, it has a computational complexity that is exponential in the action set size and number of players. Regret matching is computationally less demanding but results in solutions that are not necessarily pareto-optimal.

  Fig 2b compares the solution of regret matching with that of correlated equilibrium when the action sets are of size 25. Fig 2b demonstrates that although the regret matching algorithm is computationally fast, the solutions are not pareto-optimal.

  [START_REF] Lasaulce | Introducing Hierarchy in Energy Games[END_REF][START_REF] Lasaulce | Introducing Hierarchy in Energy Games[END_REF][START_REF] Lasaulce | Introducing Hierarchy in Energy Games[END_REF][START_REF] Lasaulce | Introducing Hierarchy in Energy Games[END_REF], g min = 0.01, g max = 30, |G| = 50A1 = A2 = {0.01, 80}, α = 10 -2 , L = 1, σ 2 = 1.(b) ..

Fig. 1 :

 1 Fig.1: Sets of: possible payoffs, correlated equilibrium payoffs and Nash equilibrium payoffs for two settings of parameters.

P 2

 2 g min = 0.01,g max = 3,|Gi| = 3 P min =-20(dB),P max =20(dB), L = 100, α = 10 -2 ,σ 2 = 1, K = 2.(a) Expected sum-utility of correlated and communication equilibria against different number of actions |A i |. Larger action sets are not plotted due to complexity issues. min = -20(dB), P max = 20(dB), |Ai| = 25 α = 10 -2 , L = 50, σ g min = 0.01, g max = 3, |Gi| = 10 (b) Sets of correlated equilibrium payoffs and Nash equilibrium payoffs. The correlated equilibrium provided by regret matching is also indicated.

Fig. 2 :

 2 Fig. 2: Equilibrium payoffs for larger action sets.

Research Unit LaMOS (Modeling and Optimization of Systems), Faculty of Exact Sciences, University of Bejaia, Bejaia, 06000, Algeria radjefms@gmail.com

• 2 L2S (CNRS-CentraleSupelec-Univ. Paris Sud), Gif-sur-Yvette, France {sara.berri, samson.lasaulce}@l2s.centralesupelec.

fr • 3 CRAN-ENSEM, Nancy, France vineeth.satheeskumar-varma@univ-lorraine.fr


Chao Zhang

Vineeth Varma

Samson Lasaulce

Robust power modulation for channel state information exchange

This contribution is on the framework of a wireless network with multiple interfering transmitter-receiver pairs. We study the case where there is no direct communication channel between the multiple transmitters and the system is completely distributed (in terms of information available and decision making). Exchanging local channel state information (CSI) among the transmitters is one solution to the problem of improving the efficiency of the network. We introduce a novel power modulation scheme that facilitates reliable exchange of local CSI through the signal power strength assuming feedback of received signal strength indicator from each receiver to its transmitter. Numerical results demonstrate the value of our approach.

Introduction

This contribution is in the framework of an interference network, i.e., a wireless network with multiple transmitter-receiver pair that cause interference with each other. Typically, these transmitting and receiving devices are distributed in decision making and possibly also in terms of the information available. One of the major challenges in such networks is designing algorithms or procedures for controlling the transmit power levels of the wireless signals in order to improve the data rate (or other relevant performance metrics). The iterative water-filling algorithm (IWFA) [START_REF] Yu | Distributed multiuser power control for digital subscriber lines. Selected Areas in Communications[END_REF][2] [START_REF] Mertikopoulos | Distributed learning policies for power allocation in multiple access channels. Selected Areas in Communications[END_REF] is a state-of-the art technique which relies only on locally availble information, the individual signal-to-interference plus noise ratio (SINR), and has a low computational complexity. On the other hand, one drawback of IWFA and many other similar distributed and learning algorithms (see e.g., [START_REF] Rose | Learning equilibria with partial information in decentralized wireless networks[END_REF][5] [START_REF] Bacci | Game theory for networks: A tutorial on gametheoretic tools for emerging signal processing applications[END_REF]) is that convergence is not always ensured [START_REF] Mertikopoulos | Distributed learning policies for power allocation in multiple access channels. Selected Areas in Communications[END_REF] and the result is typically globally inefficient.

To exploit the available feedback signal efficiently, a novel technique is given by [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF]. In [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], instead of using local observations (namely, the SINR feedback realizations or received signal strength indicator (RSSI) realizations) to allow the transmitters to converge to a Nash point, one can use them to acquire global channel state information (CSI). To obtain global CSI, the key idea in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF] is to exploit the transmit power levels as information symbols and the interference as a communication mean or channel for the transmitters.

However, the power modulation scheme proposed in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF] results in non-negligible decoding errors. A novel technique is proposed in this paper which allows to improve the decoding efficiency, i.e., to estimate the symbol coded in the transmitted power of interfering users more precisely. Furthermore, it is interesting to note that the mapping method in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF] is only a special case of the novel technique in this paper.

The chief contributions and novelty of this work are as follows: ◮ We introduce the novel idea of power modulation by considering the property of the RSSI measurements and the decoding rule. This allows the transmitter to estimate the transmitted power of interfering transmitters more easily. ◮ If local CSI is perfectly known by each user, the conditions to completely reconstruct the transmitted power of the interfering transmitters are given. ◮ The proposed technique accounts for noise in the RSSI measurements (the corresponding modeling is provided in Sec. II), while IWFA-like algorithms typically assume noiseless measurements (with a few exceptions [START_REF] Amod | Robust rate maximization game under bounded channel uncertainty[END_REF][3][9]), . ◮ We also provide numerical simulations that justify the proposed scheme.

System model

We consider a wireless network with K ≥ 2 pairs of interfering transmitters and receivers (each pair will be referred to as a user). The technique is described for the case of two users and for a single band for ease of exposition, but an extension of this to a larger number of users or bands is possible as explained in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF].

We denote the channel gain of the link between Transmitter i ∈ {1, ..., K} and Receiver j ∈ {1, ..., K} by g i j ∈ R ≥0 (this is the fading in signal power domain and not in the amplitude domain). We use a K × K channel matrix G whose entries are given by the channel gains g i j . Each channel gain, and therefore the channel matrix itself, is assumed to obey a classical block-fading variation law, i.e., the channels are assumed to be constant over a certain frame, where a frame comprises T I + T II + T III consecutive time-slots where T m ∈ N, m ∈ {I, II, III}, corresponds to the number of time-slots of Phase m of the proposed procedure described in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF]. The first phase (with duration of T I ) is reserved for estimating the local channels, i.e. each user i will estimate g ji for all j. A technique to estimate this with RSSI feedback is described in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], but the focus of this contribution is to improve the second phase, which corresponds to the time in which each user broadcasts information of his local channel gains to all the other users via his transmit power level 1 .

Transmitter i, i ∈ {1, ..., K}, can control its power from time-slot to time-slot and the corresponding power level at time slot t is denoted by p i (t) ∈ [0, P max ], with P max being the maximum transmit power. p(t) = (p 1 (t), ..., p K (t)) T represents the K-dimensional column vector formed by the transmit power levels, T standing for the transpose operator.

We assume the existence of a feedback mechanism which provides each transmitter, a noisy version of the signal power received at its intended receiver for each time-slot. The signal strength observed by Receiver i at time-slot t is expressed as

ω i (t) = g ii p i (t) + σ 2 + ∑ j =i g ji p j (t) (1)
where σ2 is the noise variance. Receiver i measures the received signal (RS) power ω i (t) (for each time slot) and quantizes it using N bits (the RS power quantizer is denoted by Q RS ). The Receiver i then sends the quantized RS power ω i (t) as feedback to Transmitter i via a noisy feedback channel. Quantization yields ω i (t) ∈ W , for all i ∈ {1, ..., K}, where

W = {w 1 , w 2 , . . . , w M } such that 0 ≤ w 1 < w 2 < • • • < w M and M = 2 N .
The feedback channel is represented by a discrete memoryless channel (DMC) whose conditional probability is represented by Γ . The distorted and noisy version 2 of ω i (t), which is available at Transmitter i, is denoted by ω i (t) ∈ W ; the quantity ω i (t) will be referred to as the received signal strength indicator (RSSI). Thus, the probability that Transmitter i decodes the symbol w ℓ if Receiver i sent the quantized RS power w k equals Γ (w ℓ |w k ).

Each user i has its own individual utility u i which is a function of the global channel matrix and all the transmit powers. Therefore, main objective of each user is to maximize its individual utility u(p; G) (for example, the data rate) by controlling its signal power p i . However due to the distributed nature of the system, user i does not have access to G and therefore can not necessarily make the optimal decision. In [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], the authors propose a novel method of exchanging local CSI resulting in each user acquiring G. We detail our contribution to this method in the following section.

Exchanging local CSI via power modulation

In this contribution, we focus on improvements to Phase II of the scheme proposed in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], specifically with regards to power modulation. The technique in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF] is to exploit the transmit power levels as information symbols and exploit the observed interference (which is observed through the RSSI or SINR feedback) for inter-transmitter communication. The corresponding implicit communication channel is exploited to acquire global CSI knowledge namely, the matrix G and therefore to perform operations such as the maximization of u(p; G).

The process of achieving the desired power control vector is divided into three phases . In Phase I, a sequence of power levels which is known to all the transmitters is transmitted (similar to a training sequence in classical channel estimation but in the power domain), and Transmitter i estimates its own channel gains (i.e., g 1i , g 2i , ..., g Ki ) by exploiting the noisy RSSI feedback; we refer to the corresponding channel gains as local CSI. In Phase II, each transmitter informs the other transmitters about its local CSI by using power modulation. By decoding the modulated power, each transmitter can estimate the channel gains of the other users and thus, at the end of Phase II each transmitter has its own estimate of the global CSI G; the situation where transmitters have a non-homogeneous or different knowledge of global CSI is referred to as a distributed CSI scenario in [START_REF] De | Degrees of freedom of the network mimo channel with distributed csi[END_REF]. In Phase III, each transmitter can then exploit global CSI to maximize (possibly in a sub-optimal manner) the network utility of interest.

The key idea of this paper is to modify the basic power modulation method to enhance the estimation quality. Firstly, we introduce the power level decoding rule and the basic power modulation method proposed in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF].

Basic power modulation scheme

As described before, only 2 users are active at any given time-slot in Phase II for our study. To inform the other transmitters about its knowledge of local CSI, Transmitter i maps the K labels of N II bits produced by the quantizer Q II i to a sequence of power levels (p i (T I + 1), p i (T I + 2), . . . , p i (T I + T II )). As described in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], Phase II comprises T II = 2 time-slots, K = 2 users, and that the users only exploit L = 2 power levels during Phase II say P = {P min , P max }. Further assume 1-bit quantizers, which means that the quantizers Q II ji produce binary labels. For simplicity, we assume the same quantizer Q is used for all the four channel gains g 11 , g 12 , g 21 , and g 22 : if g i j ∈ [0, µ] then the quantizer output is denoted by g min ; if g i j ∈ (µ, +∞) then the quantizer output is denoted by g max . Therefore a simple mapping scheme for Transmitter 1 (whose objective is to inform Transmitter 2 about (g 11 , g 21 )) is to choose p 1 (T I + 1) = P min if Q(g 11 ) = g min and p 1 (T I + 1) = P max otherwise; and p 1 (T I + 2) = P min if Q(g 21 ) = g min and p 1 (T I + 2) = P max otherwise. Therefore, depending on the p.d.f. of g i j , the value of µ, the performance criterion under consideration, a proper mapping can chosen. For example, to minimize the energy consumed at the transmitter, using the minimum transmit power level P min as much as possible is preferable; thus if Pr(Q(g 11 ) = g min ) ≥ Pr(Q(g 11 ) = g max ), the power level P min will be associated with the minimum quantized channel gain that is Q(g 11 ) = g min .

At every time-slot t ∈ {T I + 1, ..., T I + T II }, the power levels of the interfering transmitter are estimated by Transmitter i as follows

p -i (t) ∈ arg min p -i ∈P K-1 ∑ j =i p j g ji -( ω i (t) -p i (t) g ii -σ 2 ) (2)
where p -i = (p 1 , .., p i-1 , p i+1 , .., p K ). As for every j, g ji is known at Transmitter i, the above minimization operation can be performed. When the number of users is higher, each transmitter needs to estimate K -1 power levels with only one observation equation, which typically induces a non-negligible degradation in terms of symbol error rate. In this situation, Phase II can be performed by scheduling the activity of all the users, such that only 2 users are active at any given time-slot in Phase II. Once all pairs of users have exchanged information on their channel states, Phase II is concluded.

An adaptive modulation scheme with local CSI perfectly known

Without loss of generality, we consider the case with K=2. For consistency, we assume 1-bit quantizer of g ji , which means that the quantizer Q II ji produces binary labels. For transmitter i, the objective is to estimate the power emitted by another transmitter i ′ with i ′ = i. With the proposed decoding method in the previous section, it can be seen that the decoding error is mainly brought by the consequence that with different power levels emitted, the same quantized symbols (RSSI feedback) will be achieved under some extreme conditions, such as g ji is very close to 0. To avoid the occurrence that the RSSI feedbacks lie in the same quantization intervals with different transmitted power, we propose a novel power modulation scheme here. For simplicity, we analyze the received signal of receiver 1 to introduce our novel technique. Assume Phase II comprises T II = 2 time-slots, where p 1 (T I + 1) depends on g 11 and p 1 (T I + 2) depends on g 21 . The basic idea here, by denoting j ∈ {1, 2}, is to transmit with power p 1 (T I + j) which is larger than P max /2 if Q(g j1 ) = g max , otherwise the power emitted by transmitter 1 p 1 (T I + j) should be less than P max /2, . The most important point is to choose the power p 1 (T I + j) such that ω 1 (T I + j) = g 11 p 1 (T I + j) + g 21 p 2 (T I + j) + σ 2 lie in different quantization intervals with p 2 (T I + j) ≥ P max /2 and p 2 (T I + j) < P max /2. Intuitively, it can be easily found that if g 11 p 1 (T I + j) + 1 2 × g 21 P max + σ 2 exactly equals to a quantization interval bound, the point mentioned above can be perfectly fulfilled. Denoting the quantization interval bounds of RSSI quantizer by {t 1 ,t 2 , ...,t 2 M +1 } (where ω ∈ [t m ,t m+1 ) is quantized as w m ), the power p 1 (T I + j) can be selected as follows:

p 1 (T I + j) =            min 1 2 P max + 1 g 11 t a -1 2 g 21 P max -1 2 g 11 P max -σ 2 , P max if Q(g j1 ) = g max min 1 g 11 (t b -1 2 g 21 P max -σ 2 ), 1 2 P max else (3)
where the index a fulfills t a-1 ≤ 1 2 g 21 P max + 1 2 g 11 P max + σ 2 < t a and the index b fulfills t b-1 ≤ 1 2 g 21 P max + σ 2 < t b .

Proposition 1. Assume the local CSI is perfectly known by each user and the effect of the largest quantization interval is ignored, if the two conditions below are always satisfied:

(i) 1 2 P max + 1 g 11 t a -1 2 g 21 P max -1 2 g 11 P max -σ 2 ≤ P max (4)

(ii) 1 g 11 (t b -1 2 g 21 P max -σ 2 ) < 1 2 P max ( 5
)
then transmitter 1 can always reconstruct the power emitted by transmitter 2 without error, i.e. p 2 can be decoded perfectly by transmitter 1.

Proof. If g 11 p 1 (T I + j) + g 21

P max
2 + σ 2 can always achieve the boundary of the quantization interval, the p 2 can be perfectly decoded since 2 possible values of p 2 are belongs to (0, P max 2 ) and ( P max 2 , P max ) respectively. The first condition above indicates that there exists p 1 (T I + j) ∈ ( P max 2 , P max ) such that g 11 p 1 (T I + j) + g 21 P max 2 + σ 2 equals to a boundary of the quantization interval, and the second condition above indicates that there exists p 1 (T I + j) ∈ (0, P max 2 ) such that g 11 p 1 (T I + j) + g 21

P max
2 + σ 2 equals to a boundary of the quantization interval.

Considering the 2 conditions above, it can be observed that the conditions are difficult to be fulfilled if the channel gain g 11 is too small. Even it is known that g 11 follows the exponential distribution in our scenario, it is reasonable to set a minimum value g min for the channel gain since it is usually met in practical scenario.

Proposition 2. Assume the local CSI is perfectly known by each user and the effect of the largest quantization interval is ignored, if each realization of the channel gain g 11 is larger than g min = 2∆

P max where the longest length of the quantization interval(except the interval

(t 2 N ,t 2 N +1 )) is denoted as ∆ = max i∈{1,...,2 N -1} (t i+1 -t i ),
then transmitter 1 can always reconstruct the power emitted by transmitter 2 without error. Specially, if the RSSI is uniform quantized with step d, then the minimum value can be expressed as g min = 2d P max .

Proof. Omitted.

An adapted modulation scheme with noisy local CSI

However, the scheme detailed in the previous section is not robust to noisy local CSI. Since the local CSI can not be perfectly acquired (g i j is not perfectly known to the transmitters), it would be difficult to calculate the p 1 (T I + j) such that

g 11 p 1 (T I + j) + 1 2 g 21 P max + σ 2 = t a ′ (6)
where a ′ ∈ {2, ..., 2 N }. Using the same rule in Sec 3.2, with noisy (estimated) g 11 and g 21 , the p 1 (T I + j) can be only obtained by:

g 11 p 1 (T I + j) + 1 2 g 21 P max + σ 2 = t a ′ (7)
Therefore, it can not be guaranteed that the real RSSI when the other transmiiter uses P max /2 exactly equals to the interval bound due to the noise in the local CSI. To compromise for this, we propose a heuristic technique to improve the decoding rate by combining the schemes used in Sections 3.1 (which is robust to the noise in local CSI estimation but has lower decoding success rate even with perfect local CSI) and Section 3.2 (which has a very high decoding success rate with perfect local CSI, but is not robust to noise). We are able to achieve this by introducing an "offset" x, which we will describe in this section. Define g 11 = g 11 + z 11 and g 21 = g 21 + z 21 (the noisy estimate of local CSI available at the transmitter), ( 7) can be rewritten as:

g 11 p 1 (T I + j) + 1 2 g 21 P max + σ 2 + z 11 p 1 (T I + j) + 1 2 z 12 P max = t a ′ (8)
Nevertheless, no matter what is the distribution of the noise, it can be found that with a probability q that Pr |z 11 p 1 (T I + j) + 1 2

z 12 P max | ≥ c ≤ q (9)
where c is a constant related to the probability q and the noisy local CSI. Combine (8)(9), it can be concluded Pr(A) ≥ q where the event A is defined as

t a ′ -c ≤ g 11 p 1 (T I + j) + 1 2 g 21 P max + σ 2 ≤ t a ′ + c (10)
Equivalently, the (10) can be rewritten as:

g 11 p 1 ((T I + j) + g 21 ( 1 2 P max + c g 21 ) + σ 2 ≥ t ′ a ( 11
)
g 11 p 1 ((T I + j) + g 21 ( 1 2 Pmax - c g 21 ) + σ 2 ≤ t ′ a ( 12
)
(11) and ( 12) imply that even though the noisy local CSI reduce the accuracy of the emitted power, the received signal ω 1 (T I + j) = g 11 p 1 (T I + j) + g 21 p 2 (T I + j) + σ 2 lie in different quantization intervals with probability q ′ ≥ q for p 2 (T I + j) ≥ ( 12 Pmax + c g 21 ) and p 2 (T I + j) ≤ ( 1 2 Pmax -c g 21 ). To improve the robustness of the transmission, the offset x (x ≤ 0.5) will be implemented such that transmitting with power p 1 (T I + j) which is larger than ( 1 2 + x)P max if Q(g j1 ) = g max , transmitting with power p 1 (T I + j) which is less than ( 1 2x)P max if Q(g j1 ) = g min . Induced by noisy local CSI, the power p 1 (T I + j) can be redetermined as follows: Transmitting with 0 and P max is a special case of the new scheme corresponding to x = 0.5. The optimal x can be calculated by performing numerical simulations.

p 1 (T I + j) =            min ( 1 2 + x)P max + 1 g 11 t * a -1 2 g 21 P max -( 1 2 + x)g 11 P max -σ 2 , P max if Q( g j1 ) = g max min 1 g 11 (t * b -1 2 g 21 P max -σ 2 ), (
Remark 1. In the Section 3.2 and 3.3, it is assumed that the error of DMC is negligible, i.e. ω i = ω i . The approach to reduce the influence of the error induced by DMC will be studied in future works.

Remark 2. (extension to the multi-band scenario). As explained in the beginning of this paper, each band performs in parallel like the single-band case. Since there are power constraints for each transmitter, the modulated power should satisfy ∑ S s=1 p s i (t) ≤ P max .

Numerical analysis

In this section, as a first step, we start with providing simulations which result from the adapted power modulation when local CSI is perfectly known. To make a coherent comparison with the simple power modulation in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], the quality of decoding will be evaluated by considering the distortion and sum-rate utility. As a second step, we study the scenario with imperfect local CSI and analyze the influence of the offset, which is introduced to assure the robustness of the system. Firstly, we focus on the scenario with perfect CSI. The novel modulation scheme is designed to reduce the error probability of decode. To estimate the decoding quality, the decoding success rate (DSR) is introduced, which represent the probability that the transmitter i correctly reconstruct the transmitted power by all the other transmitters in all time-slots. For the sake of simplicity, we consider the K=2, S=1 scenario. The channel gain of direct channel g ii follows the exponential distribution with expectation 1, where the channel gain of cross channel g ji ( j = i) follows the exponential distribution with expectation 0.1, i.e. f (g ii ) = exp(-g ii ) and f (g ji ) = 10 exp(-10g ji ). Furthermore, each channel is assumed to be independent. We use the uniform quantizer for the RSSI feedback and Maximum entropy quantizer (see [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF]) for the channel gain in this part. In Fig. 1, it can be seen that the adapted power modulation achieves a much better performance in terms of DSR with both low and high quantization bits. Moreover, with the adapted modulation scheme, Fig. 1 shows that the DSR increases when we have more quantization bits, and it tends to 1 when the number of quantization bit is large.

Secondly, we take into account the noisy local CSI scenario. As described in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF], the local CSI can be estimated by RSSI feedback with training matrix. To distinguish the noise level here, we define the high noise level when we use 4 quantization bits of RSSI estimating the local CSI and the low noise level when we use 8 quantization bits of RSSI estimating the local CSI. Fig. 2 illustrate the offset is useful to reduce the estimation distortion, which is defined by

Distortion = E[ G -G i 2 ]
Also, it can be observed that different noise level has different optimum offset value. When noise level is low, the optimum offset value is between 0.2 and 0.3, where 0 and P max becomes a good solution when noise level is high since less choices can avoid the decoding error. Meanwhile, as the objective of the adapted modulation scheme in perfect local CSI scenario, the red curve illustrates our scheme (offset 0) beats the simple modulation scheme (offset 0.5).

Conclusion

In this contribution, we have proposed a novel power modulation scheme that improves the inter-transmitter communication efficiency of the technique introduced in [START_REF] Vineeth | Power modulation: Application to inter-cell interference coordination[END_REF]. Under the assumption of noiseless feedback of the quantized received signal strength, our simulations show that a significant improvement in the distortion or decoding success rate can be achieved both for the case where perfect local CSI and noisy local CSI is available. A relevant extension to this scheme would be to account for the noise in feedback (in the DMC).

Fig. 1

 1 Fig.1K=2, S=1, SNR=30dB,ε = 0, perfect local CSI, 1 quantization bit for channel gain, we observe that the adapted modulation scheme brings higher DSR in low and high resolution case.

  fulfills t a * -1 < 1 2 g 21 P max + ( 1 2 + x) g 11 Pmax + σ 2 ≤ t * a and the index b * fulfills t b * -1 < 1 2 g 21 P max + σ 2 ≤ t * b .

	1 2 -x)P max	else
		(13)
	where the index a

*

Note that "normal" or standard communication between the transmitter and receiver of every user occurs at every phase independent of what is done in Phase I and II.

Note that, for the sake of clarity and ease of exposition,we assume the RS power quantizer and DMC to be independent of the user index, but the proposed approach also holds in the general case.


François Lemaire
email: francois.lemaire@univ-lille1.fr

Alexandre Temperville
email: a.temperville@ed.univ-lille1.fr

On Defining and Computing "Good" Conservation Laws

Keywords: conservation analysis, sparse conservation laws, biological models, sparse null space, greedy algorithm. 1 index j

Conservation laws are a key-tool to study systems of chemical reactions in biology. We address the problem of defining and computing "good" sets of conservation laws. In this article, we chose to focus on sparsest sets of conservation laws. We present a greedy algorithm computing a sparsest set of conservation laws equivalent to a given set of conservation laws. Benchmarks over a subset of the curated models taken from the BioModels database are given.

Introduction

Many biological processes can be modelled by systems of chemical reactions. In order to study such systems, one usually computes its (linear) conservation laws (i.e. linear combinations of number of species) which have the property of being constant along the time. In this article, we only consider (linear) conservation laws. For a given chemical reaction system, a complete set of conservation laws is easily computed by computing a basis of the kernel of the transpose of the stoichiometry matrix of the system [START_REF] Herbert | Conservation analysis in biochemical networks: computational issues for software writers[END_REF].

This paper tries to answer to the difficult question: "what is a good conservation law ?". Consider for example the well known enzymatic degradation given by the reactions E + S ↔ C and C → E + P . It admits for example E + C and S + C + P as conservation laws. One could as well consider their sum (i.e. E + 2C + S + P ), their difference (i.e. E -S -P ), . . . On the example, the laws E + C and S + C + P seem less artificial and closer to the physics of the system, than the two laws E + 2C + S + P and E -S -P . Indeed the law E + C corresponds to the conservation of the enzyme E, and S + C + P corresponds to the conservation of the substract S (which is either in the form S, C or P ).

In an attempt to define "good conservation laws", we think that a good conservation law should have many zero coefficients (i.e. sparse laws) and many positive coefficients. Concerning the sparse property, we think that a practitioner would understand better sparse laws than dense laws, since sparse laws are shorter and thus easier to read. Moreover, a sparse conservation law can also be useful when doing substitutions in differential equations to preserve the sparsity of the differential equations. For example, if one has a sparse system Ẋ = F (X) (where X is a vector of species X 1 , X 2 , . . . ), one can use a conservation law involving X 1 (say X 1 + X 5 -X 8 = c 0 ) to discard the variable X 1 by substituting X 1 by an expression in the other species (c 0 -X 5 + X 8 on the example). Consequently, a sparse conservation law will more likely preserve the sparsity of the differential equations. Concerning positive coefficients, we think that conservation laws with positive coefficients are more likely to represent a conservation of matter.

Those two criteria (sparsity and positiveness) are sometimes impossible to satisfy at the same time. For example, if we have a basis of two conservation laws X 1 + X 2 + X 3 and X 2 + X 3 + X 4 , then the difference X 1 -X 4 is sparser than any of the two laws but involves a negative coefficient. Moreover, in some particular examples, there are no conservation laws with positive coefficients only (like in A + B → ∅ which only admits A -B as a conservation law).

In this paper, we have chosen to compute a sparsest basis of conservation laws, leaving the positivity argument for a further work. As a consequence, our approach differs from computing minimal semi-positive P-invariants (i.e. conservation laws with non-negative coefficients with minimal support [START_REF] Soliman | Invariants and Other Structural Properties of Biochemical Models as a Constraint Satisfaction Problem[END_REF]).

Our approach corresponds to the well known Nullspace Sparsest Problem (NSP) which is proven to be NP-hard in [START_REF] Thomas | The null space problem i. complexity[END_REF]. NSP consists in finding a matrix with the fewest nonzeros, whose columns span the null space of a given matrix. Approximate algorithms to solve NSP are given in [START_REF] Thomas | The null space problem i. complexity[END_REF][START_REF] Thomas | The null space problem ii. algorithms[END_REF].

We chose to develop our method by testing it on the Biomodels database [START_REF]Biomodels database[END_REF]. Our hope was that biological models might have special properties and might be solved easily even if the problem is NP-hard. Even if we could not exhibit special properties of the biological models, our method computes the sparsest basis of conservation laws for most curated models of the Biomodels database (see Section 5), thus validating our approach.

Some usual linear algebra algorithms can sometimes produce a sparser basis, with no guarantee it is a sparsest one. The Hermite normal form is such a technique, the (reduced or not) row echelon form is another. In the context of Z-lattices, [START_REF] Fieker | Short bases of lattices over number fields[END_REF] introduces and computes "short" (pseudo-) bases using the LLL algorithm [START_REF] Lenstra | Factoring Polynomials with Rational Coefficients[END_REF] and a variant of the Hermite normal form. In a numerical context (i.e. using floating point coefficients), there are methods to compute sparser basis (as in [START_REF] Berry | An algorithm to compute a sparse basis of the null space[END_REF] where the turnback algorithm computes a sparse and banded null basis of a large and sparse matrix).

[14, 10] present method based on numerical computations (QR decomposition, SVD,. . . ) to compute exact conservation laws of large systems. Finally, we wonder if the techniques used in the extreme pathways [START_REF] Schuster | On elementary flux modes in biochemical reaction systems at steady state[END_REF][START_REF] Schilling | Theory for the Systemic Definition of Metabolic Pathways and their use in Interpreting Metabolic Function from a Pathway-Oriented Perspective[END_REF] could be used to compute sparse conservation laws.

The paper is organized in the following way. Section 2 presents, on an example, the idea of our algorithm ComputeSparsestBasis(B), which computes a sparsest basis equivalent to B. Section 3 presents the algorithm ComputeSparsestBasis and its sub-algorithms, and Appendix A details their proofs. Section 4 details the implementation and improvements. Finally, Section 5 shows our benchmarks.

We illustrate our method on the model number 361 (BIOMD0000000361) of the BioModels database [START_REF]Biomodels database[END_REF][START_REF] Panteleev | Spatial propagation and localization of blood coagulation are regulated by intrinsic and protein c pathways, respectively[END_REF]. For clarity reasons, one renames the species in the following way :

V IIa T F → V , V IIa T F X → V X, V IIa T F Xa → V X a , T F P I → T , Xa T F P I → X a T , Xa T F P I V IIa T F → X a T V , Xa → X a .
The model contains the five following chemical reactions:

X + V ↔ V X, V X → V X a , V X a ↔ X a + V X a + T ↔ X a T , V + X a T ↔ X a T V .
By choosing the vector of species t X a T V, X a , X, V X a , V X, T, V, X a T , one can compute the stoichiometry matrix M and a basis of conservation laws (written row by row)

B =   1 1 1 1 1 0 0 1 1 0 0 1 1 0 1 0 0 -1 -1 -1 -1 1 0 0 
 by computing a basis of the nullspace of the transpose of M (denoted t M ). Thus, the matrix B represents the three conservation laws:

1. X a T V + X a + X + V X a + V X + X a T 2. X a T V + V X a + V X + V 3. -X a -X -V X a -V X + T
Our method for decreasing the number of nonzeros consists in finding a linear combination w = t vB of the rows of B such that w contains less nonzeros than one row B i of B. If one can find such a combination w, it feels natural to replace the row B i by w in order to decrease the total number of nonzeros in B. Repeating this process until no such linear combination can be found, one obtains a sparsest basis in terms of nonzeros. This approach is greedy and is justified in Section 3.

However, replacing a row of B by w should only be done if one maintains a basis. This last requirement is fulfilled by replacing the row B i of B by w only if v i = 0 (which loosely speaking means that the information in the row B i has been kept).

Consider the linear combination w = t vB with v = (α, β, γ), so one has w = (α + β, α -γ, α -γ, α + β -γ, α + β -γ, γ, β, α). The number of nonzeros of w clearly depends on the values of α, β and γ. In order to compute the number of nonzeros of w, one considers all possible cases, corresponding to the cancellation or the non cancellation of each element of w. In theory, if w has n components, one has 2 n cases to consider. For example, if we request w to have the form ( = 0, 0, 0, 0, 0, = 0, 0, = 0) one considers the following system of equations and inequations:

               α + β = 0 (column 1)
α -γ = 0 (columns 2 and 3) α + β -γ = 0 (columns 4 and 5) γ = 0 (column 6) β = 0 (column 7) α = 0 (column 8)

This last system admits for example the solution α = γ = 1 and β = 0. The corresponding linear combination w = t vB is w = B 1 + B 3 = (1, 0, 0, 0, 0, 1, 0, 1) which contains 3 nonzeros, and is thus better than the rows B 1 , B 2 and B 3 (which respectively involves 6, 4, and 5 nonzeros). Since w involves the rows B 1 and B 3 (i.e. α = 0 and β = 0), one can replace B 1 or B 3 by w. Note that replacing B 2 by w would lead to a matrix of rank 2, meaning that our basis has been lost. For example, replacing B 1 by w, one obtains an equivalent basis:

B =   1 0 0 0 0 1 0 1 1 0 0 1 1 0 1 0 0 -1 -1 -1 -1 1 0 0   .
In practice, one does not enumerate all the possible patterns of zeros and nonzeros for the vector w. Instead, one considers the columns of B one by one from left to right, and one builds systems of equations (corresponding to the zeros in w) and inequations (corresponding to the nonzeros in w). Since each column of B yields two cases, one builds a binary tree of systems of equations and inequations. By doing this, many branches are hopefully cut before all the columns of B have been treated. For example, if one tries to cancel the first five columns of w, one gets the system of equations

   α +β = 0 (1 st column of B) α -γ = 0 (2 nd , 3 rd columns of B) α +β -γ = 0 (4 th , 5 th columns of B)
. which only admits the useless solution α = β = γ = 0.

Let us continue the improvement of B . Following the same ideas as above, one can find the following linear combinations of the rows of B :

w = B 1 -B 2 -B 3 = (0, 1, 1, 0, 0, 0, -1, 1)
which has less nonzeros than B 3 . By replacing B 3 by w, one gets the basis

B =   1 0 0 0 0 1 0 1 1 0 0 1 1 0 1 0 0 1 1 0 0 0 -1 1   .
Finally, further computations would show that there does not exist a linear combination of the rows allowing to improve our basis B . In that case, one knows that our basis is a sparsest one. If one replaces the third line of B by 1 1 1 1 1 0 0 1 , then one gets three conservation laws with a total of 13 nonnegative coefficients (instead of 11 nonzero coefficients for B ). One can show that computing the minimal semi-positive P-invariants (as done in [START_REF] Soliman | Invariants and Other Structural Properties of Biochemical Models as a Constraint Satisfaction Problem[END_REF]) would retrieve these three conservation laws with 13 nonzeros.

To summarize, our method adopts a greedy approach by successively improving the initial basis (each time by only changing one row), until it reaches a sparsest basis. Each improvement consists in a binary tree exploration where each node is a system of equations and inequations.

The Algorithm ComputeSparsestBasis

Sparsest Basis of Conservation Laws

As mentioned in the introduction, we have chosen to compute sparsest basis of conservation laws. We define that notion precisely in this part.

Let B be a matrix of dimensions m × n over Q, with m ≤ n. The matrix B is called a basis matrix if B is a full rank matrix i.e. Rank(B) = m. In this paper, a basis matrix contains a basis of conservation laws written row by row. Let B and B be basis matrices of same dimensions. B and B are equivalent if and only if there exists an invertible matrix Q such that B = QB . Let M (resp. v) be a matrix (resp. vector), we denote N (M ) (resp. N (v)) the number of nonzero coefficients of M (resp. v). Let v be a vector, we denote N k (v) the number of nonzero coefficients in the first k coefficients of v. A basis matrix B is a sparsest basis matrix if and only if for any basis matrix B equivalent to B, one has N (B ) ≤ N (B).

For any basis matrix B, it is clear that there exists a sparsest matrix B equivalent to B. Indeed, consider the set of all equivalent matrices to B, and pick one matrix B in that set such that N (B ) is minimal.

A Greedy Approach

Our method follows a greedy approach. Given a basis matrix B, one looks for a vector v and an index i such that N ( t vB) < N (B i ) and v i = 0. If such v and i exists, one can decrease the number of nonzeros of B by replacing the row B i by t vB. Moreover, the rank of B does not drop since one has v i = 0. When such suitable v and i do no exist, our method stops and claims that our basis has become a sparsest one. This last claim is not obvious, since one could have fallen in a local minimum. The following theorem justifies our greedy approach. Theorem 1. A basis matrix B is not a sparsest one if and only if there exist a vector v and an index j such that N ( t vB) < N (B j ) and v j = 0.

Proof. ⇐: Taking B = B and replacing the row B j by t vB, one gets a matrix B equivalent to B such that N (B ) < N (B), which proves that B is not a sparsest one.

⇒: Assume B has dimensions m × n. There exists B equivalent to B such that N (B ) < N (B). By permuting the rows of B and the rows of B , one can assume

N (B 1 ) ≥ N (B 2 ) ≥ • • • ≥ N (B m ) and N (B 1 ) ≥ N (B 2 ) ≥ • • • ≥ N (B m ). As N (B) = m i=1 N (B i ) > m i=1 N (B i ) = N (B ), there exists an index k such that N (B k ) > N (B k ). Since B and B are equivalent, each row of B is a linear combination of rows of B. If all the m -k + 1 rows B k , B k+1 , • • • , B m
were linear combinations of the m -k rows B k+1 , . . . , B m , then B would not be a full rank matrix. Thus, there exist a vector v and indices j, l with j ≤ k ≤ l such that

B l = t vB with v j = 0. Since N (B l ) ≤ N (B k ) < N (B k ) ≤ N (B j ), one has N ( t vB) = N (B l ) < N (B j ) with v j = 0.

Description of a Task

As explained in Section 2, our method builds a (binary) tree of systems of equations and inequations. In practice, one only stores the leaves of the tree in construction. One introduces the notion of task which basically represents one leaf of the tree. In order to cut useless branches as soon as possible, one also requires a task to satisfy the extra properties LCP and IZP of Definition 1. Let v be a vector of dimension n. The notation v ≡ 0 means that ∀i ∈ 1, n , v i = 0. Definition 1. Let B be a basis matrix. Let A, Λ be matrices with m columns. Let (S) :

Ax = 0 Λx ≡ 0 be a system in the variable x. Let c be the number of rows of Λ.

Let k be the sum of the number of rows of A and Λ.

A task t = task[A, Λ, c, k],
stemming from B, is defined as follow :

• the union of the rows of A and Λ coincides with the first k columns of B (up to the order), • (A, Λ) satisfies the so-called LCP property (Linear Combination Property)

i.e. there exists a solution v of (S) with at least two nonzero coefficients, • (A, Λ, c, k) satisfies the so-called IZP property (Increase Zeros Property) i.e.

there exist a solution v of (S) and an index j such that v j = 0 and N k ( t vB) < N (B j ).

Proposition 1. Consider a task t = task[A, Λ, c, k] stemming from a basis matrix B, and the system (S), as defined in Definition 1. Consider U = {i ∈ 1, m , c < N (B i )}. Then one has the following properties:

1. A ∈ Q (k-c)×m and Λ ∈ Q c×m with 0 ≤ c ≤ k ≤ n, 2.
For each nonzero solution v of (S), N k ( t vB) = c i.e. c is the number of nonzeros in the first k coefficients of any solution of (S), 3. There exist a solution v of (S) and j ∈ U such that v j = 0 and c < N (B j ).

Proof. 1. Trivial.

2. Take a nonzero solution v of (S) and consider w = t vB. Consider a column j of B with j ≤ k. Then the transpose of this column is either a row of A or a row of Λ. If it is a row A i of A (resp. a row Λ i of Λ), then the j th coefficient of w equals zero (resp. is nonzero) since

A i v = 0 (resp. Λ i v = 0).
Consequently, the number of nonzero elements among the first k coefficients of w = t vB equals c (i.e. the number of rows of Λ). 3. It is a consequence of the IZP property and N k ( t vB) = c.

The task t 0 = task[the 0 × m matrix, the 0 × m matrix, 0, 0] is called the initial task. A task t = task[A, Λ, c, k] stemming from a basis matrix B of dimensions m × n is called a solved task if k = n.
Using Proposition 1, a solved task ensures the existence of a vector v and an index i such that N ( t vB) < N (B i ) and v i = 0, allowing the improvement of B.

The Algorithms

In this section, one presents the pseudo codes of the algorithms, and gives some hints on the way they work. The rigorous proofs of the algorithms are given in Appendix A.

Algorithm Algorithm EnhanceBasisUsingSolvedTask(t, B). It basically finds a vector v and an index i such that N ( t vB) < N (B i ) and v i = 0, which necessarily exist since t is a solved task. It then builds an improved basis B by making a copy of B and replacing the row B i by t vB.

Algorithm IsTask(t). It checks where a candidate task t is indeed a task, by checking if t satisfies the LCP and IZP properties. The goal of this function is to detect as soon as possible useless tasks (i.e. tasks that will not help improving our basis).

Algorithm NextVector(u). The goal of NextVector is to iterate the p-tuples of Z p . This is needed in Algorithm EnhanceBasisUsingSolvedTask to obtain a solution v of the system (S), which is composed of equations (i.e. Ax = 0) and inequations (i.e. Λx ≡ 0). Indeed, the only way we have found to obtain solutions of Ax = 0 also satisfying Λx ≡ 0 consists in iterating some solutions of Ax = 0 until Λx ≡ 0 is satisfied.

Input: an integer vector u of dimension p Output: the vector following u for some fixed ordering on Z p Algorithm 6: NextVector(u)

4 Complexity and Improvements

Complexity

The bottleneck of ComputeSparsestBasis is located in BasisToSolvedTask. Indeed, the number of while loops performed in BasisToSolvedTask can be close to 2 n+1 -1 in the worst case (i.e. when the binary tree is almost completely explored). It is the only place where an exponential complexity occurs, since all other operations rely on linear algebra operations (such as computing a nullspace, a RREF, . . . ). However, many branches are cut thanks to the line 11 in the algorithm IsTask. The sparser the matrix B is, the more branches are cut. Indeed, let us denote d = max{N (B i ), i ∈ 1, m }. Suppose that, in our binary tree, the left (resp. right) child corresponds to adding an equation (resp. an inequation). If the number of inequations c of some task is greater than or equal to d, the set U at line 10 in IsTask is empty, thus IsTask returns false. This implies that only the branches starting from the root and going through the right children at most d times will be explored. The number of processed nodes at depth k is equal to n + 1 i + 1 ≤ 2(n + 1) d+1 which can be much smaller than 2 n+1 -1 (for example when d is much smaller than n). Experimentally, we observed that models with small values d were easily solved. Finally, the number of calls to EnhanceBasis is bounded by the number of nonzeros of the initial basis B, which is bounded by nd (since the number of nonzeros decreases at least by 1 at each call of EnhanceBasis, except for the last call).

Implementation

We chose to implement algorithms given in Section 3 using the Computer Algebra software Maple, which natively handles long integers and contains many linear algebra routines with exact coefficients. With no surprise, those algorithms can be improved because many useless computations are performed. For example, many useless rank computations are done in Algorithm IsTask. The next section describes the improvements of the algorithms given in Section 3.

Improvements

Computation and Choice of the Solved Task Algorithm BasisToSolvedTask stops when it first encounters a solved task. This solved task may change if one pushes t 2 before t 1 in the Algorithm BasisToSolvedTask. This change has no real impact since it speeds up some examples, and slows down others.

We have experimented another strategy consisting in computing the set of all the solved tasks stemming from B instead of stopping at the first solved task. Once this set is computed, one can choose the solved task that leads to the biggest decrease of the number of nonzeros, and only keep solved tasks with Rank(A) = m -1 if one encounters a solved task with Rank(A) = m -1 at some point. Indeed, solved tasks with Rank(A) = m -1 lead to easy computations in Algorithm EnhanceBasisUsingSolvedTask since p = 1.

It is not clear whether the strategy of computing all the solved tasks is better or worst than stopping on the first solved task. Indeed, searching all solved tasks is obviously more costly, but choosing a suitable solved task might decrease the number of subsequent calls to Algorithm EnhanceBasis.

Using Reduced Row Echelon Forms It is possible to request more properties for a task task[A, Λ, c, k]. For example, one can request A to be in reduced row echelon form with no zero rows. Moreover, one can request the rows of Λ to be reduced w.r.t. the matrix A in the following sense: a row b is reduced w.r.t. A if the row contains a zero at the location of the pivots of A (reducing a row Λ i by A can be done by subtracting multiple of rows of A to Λ i to get zeros at the positions of the pivots of A).

Those requirements have several advantages, especially in Algorithms 4 and 5. In Algorithm 4, the computation of Rank(A) is immediate since it equals the number of rows of A. Moreover the condition Rank A Λ i = Rank(A) + 1 can be checked immediately: indeed, since Λ i is reduced w.r.t. A, the condition is true if and only if the row Λ i is not the zero row. In Algorithm 5, the computation of the matrix K is immediate and can be done by simply rearranging the entries of A in a new matrix K.

Finally, since the rows Λ i are reduced w.r.t. A, it is easier to detect that two rows Λ i and Λ j are equal modulo a linear combination of lines of A. Indeed, if this is the case, both lines are necessarily proportional (and then one can discard one of them).

Benchmarks

We have tested our methods on some models taken from the BioModels database [START_REF]Biomodels database[END_REF] (accessed on June 16th, 2014). Among the curated models, we have selected all the models involving only one compartment, with rational and integer stoichiometric coefficients. After rejecting models without conservation laws, we ended up with 214 models.

After computing a basis of linear conservation laws for each of the 214 models, our method detected that 141 bases were already sparsest ones. In the rest of the section, one only considers the remaining 73 models.

In this section, the CSB version denotes the non-improved version of the algorithm ComputeSparsestBasis, as described in Section 3, and the CSB' version denotes the version based on RREF computations as described in Section 4.3. Timings were measured on a Pentium Xeon 3.40GHz with 32Gb of memory. Checking that the 141 bases are indeed sparsest ones takes at most 1s for each model. The improved CSB' version takes less than 3000s for 69 models out of the 73 models. The remaining 4 models involve heavier computations: model 332 takes around 4000s, model 175 takes around 1 day, computation of models 205 and 457 were stopped after 2 days. These 4 models involve between 15 and 40 conservation laws (i.e. B has between 15 and 40 rows), and between 50 and 200 species (i.e. B has between 50 and 200 columns). As shown in [START_REF] Thomas | The null space problem i. complexity[END_REF], finding a sparsest basis for the null space of a matrix is NP-hard, so it is not surprising that some models are challenging.

The CSB' version is faster than the basic version, usually by a factor of at least 10. Figure 1 shows some timings for a sample of the 69 models. Note that the timings can be different between models with bases of similar sizes (like models 153 and 152).

For each basis B, one can define the ratio x = N (B ) N (B) where B is a sparsest basis equivalent to B. For a non sparsest basis B, this ratio satisfies 0 < x < 1. Figure 2 represents the frequency of the 69 models, plus the models 332 and 175 involving heavier computations, as a function of the ratio x.

Appendix B shows a comparison between the number of nonzeros obtained after some usual linear algebra methods and our algorithm.

Conclusion

We plan to implement our method in Python, for different reasons: it is a free software, it handles large integers natively, and it is easy to interface with other software. One could also implements our method in C (by using the multiprecision arithmetic library GMP [START_REF] Granlund | GNU MP: The GNU Multiple Precision Arithmetic Library[END_REF]) for performance reasons. Looking back to the problem of getting "good" sets of conservation laws we discussed in Section 1, it is not always possible to have a basis with the less possible negative values and a sparsest basis at the same time. There probably exists a compromise between these two properties, which is left for further work.

2. K j = 0.

Proof. (1) ⇒ (2): There exists a row A i = 0 • • • 0 1 0 • • • 0 where the 1 is at position j, for some i in 1, r . Thus, any solution v (in particular all elements of the basis K) of A v = 0 must satisfy v j = 0, hence K j = 0.

(2) ⇒ (1): From A K = 0, one has t K t A = 0. Since K j = 0, the column j of t K is zero, which implies that the canonical vector e j belongs to Ker( t K) using Proposition 2. Hence, the row l = t e j = 0 • • • 0 1 0 • • • 0 , where the 1 is at position j, is a linear combination of the rows of A . Since A is in reduced row echelon form, if the combination l of rows of A involved strictly more than one row of A , l would at least involve two nonzero coefficients (corresponding to the pivots). Thus, l is a row of A and A is row-unit of index j.

Halting Trivial. (S) has nonzero solutions according to Theorem 2. 2. If the second condition (line 6) is true, since (S) admits nonzero solutions and thanks to Theorem 3, (S) does not have solutions with at least two nonzero coefficients, so t does not satisfy LCP, and one returns false. Otherwise, t satisfies LCP. 3. If the third condition (line 11) is true, then for any j ∈ U, one has K j = 0 thanks to Theorem 4. Consequently, any solution v of (S) satisfies v j = 0 for all j ∈ U. Therefore, t does not satisfy IZP and one returns false.

Otherwise, there exists a j ∈ U such that K j = 0 thanks to Theorem 4, so there exists a column w of K such that w j = 0. It may happen that w j is not solution of (S). In that case, consider a nonzero solution u of (S) with at least two nonzero coefficients. By a topology argument, the vector ū = u + ε 1 w is also solution of (S) for any rational ε 1 satisfying |ε 1 | < ε (where ε is a small rational). Finally, one can choose a suitable ε 1 to obtain the condition ūj = 0, and (S) satisfies IZP. One then returns true.

B Comparison with some Matrix Algorithms

We consider a sample of 10 models taken from the 61 models where computations end in less than 3000s. We present here the number of nonzeros of the initial bases and of the bases obtained after using our algorithm (CSB), the Reduced Row Echelon Form (RREF), the LLL algorithm [START_REF] Lenstra | Factoring Polynomials with Rational Coefficients[END_REF] and the Hermite Normal Form (HNF). When comparing our algorithm with the usual linear algebra algorithms we used, one sees that the sparsest bases are not always reached by these linear algebra algorithms and that they sometimes make it worst.

2 B ← B ; 3 5 if p = 1 then 6 v 8 else 9 i 10 while i = 0 do 11 u

 2356891011 Input: a solved task t = task[A, Λ, c, n], stemming from B of size m × n Output: B a basis matrix equivalent to B such that B and B only differ by one row and N (B ) < N (B) 1 begin Compute a basis of Ker(A) and store it columnwise in the matrix K of dimensions m × p, where p = m -Rank(A) ; 4 Compute U = {i ∈ 1, m , c < N (Bi)} ; ← the unique column of K ; 7 Choose i ∈ U such that vi = 0 ; ← 0 ; u ← 0 ;// u is the zero vector of dimension p ← NextVector(u) ;

3 if 4 return false ; 5 // 7 return false ; 8 // 9 Compute

 345789 ∃i ∈ 1, c , Rank A Λi = Rank(A) or (Rank(A) = m) then One has ∀i, Rank A Λi = Rank(A) + 1 and (Rank(A) ≤ m -1) 6 if (Rank(A) = m -1) and A contains at least one zero column then One has (Rank(A) ≤ m -2) or A does not have any zero column ARREF , the RREF form of A ; 10 Compute U = {i ∈ 1, m , c < N (Bi)} ; 11 if ∀j ∈ U , ARREF is row-unit of index j then 12

Algorithm 5 :

 5 IsTask(t, B)

Fig. 1 .

 1 Fig. 1. Timing comparison between the basic version and the improved version (with B basis matrices of the models and d = max{N (Bi), i ∈ 1, m } as described Section 4.1)

Fig. 2 .

 2 Fig. 2. Number of bases given some proportion x = N (B ) N (B)

Correction 1 .

 1 If the first condition (line 3) is true, then (S) has no nonzero solutions in Q m , so t does not satisfy condition LCP, and one returns false. Otherwise, one has (∀i ∈ 1, c , Rank A Λ i = Rank(A) + 1) and (Rank(A) ≤ m -1) so

Fig. 3 .

 3 Fig. 3. Number of nonzeros after some matrix algorithms

  ComputeSparsestBasis(B). It is the main algorithm. It takes a basis B as an input and returns a sparsest basis B equivalent to B. It relies on the algorithm EnhanceBasis(B) which either returns an equivalent basis B with N (B ) < N (B) or proves that B was a sparsest basis. Thus, ComputeSparsestBasis(B) iterates calls to EnhanceBasis(B) until the basis is a sparsest one. : B a basis matrix of dimensions m × n Output: B , a sparsest basis matrix equivalent to B 1 begin It relies on the algorithms BasisToSolvedTask and EnhanceBasisUsingSolvedTask. The algorithm BasisToSolvedTask(B) builds a solved task stemming from B if it exists, or returns the empty set if no such solved task exists. If such a solved task can be computed, Algorithm EnhanceBasisUsingSolvedTask is used to improve the basis B. Algorithm BasisToSolvedTask(B). It looks for a solved task by exploring a binary tree. It makes use of a stack, initially filled with the initial task. At each step of the loop, a task t (where k columns of B have been processed) is popped and two new candidate tasks t 1 and t 2 are built by processing the column k +1 of B. These candidates are pushed onto the stack if they are actually tasks (which is checked by Algorithm IsTask).

		Input: a basis matrix B
		Output: a solved task t, stemming from B, or ∅ if no such solved task exists
		(i.e. B is a sparsest basis)
	1 begin
	2	Let St be an empty stack ;
	9	InputB) A ← A w ; Λ ← Λ w ;
		otherwise
	1 begin
	2	t ← BasisToSolvedTask(B) ;
	3	if t = ∅ then
	4	B ← EnhanceBasisUsingSolvedTask(t, B) ;
	5	return true, B ;
	6	else
	7	return false, B ;
		Algorithm 2: EnhanceBasis(B)

2 B ← B ; a ← true ; 3 while a do 4 a, B ← EnhanceBasis(B ) ; 5 return B ; Algorithm 1: ComputeSparsestBasis(B)

Algorithm EnhanceBasis(B). Input: B a basis matrix of dimensions m × n Output: One of the two cases: false and B if B is a sparsest basis matrix ; true and a basis matrix B equivalent to B such that N (B ) < N (3 Push the initial task t0 onto St ; 4 while St = ∅ do 5 Pop a task t = task[A, Λ, c, k] from St ; 6 if k < n then 7 // The task t is not solved 8 Let w be the transpose of the (k + 1) th column of B ; 10 t1 ← [A , Λ, c, k + 1] ; // t1 may be a task 11 if IsTask(t1, B) then Push t1 onto St ; 12 ; 13 t2 ← [A, Λ , c + 1, k + 1] ; // t2 may be a task 14 if IsTask(t2, B) then Push t2 onto St ; 15 ; 16 else 17 return t ; 18 return ∅ ; Algorithm 3: BasisToSolvedTask(B)

  LCP (resp. IZP) is true if the tests lines 6 (resp. 11) and 3 are false

12 v ← Ku ; 13 if Λv ≡ 0 then 14 // v is a nonzero solution of (S) 15 Choose i ∈ U such that vi = 0 if it exists ; 16 B i ← t vB ; 17 Multiply B i by the LCM of the denominators of the elements of B i ; 18 Divide B i by the GCD of the elements of B i ; 19 return B ; Algorithm 4: EnhanceBasisUsingSolvedTask(t, B) Input: t = [A, Λ, c, k], satisfying all conditions of a task stemming from a basis matrix B of dimensions m × n, except LCP and IZP properties Output: true if t satisfies LCP and IZP (i.e. t is a task), false otherwise 1 begin 2 //

Acknowledgements. We would like to thank the reviewers for their helpful comments and Sylvain Soliman for his help on the Nicotine software (http: //contraintes.inria.fr/ ~soliman/nicotine.html).

A Proof of Algorithms

Each subsection of this section proves an algorithm, by proving it stops and returns the correct result. When needed, some mathematical properties are proved at the beginning of the subsections.

A.1 ComputeSparsestBasis(B)

Halting The algorithm halts because EnhanceBasis returns either false, B or true, B and this last case cannot happen indefinitely as the number of nonzeros in B is strictly decreasing and bounded by 0.

Correction The algorithm halts when a = false, i.e. when EnhanceBasis detects that B is a sparsest basis.

A.2 EnhanceBasis(B)

Halting Trivial.

Correction BasisToSolvedTask returns either a solved task or ∅ if no such solved task exists. If t = ∅, EnhanceBasisUsingSolvedTask computes a new basis matrix B with N (B ) < N (B), otherwise no solved task exists, which proves that B is a sparsest basis.

A.3 EnhanceBasisUsingSolvedTask(t, B)

Lemma 1 (Finite union of vector subspaces). Let K be an infinite field. If E is a K-vector subspace, then every finite union of proper subspaces of E is strictly included in E.

Proof. See [START_REF] Roman | Advanced Linear Algebra[END_REF].

Theorem 2 (Characteristic theorem of solutions of (S)).

matrix representing a basis of Ker(A) stored columnwise, and the system (S) :

The following assertions are equivalent:

Proof.

(1) ⇒ (3): Take a nonzero solution v in Q m of (S). Thus, there exists u in Q q \ {0} such that v = Ku (since the columns of K are a basis of Ker(A)).

Taking u = λu (with a suitable integer λ such that u belongs to Z q \ {0}), one has Ku = λv. Since λv is also a nonzero solution of (S), one has ΛKu ≡ 0.

(3) ⇒ (1): Take v = Ku. Then Av = 0 and Λv ≡ 0, so v is a nonzero solution of (S).

(1) ⇒ (2): Consider a nonzero solution v of (S). Since v is nonzero and satisfies Av = 0, one has Ker(A) = {0} thus Rank(A) ≤ m -1 (according to the rank-nullity theorem). Let us suppose that Rank A Λ i = Rank(A) for some i.

Then, Λ i is a linear combination of rows of A so Λ i v = 0, hence v is not a solution of (S), contradiction. We conclude that ∀i ∈ 1, c , Rank

Halting If p = 1, it is trivial that the algorithm halts. If p = 1, one needs to check that the while loop halts:

• As t is a solved task, there exists a nonzero solution v in Q m of (S) and j ∈ U such that v j = 0 and c < N (B j ) according to Proposition 1. • Using the third point of Theorem 2, there exists a vector u in Z q such that the vector v = Ku is a nonzero solution of (S). Hence, the while loop will eventually reach such a u since NextVector enumerates all elements of Z q .

Correction At line 2, B ← B. At the end of the algorithm, B i is modified and contains more zeros than

Finally, B is also a basis matrix since B i = t vB with v i = 0.

A.4 NextVector

There is some freedom in coding this algorithm (which is not given in this paper), in particular the ordering in which the p-tuples are output. However, to ensure the termination of Algorithm EnhanceBasisUsingSolvedTask, one requires the Algorithm NextVector(u) to iterate all the p-tuples of Z p . This can be achieved for example by starting from the zero tuple, and enumerating the tuples of Z p by increasing 1-norm (where the 1-norm of a tuple is the sum of the absolute values), and by lexicographic order for tuples of the same 1-norm. For example, when p = 2, the p-tuples can be enumerated in the following way: (0, 0), (-1, 0), (0, -1), (0, 1), (1, 0), (-2, 0), (-1, -1), (-1, 1), (0, -2), (0, 2), . . . One could also rely on a random number generator, provided it has the property to eventually generate any p-tuple with a nonzero probability (in order to ensure the halting of Algorithm EnhanceBasisUsingSolvedTask).

A.5 BasisToSolvedTask(B)

Halting Consider a current task t. The algorithm stops if t is solved (i.e. k = n). Otherwise, t is not a solved task and generates new tasks, obtained from t, with k + 1 columns processed (instead of k). This last case cannot occur indefinitely.

Correction Consider a non solved task t inside the while loop. One creates the object t 1 (resp. t 2 ) corresponding to the cancellation (resp. the non cancellation) of the coefficient number k + 1 of the linear combinations of the lines of B. The objects t 1 and t 2 are possibly discarded thanks to the function IsTask if they are not tasks (i.e. if they cannot be used to increase the number of zeros). Consequently, all cases are considered and the function will return ∅ if and only if there does not exist any solved task stemming from B.

A.6 IsTask(t, B) Proposition 2. Let M be a matrix. Then M has a zero column if and only if Ker(M ) contains at least one vector with exactly one nonzero coefficient.

Proof. Trivial

Corollary 1. Let M be a matrix. Then M does not have any zero column if and only if Ker(M ) \ {0} only contains vectors with at least two nonzero coefficients. Proposition 3. Take A ∈ Q r×m , Λ ∈ Q c×m and the system (S) : Definition 2. A matrix A is row-unit of index j if there exists a row of A with only one nonzero coefficient, which is at position j. Theorem 4. Consider a matrix A in Q r×m under reduced row echelon form. Let K be a matrix containing a basis of Ker(A ) stored columnwise. Suppose that Ker(A ) = {0}. For any index j, the following assertions are equivalent: