Mehdi Maamar
email: maamar.mehdi@edu.univ-oran1.dz

Nadjib Lazaar
email: nadjib.lazaar@lirmm.fr

Samir Loudni
email: samir.loudni@unicaen.fr

Yahia Lebbah
email: lebbah.yahia@univ-oran.dz

Fault localization using itemset mining under constraints

Keywords: Fault localization, Itemset mining, Constraint programming, Test case coverage 1 Introduction

We introduce in this paper an itemset mining approach to tackle the fault localization problem, which is one of the most difficult processes in software debugging. We formalize the problem of fault localisation as finding the k best patterns satisfying a set of constraints modelling the most suspicious statements. We use a Constraint Programming (CP) approach to model and to solve our itemset based fault localization problem. Our approach consists of two steps: i) mining top-k suspicious suites of statements; ii) fault localization by processing top-k patterns. Experiments performed on standard benchmark programs show that our approach enables to propose a more precise localization than a standard approach.

bugs are [START_REF] Vessey | Expertise in debugging computer programs: A process analysis[END_REF], this is known as fault localization problem. Nonetheless, locating a fault is still an extremely time-consuming and tedious task. Over the last decade, several automated techniques have been proposed to tackle this problem.

Most of automated techniques for fault localization compare two kinds of execution traces, namely the passed and failed executions [START_REF] Jones | Visualization of test information to assist fault localization[END_REF], such as Tarantula [START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF] which is one of the most popular fault localization technique. These methods are based on a scoring function to evaluate the suspiciousness of each statement in the program by exploiting the occurrences of the considered statements in passing and failing test cases. Since, the statements with high scores correlate highly with faults and that brings us to a total order on statements from highly suspicious to guiltless statements. It is important to stress that this technique does not differentiate between two failing (resp. passing) test cases, and consequently it ignores the dependencies between statements that can help us to locate the fault.

Other techniques take in account the cause effect chains with a dependence analysis by using, for instance, program slicing [START_REF] Agrawal | Debugging with dynamic slicing and backtracking[END_REF]. The disadvantage here is that fault can be located in a quite large slice (static slicing) and/or can be time/space consuming (dynamic slicing).

Recently, the problem of fault localization was abstracted as a data mining problem. Cellier andal. (Cellier et al, 2009, 2011) propose a combination of association rules and Formal Concept Analysis (FCA) to assist in fault localization. The proposed methodology tries to identify rules between statement execution and corresponding test case failure. The extracted rules are then ordered as a lattice and explored bottom up to detect the fault.

In the data mining community, many approaches have promoted the use of constraints to focus on the most promising knowledge according to a potential interest given by the final user. As the process usually produces a large number of patterns, a large effort is made to a better understanding of the fragmented information conveyed by the patterns and to produce pattern sets i.e. sets of patterns satisfying properties on the whole set of patterns [START_REF] De Raedt | Constraint-based pattern set mining[END_REF][START_REF] Rojas | Mining (soft-) skypatterns using dynamic CSP[END_REF]. Discovering top-k patterns (i.e. the k best patterns according to a score function) is a recent trend in constraintbased data mining to produce useful pattern sets [START_REF] Crémilleux | Discovering knowledge from local patterns with global constraints[END_REF].

In this spirit of this promising avenue, we propose in this paper an itemset mining approach to tackle the fault localization problem. We formalize the problem of fault localization as finding the k best patterns satisfying a set of constraints modeling the most suspicious statements. We use the test cases coverage of a program collected during the testing phase for finding the location of program fault. Our approach, which benefits from the recent progress on cross-fertilization between data mining and Constraint Programming [START_REF] Guns | Itemset mining: A constraint programming perspective[END_REF][START_REF] Khiari | Constraint programming for mining n-ary patterns[END_REF], is achieved in two steps:

1. Mining top-k suspicious patterns (i.e., set of statements) according to dominance relations and using constraints dynamically posted, 2. Ranking the statements by processing the top-k patterns in an ad-hoc ranking algorithm to locate the fault.

Our localization approach is based on two working hypothesis. The first hypothesis is usually called the competent developer hypothesis [START_REF] Demillo | Hints on test data selection: Help for the practicing programmer[END_REF] . It states that, even some faults are introduced, the resulting program would certainly fulfill almost all of its specifications. In other words, the resulting program may contain faults, but it will certainly tackle the problem it has been designed for. The second requirement on which our approach is based is the single fault hypothesis, i.e., there is only one faulty statement in the program. This requirement might appear as being restrictive but it has been shown that complex faults usually result from the coupling of single faults (a.k.a. the coupling effect [START_REF] Jones | Visualization of test information to assist fault localization[END_REF]).

Experiments performed on several benchmark of single fault programs (Siemens suite) show that our approach enables to propose a more precise localization as compared to the most popular fault-localization technique Tarantula [START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF]. To the best of our knowledge, our proposal is the first data mining approach that exploits pattern sets to fault localization.

This paper is organized as follows. Section 2 sketches definitions and presents the context. Section 3 gives a background on pattern discovery and describes the CP modeling for itemset mining. Section 4 presents the detail of our approach for fault localization. Section 5 presents an illustrative example of fault localization using our approach. Section 6 reports experimental results and a complete comparison with Tarantula. Section 7 presents the related faultlocalization methods in the area of data mining and using failing/passing executions. Section 8 concludes the paper and gives some directions for the future works.

Background and Motivation

This section presents background knowledge about the problem of fault localization and constraint satisfaction problems.

Fault localization problem

In software engineering, a failure is a deviation between expected and actual result. An error is the part of the program that is liable to lead to a subsequent failure. Finally, a fault in the general sense is the adjudged or hypothesized cause of an error (Laprie et al, 1992). The purpose of fault localization is to pinpoint the root cause of observed symptoms under test cases.

Given a faulty program P having n lines, labeled L = {e 1 , e 2 , ..., e n }. For instance, for the program "Character counter" given in Figure 1, we have L = {e 1 , e 2 , ..., e 10 }.

(c = s[i++]) { 1 1 1 1 1 1 1 1 e2: if('A'<=c && 'Z'>=c) 1 1 1 1 1 1 0 1 e3: let += 2; //-fault - 1 1 1 1 1 1 0 0 e4: else if ( 'a'<=c && 'z'>=c ) 1 1 1 1 1 0 0 1 e5: let += 1; 1 1 0 0 1 0 0 0 e6: else if ( '0'<=c && '9'>=c ) 1 1 1 1 0 0 0 1 e7: dig += 1; 0 1 0 1 0 0 0 0 e8: else if (isprint (c)) 1 0 1 0 0 0 0 1 e9: other += 1; 1 0 1 0 0 0 0 1 e10: printf("%d %d %d\n", let, dig, other);} 1 1 1 1 1 1 1 1 Passing/Failing F F F F F F P P
Fig. 1: Example of program and its associated transactional dataset. Definition 4 (Test case coverage) Given a test case tc i and a program P , the set of executed (at least once) statements of P with tc i is a test case coverage I i = (I i,1 , ..., I i,n ) where I i,j = 1 if the j th statement is executed, 0 otherwise.

A test case coverage indicates which parts of the program are active during a specific execution. For instance, the test case tc 4 in Figure 1 covers the statements e 1 , e 2 , e 3 , e 4 , e 6 , e 7 , e 10 . The according test case coverage is then I 4 = (1, 1, 1, 1, 0, 1, 1, 0, 0, 1).

Constraint Satisfaction Problem

A CSP consists of a finite set of variables X = {x 1 , . . . , x n } with finite domains D = {D 1 , . . . , D n } such that each D i is the set of values that can be assigned to x i , and a finite set of constraints C. Each constraint C(Y ) ∈ C express a relation over a subset Y of variables X. The objective is to find an assignment (x i = d i ) with d i ∈ D i for i = 1, . . . , n, such that all constraints are satisfied.

Example 1 Let be the following CSP:

X = {x 1 , x 2 , x 3 } D = {D 1 , D 2 , D 3 } where, D 1 = {1, 2}, D 2 = {0, 1, 2, 3}, D 3 = {2, 3} C = {C 1 (x 1 , x 2 ), C 2 (x 1 , x 2 , x 3 ), C 3 (x 1 )} where, C 1 (x 1 , x 2 ) : x 1 > x 2 C 2 (x 1 , x 2 , x 3 ) : x 1 + x 2 = x 3 C 3 (x 1 ) : x 1 = 0
Here, the current CSP admits two solutions :

(x 1 = 2, x 2 = 0, x 3 = 2) and (x 1 = 2, x 2 = 1, x 3 = 3).
In Constraint Programming (see [START_REF] Rossi | Handbook of Constraint Programming (Foundations of Artificial Intelligence)[END_REF]), the resolution process consists of iteratively interleaving search phases and propagation phases. The search phase essentially consists of enumerating all possible variable-value combinations, until we find a solution or prove that none exists. It is generally performed on a tree-like structure. In order to avoid the systematic generation of all the combinations and reduce the search space, the propagation phase shrinks the search space: each constraint propagation algorithm (also called propagator) removes values that a priori cannot be part of a solution w.r.t. the partial assignment built so far. The removal of inconsistent domain values is called filtering. If all inconsistent values are removed from the domains with respect to a constraint C, we say that C is domain consistent.

Consider for example the constraint C 1 (x 1 , x 2 ) : x 1 > x 2 . The propagator for this constraint will remove values 2 and 3 from D 2 . The repeated application of propagators can lead to a successive reduction of domains until reaching a fixed point where no value can be pruned. At this point, the search assigns a variable one of its values. Whenever the domain of one of the variables becomes empty, the search backtracks to explore alternatives.

Constraint programming provides a very expressive type of constraints. One can denote the predefined constraints (i.e., arithmetic constraints), constraints given in extension (list of allowed/forbidden combinations of values), and logical combination of constraints. Another kind of constraints are reified constraints, also known as meta constraints. A reified constraint b ↔ c involves a boolean variable b and a constraint c, and it is equivalent to (b = 1 ∧ c) ∨ (b = 0∧¬c). This kind of constraints is useful to express propositional formulas over constraints, formulas that are common in data mining (see Section 3).

Frequent Itemset Mining

This section gives a brief overview of the CP approach for itemset mining [START_REF] De Raedt | Constraint programming for itemset mining[END_REF][START_REF] Guns | Itemset mining: A constraint programming perspective[END_REF].

Context and definitions

Let I be a set of distinct literals called items and T = {1, ..., m} a set of transaction identifiers. An itemset (or pattern) is a non-null subset of I. The Table 1

: Transactional dataset D Trans. Items t 1 B E F t 2 B C D t 3 A E F t 4 A B C D E t 5 B C D E t 6 B C D E F t 7 A B C D E F Table 2: Binary matrix Trans. A B C D E F t 1 0 1 0 0 1 1 t 2 0 1 1 1 0 0 t 3 1 0 0 0 1 1 t 4 1 1 1 1 1 0 t 5 0 1 1 1 1 0 t 6 0 1 1 1 1 1 t 7 1 1 1 1 1 1
language of itemsets corresponds to L I = 2 I \∅. A transactional dataset is a set D ⊆ I×T . Table 1 presents a transactional dataset D where each transaction t i gathers articles described by items denoted A,. . . ,F . The traditional example is a supermarket dataset in which each transaction corresponds to a customer and every item in the transaction to a product bought by the customer.

Definition 5 (Coverage and Frequency) The coverage of an itemset x is the set of all identifiers of transactions in which x occurs: cover D (x) = {t ∈ T |∀i ∈ x, (i, t) ∈ D}. The frequency of an itemset x is the size of its cover:

freq D (x) = |cover D (x)|.
Example 2 Consider the transactional dataset in Table 1. We have for x = BEF that cover D (x) = {t 1 , t 6 , t 7 } and freq D (x) = 3.

Constraint-based pattern mining aims at extracting all patterns x of L I satisfying a query q(x) (conjunction of constraints), which usually define what we call a theory [START_REF] Mannila | Levelwise search and borders of theories in knowledge discovery[END_REF]

: T h(q) = {x ∈ L I | q(x) is true}.
A common example is the frequency measure leading to the minimal frequency constraint. The latter provides patterns x having a number of occurrences in the dataset exceeding a given minimal threshold min f r : freq D (x) ≥ min f r . Another usual constraint is the size constraint which constrains the number of items of a pattern x.

Example 3 Let us consider the following query q(x) = freq D (x) ≥ 5 ∧ size(x) = 2. It addresses all frequent patterns (min f r = 5), having a size equal to 2. With the running example in Table 1, we get four solutions : BE, BC, BD and CD.

In many applications, it appears highly appropriate to look for contrasts between subsets of transactions, such as passing and failing test cases in fault localization (see Section 4). The growth rate is a well-used contrast measure [START_REF] Novak | Supervised descriptive rule discovery: A unifying survey of contrast set, emerging pattern and subgroup mining[END_REF]. Let D be a dataset partitioned into two subsets D 1 and D 2 : Definition 6 (Growth rate) The growth rate of a pattern x from D 2 to D 1 is:

m gr (x) = |D 2 | × f req D1 (x) |D 1 | × f req D2 (x)
Emerging Patterns are those having the growth rate greater than some given threshold. Jumping Emerging Patterns are those which do not occur in D 2 . They are at the core of a useful knowledge in many applications involving classification features such as the discovery of suspicious statements in the program (see Section 5).

The collection of patterns contains redundancy w.r.t. measures. Given a measure m, two patterns x i and x j are said to be equivalent if m(x i ) = m(x j ). A set of equivalent patterns forms an equivalent class w.r.t. m. The largest element w.r.t. the set inclusion of an equivalence class is called a closed pattern.

Definition 7 (Closed pattern) A pattern x i ∈ L I is closed w.r.t. a measure m iff ∀ x j ∈ L I , x j x i ⇒ m(x j ) = m(x i ).
The set of closed patterns is a compact representation of the patterns (i.e we can derive all the patterns with their exact value for m from the closed ones).

Example 4 Consider the frequency measure (i.e. m = f req). In our running example, if we impose that freq D (x) ≥ 5, we get 9 frequent patterns which are summarized by four equivalence classes (and thus 4 closed frequent patterns). For instance, BCD 5 1 is a closed pattern for BC 5 , BD 5 , CD 5 , C 5 and D 5 . The three other equivalence classes are: B 6 , E 6 and BE 5 .

Moreover, the user is often interested in discovering richer patterns satisfying properties involving several local patterns. These patterns define pattern sets [START_REF] De Raedt | Constraint-based pattern set mining[END_REF] or n-ary patterns [START_REF] Khiari | Constraint programming for mining n-ary patterns[END_REF]. The approach that we present in this paper is able to deal with pattern sets such as the top-k patterns.

CP model for the itemset mining

As defined in Section 3.1, let D be a dataset where I is the set of its n items and T the set of its m transactions. The set of items can be indexed by consecutive integers and thus can be referenced by their indexes; consequently, without loss of generality, the set of items I is supposed to be a set of n integers I = {1, ..., n}. The transactional dataset D can be represented with a 0/1 (m,n) transactional boolean matrix d, such that

∀t ∈ T , ∀i ∈ I, (d t,i = 1) ↔ (i ∈ t).
1 Value between . indicates the frequency of a pattern.

Table 2 illustrates the transactional binary matrix of the transactional dataset given in Table 1, where 1 (resp. 0) means that an article is (resp. is not) in a transaction. In [START_REF] De Raedt | Constraint programming for itemset mining[END_REF], the authors model an unknown pattern M ⊆ I and its associated dataset T by introducing two sets of boolean variables:

-item variables {M 1 , M 2 , ..., M n } where (M i = 1) iff (i ∈ M ), -transaction variables {T 1 , T 2 , ..., T m } where (T t = 1) iff (M ⊆ t).
The relationship between M and T is modeled by reified constraints stating that, for each transaction t, (T t = 1) iff M is a subset of t:

∀t ∈ T : (T t = 1) ↔ i∈I M i × (1 -d t,i ) = 0 (1)
Using the boolean encoding, it is worth noting that some measures are easy to encode: freq D (M ) = t∈T T t and size(M ) = i∈I M i . So, the minimal frequency constraint freq D (M ) ≥ min f r (where min f r is a threshold) is encoded by the constraint t∈T T t ≥ min f r . In the same way, the minimal size constraint size(M ) ≥ α (where α is a threshold) is encoded by the constraint

i∈I M i ≥ α.
Finally, the closedness constraint closed f req (M ) ensures that a pattern has no superset with the same frequency; it is encoded in [START_REF] Guns | Itemset mining: A constraint programming perspective[END_REF] using equation ( 1) as follows:

closed f req (M ) ≡ ∀i ∈ I : (M i = 1) ↔ t∈T T t × (1 -d t,i ) = 0
(2)

Fault Localization by Itemset Mining

In this section, we first give our modeling of the fault localization problem as an itemset mining under constraints. Then, we formalize the problem of fault localization as finding the k most suspicious patterns and we detail our algorithm for mining top-k suspicious patterns. Finally, we describe how to exploit the top-k patterns to return at the end an accurate fault localization.

Modeling the fault localization as a constrained itemset mining

Let L = {e 1 , ..., e n } be a set of indexed statements composing a program P and T = {tc 1 , ..., tc m } a set of test cases. The transactional dataset D is defined as follows: (i) each statement of L corresponds to an item in I, (ii) the coverage of each test case tc i forms a transaction in T . Moreover, to look for contrasts between subsets of transactions, T is partitioned into two disjoint subsets T + and T -. T + (resp. T -) denotes the set of coverage of positive (resp. negative) test cases.

Let d be the 0/1 (m,n) matrix representing the dataset D. So, ∀t ∈ T , ∀i ∈ I, (d t,i = 1) if and only if the statement i is executed (at least once) by the test case t. Figure 1 shows the transactional dataset associated to the program Character counter. For instance, the coverage of the test case t 5 is I 5 = (1,1,1,1,1,0,0,0,0,1). As t 5 fails, thus

I 5 ∈ T -.
Let M be the unknown suspicious pattern we are looking for. As detailed in Section 3.2, we introduce two sets of boolean variables: item variables {M 1 , M 2 , ..., M n } representing statements, and transaction variables {T 1 , T 2 , ..., T m } representing test cases. So, M will represent a suspicious set of statements.

Like for T , transaction variables are partitioned into two disjoint subsets: variables {T + t } representing positive test cases T + and variables {T - t } denoting negative test cases T -. We define also the frequency measure on T + (resp.

T -) as freq + (M ) = t∈T + T + t (resp. freq -(M ) = t∈T -T - t ).
To reduce the redundancy among the extracted patterns, we have to impose closed f req (M ), which states that M must be a closed pattern w.r.t the frequency measure. In our case, this constraint is imposed on T + and T - to ensure that M has no superset with the same frequencies in the two sets. We encode this constraint using equation ( 3), which is a decomposition of the closedness constraint (2) on T + and T -.

closed f req (M ) ≡ ∀i ∈ I : (M i = 1) ↔ t∈T + T + t × (1 -d + t,i ) = 0 ∧ t∈T - T - t × (1 -d - t,i ) = 0 (3)

top-k suspicious patterns

The intuition behind the most of fault localization methods is that statements that appear in the failing test cases are more likely to be suspicious, while the statements that appear only in the traces of passed executions are more likely to be guiltless (Eric [START_REF] Wong | A family of code coverage-based heuristics for effective fault localization[END_REF][START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF][START_REF] Renieres | Fault localization with nearest neighbor queries[END_REF]. To extract the most suspicious patterns (i.e. set of statements), we define a dominance relation R between patterns.

Definition 9 (Dominance relation) Given a bipartition of T into two disjoint subsets T + and T -, a pattern M dominates another pattern M (denoted M R M ), iff:

f req -(M ) > f req -(M ) ∨ (4) [(f req -(M ) = f req -(M )) ∧(f req + (M ) < f req + (M ))]
The dominance relation states that

M R M , if M is more frequent than M in T -. If M
= R M ), iff f req -(M ) = f req -(M ) ∧ f req + (M ) = f req + (M ) ∧ (M ⊂ M ) ∧ (M ⊂ M ) (5)
The indistinct relation states that two patterns can have the same level of suspiciousness. According to the dominance relation R , we define a top-k suspicious patterns.

Definition 11 (top-k suspicious pattern) A pattern M is a top-k suspi- cious pattern (according to R ) iff ∃P 1 , . . . P k , ∈ L I , ∀1 ≤ j ≤ k, P j R M .
Thus, M is a top-k suspicious pattern if there exists no more than (k -1) patterns that dominate M . A set of top-k suspicious patterns is defined as follows:

{M ∈ L I | elementary(M ) ∧ ∃P 1 , . . . P k , ∈ L I , ∀1 ≤ j ≤ k, P j R M }
The constraint elementary(M ) allows to specify that M must satisfy a basic property. In our case, we impose that searched suspicious patterns have to satisfy the following property:

elementary(M ) ≡ closed f req (M ) ∧ f req -(M ) ≥ 1 ∧ size(M ) ≥ 1,
which means that a suspicious pattern M must be closed in positive and negative test cases, must include at least one statement and must appear at least in a negative test case.

Mining top-k suspicious patterns

This section shows how the top-k suspicious patterns can be extracted using constraints dynamically posted during search [START_REF] Rojas | Mining (soft-) skypatterns using dynamic CSP[END_REF]. The main idea is to exploit a dominance relation (noted R ) between sets of statements to produce a continuous refinement on the extracted patterns thanks to constraints dynamically posted during the mining process. Each dynamic constraint will impose that none of the suspicious patterns already extracted is better (w.r.t R ) than the next pattern (which is searched). This process stops when no better solution can be obtained.

Algorithm 1 extracts the top-k patterns (i.e., most suspicious patterns) according to the dominance relation R . It takes as input the positive and negative test cases (T + and T -), a positive integer k, and returns as output top-k suspicious patterns. The algorithm starts with a constraint store equal to elementary(M ) (line 1). First, we search for the k first suspicious patterns that are solutions of the current constraint store, using the SolveN ext(C) function (lines 2-4). The SolveN ext(C) function asks a CP solver to return a solution of C which is different from the previous returned solutions. Initially, the first call of SolveN ext(C) returns the pattern P 1 . The second call will return a pattern P 2 = P 1 , and so on. During the search, a list of top-k suspicious patterns S is maintained. Once the k patterns are found, they are sorted according to decreasing order R (line 5). Thereafter, each time a new pattern is found, we remove from S the least preferred pattern w.r.t. R (line 10), we add the new pattern to S according to the decreasing order R (line 11) and we add dynamically a new constraint (M R S k ) at line 7 stating that the new searched pattern M must be better w.r.t. R than the least pattern in the current top-k list S. Thus, the next solutions should verify both the current set of constraints store and the new constraints added dynamically. This process is repeated until no pattern is generated.

Fault localization by processing top-k suspicious patterns

Our top-k algorithm returns an ordered list of k best patterns S = S 1 , . . . , S k . Each pattern S i represents a set of statements that can explain and locate the fault. From one to another, some statements are appearing/disappearing according to their frequencies in the positive/negative datasets (f req -, f req + ).

We propose in this section an algorithm (see Algorithm 2) that takes as input the top-k patterns and returns a ranked list Loc of most accurate suspicious statements enabling to better locate the fault (line 14). The returned list Loc includes three computed ordered lists noted Ω 1 , Ω 2 and Ω 3 . Elements of Ω 1 are ranked first, followed by those of Ω 2 , then by elements of Ω 3 , which contains the least suspicious statements.

Algorithm 2 starts by merging the indistinct patterns of S (def. 10), i.e. thus having the same frequencies in failing and passing test cases, or equivalently having the same level of suspiciousness. This treatment is achieved by the function M erge(S) (line 1) which returns a new list SM. A pattern resulting from merging mutliple patterns is a single pattern including all statements contained in original patterns. Let us note that the returned list SM may be equal to the initial top-k list S if there are no pair of patterns with the same frequencies (f req -, f req + ). Then, Algorithm 2 initializes Ω 1 and Ω 3 to the empty list. Ω 2 is initialized by statements that appear in the most suspicious pattern of SM (i.e. SM 1 ) (line 3).

From this set of most suspicious statements Ω 2 , Algorithm 2 will try to differentiate these statements by taking advantage of the three following properties 1, 2 and 3.

Property 1 Given a top-k patterns S, SM 1 is an over-approximation of the fault-location:

∀S i ∈ S : (f req - 1 = f req - i ) ∧ (f req + 1 = f req + i ) ⇒ S i ⊆ SM 1
Thus, SM 1 contains most likely the faulty statement, i.e. statements that are most frequent in the negative dataset and less frequent in the positive dataset. However, it may also contains other statements that are suspicious for a certain degree (over-approximation). That is why, in Algorithm 2, Ω 2 is initialized to SM 1 at line 3.

Property 2 Given the sets of patterns SM (resulting from top-k patterns S) and an over-approximation of the fault-location SM 1 (Ω 2 ), some statements of SM 1 will disappear in the next

SM i ∈ SM : (SM 1 \ SM i ) = ∅.
In Algorithm 2, statements that disappear from Ω 2 and appear in a given SM i are noted ∆ D (line 5). According to the frequency values, we have two cases:

1. SM i has the same frequency as SM 1 in the negative dataset but SM i is more frequent in the positive dataset than SM 1 . Thus, statements of ∆ D are less frequent in the positive dataset compared with those of (Ω 2 \ ∆ D ). So, statements of ∆ D are more suspicious than the remaining statements in Ω 2 \∆ D and should be ranked first (removed from Ω 2 (line 8) and added to Ω 1 (line 7)). 2. SM i is less frequent than SM 1 in the negative dataset. Again, statements of ∆ D should be ranked first and added to Ω 1 .

Consequently, Ω 1 contains the most suspicious statements derived from the initial state of Ω 2 . The remaining statements in Ω 2 are those appearing in all patterns SM i and are ranked second in terms of suspiciousness.

Property 3 Given the sets of patterns SM, and an over-approximation of the fault-location SM 1 , some statements will appear in the next SM i ∈ SM :

(SM i \ SM 1 ) = ∅.
According to the property 3, the statements that are not in SM 1 and that appear in a given SM i (line 9) are noted ∆ A . So, ∆ A should be added to Ω 3 (line 13) and ranked after Ω 2 as the least suspicious statements in the program (line 14).

In summary, the first ordered list Ω 1 ranks the statements of Ω 2 according to their order of disappearance in SM 2 ..SM l (line 4). At the end of the algorithm, Ω 2 will contain the remaining less suspicious statements as compared to the ones of Ω 1 . Finally, the ordered list Ω 3 will contain statements, which do not belong to Ω 2 but appears gradually in SM 2 to SM l .

Algorithm 2: Fault localization using top-k patterns Input top-k patterns S = S 1 , . . . , S k , frequecies of each S i :

(f req + i , f req - i ) Output an ordered list of suspicious statements Loc = Ω 1 , Ω 2 , Ω 3 1 SM ← M erge(S); 2 Ω 1 ← ; Ω 3 ← ; Loc ← ; 3 Ω 2 ← SM 1 ; l ← |SM|; 4 foreach i ∈ 2..l do 5 ∆ D ← Ω 2 \ SM i ; 6 if ∆ D = ∅ then 7 Ω 1 .add(∆ D ); 8 Ω 2 .removeAll(∆ D ); end 9 ∆ A ← SM i \ SM i-1 ; 10 ω ← ∅; 11 foreach b ∈ ∆ A do 12 if (∀ω ∈ Ω 3 , ∀ω" ∈ Ω 1 : b / ∈ ω ∧ b / ∈ ω") then ω ← ω ∪ {b} end 13 if ω = ∅ then Ω 3 .add(ω) end 14 Loc.addAll(Ω 1 ); Loc.add(Ω 2 ); Loc.addAll(Ω 3 );
return Loc;

Running Example

In this section, we give an illustrative example to show the result of our method through a simple program named Character counter, introduced in (Gonzalez-Sanchez et al, 2011) and given in Figure 1. The program contains ten executable statements, executed on eight test cases noted from tc 1 to tc 8 with failing (tc 1 to tc 6 ), and passing test cases (tc 7 and tc 8 ). In this example, the fault is introduced at line 3 where the correct statement "let += 1" is replaced by "let += 2". Figure 1 reports the coverage of each statement with the value 1 if the statement is executed at least once by the test case, 0 otherwise. According to our approach, the coverages of failing (resp. passing) test cases form transactions of the negative (resp. positive) dataset T -= {tc 1 , tc 2 , tc 3 , tc 4 , tc 5 , tc 6 } (resp. T + = {tc 7 , tc 8 }). Our constraint based model presented in Section 4 aims at extracting the k most suspicious patterns. We recall that the meaning given to the notion of suspiciousness is related to the frequency of a given statement in the negative and/or positive datasets. This means that the most suspicious statements are the ones with the highest negative frequency and/or the lowest positive frequency. In this example, we select k equals to the number of statements (i.e., k = 10) which is sufficient Table 5: The most suspicious statements with their ranks returned by Tarantula and Growth rate measure.

to reach a good accuracy (see Section 6) and then we give a ranking on all statements in P from the most suspect to the guiltless ones by processing the top-k patterns. Table 3 reports the top-k ranking of the different patterns computed by the first step of our approach with their respective frequencies. Table 4 gives the results obtained by Algorithm 2 taking as input Table 3 as a second step of our approach.

Tarantula and Growth rate measures: Tarantula's formula, which computes the suspiciousness of a statement e in the considered program, is defined as follows:

suspiciousness(e) = f ailed(e) totalf ailed passed(e) totalpassed + f ailed(e) totalf ailed (6) Formula ( 6) is very similar to the growth rate formula (given in definition 6) which measures the emergence of a given pattern from a dataset to another one (i.e., in our case, from negative test cases to positive test cases). The two formulas evaluate the amount of the negative frequency compared with the positive frequency; but Tarantula avoids dividing by a null positive frequency. In fact, the Tarantula values have the same increase as the growth rate, and consequently, they both give the same ranking of the statements suspiciousness, as illustrated in Table 5.

According to the Algorithm 2 (line 1), as a first step and before ranking the statements, we merge all patterns with the same frequencies (i.e., the same level of suspiciousness). From Table 3, we have ∀i < 9, SM i = S i , and the two last patterns S 9 , S 10 will be merged in one pattern SM 9 that will contain all statements in S 9 and in S 10 .

By comparing the results and the ranking strategy of our method (Table 4 ) with those given by Tarantula in Table 5 (2 nd column) and the growth rate (3 th column) we draw the following observations:

• First, in our results, the fault, that is introduced in e 3 , is ranked first and is the most suspicious statement. According to our Algorithm 2, this statement is in the most suspicious pattern SM 1 having frequencies (f req + = 0, f req -= 6). We point out that e 3 disappears in the next pattern SM 2 having frequencies (f req + = 1, f req -= 6).

Here, e 3 appears in SM 1 and disappears in the next SM i , which means that e 3 must be added as a disappeared statement to ∆ D (at line 5) and removed from Ω 2 (line 8). But, Tarantula considers the three statements e 3 , e 5 , e 7 as the most suspicious with the same rank although these statements have different frequencies in negative test cases. In fact, Tarantula is unable to distinguish between them.

• By using Tarantula, the statements {e 5 , e 7 } are ranked before {e 2 , e 4 , e 6 }, whereas statements {e 2 , e 4 , e 6 } are more suspicious since they are more present in failing test cases than {e 5 , e 7 }. This weak ranking is due to the low precision provided by Tarantula's formula (6). This weakness is offset by the notion of appearing statements. Let us come back to our example, the statements e 4 and e 6 are not in the most suspicious pattern SM 1 and they appear in the next patterns. Consequently, these statements are added as appearing statements to ∆ A (line 9) and then added to Ω 3 . • Concerning the statements {e 8 , e 9 } and {e 1 , e 10 }, Tarantula reveals that these statements have the same suspiciousness without distinguishing them.

Our approach gives the same results and the same suspiciousness because we got them in the same time (same group). In our approach the statement e 7 has also the same suspiciousness as {e 8 , e 9 } because these statements appear in different patterns but these patterns have the same frequencies, so these statements will be in the same pattern after the merging step.

This comparison between the effectiveness of Tarantula and our constrained mining method reveals clearly the accuracy of our method which exploits better the relationships between the statements of the program.

Experiments

This section reports several experimental studies. First, we present the benchmark programs we used for our experiments (see Sect. 6.1). Second, we detail our experimental protocol for evaluating our approach (see Sect. 6.2). Third, we study and discuss the influence of parameter k for mining top-k suspicious patterns (see Sect. 6.3). Fourth, we compare our results with those obtained with Tarantula [START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF] and we study the impact of varying the size of the test cases on performances of the two approaches (see Sects. 6.4 and 6.5). Finally, we study and analyze the CPU times of our approach (see Sect. 6.6).

Benchmark programs

As benchmark programs, we considered the Siemens suite, which is the most common program suite used to compare the different fault localization techniques. The suite contains seven C programs, faulty versions of these programs, and test suites for each category of program2 .

We recall here that the Siemens programs suite is assembled for fault localization studies as well as for fault detection capabilities of control-flow and data-flow coverage criteria. In our experiments, we exclude 21 versions that are, roughly speaking, out of scope of localization task (e.g., segmentation faults3 ). On the whole, we have 111 programs with faults summarized in Table 6.

Experimental protocol

First of all, we need to know which statement is covered by a given execution. For that, we use the Gcov4 profiler tool to find out the statements that are actually executed by a given test case. Gcov can be used also to know how often each statement is executed. This being said, in our approach we need only the coverage matrix. Thus, we run the program P with a test case and we compare the returned result to the expected one. If the returned result matches with the expected one, we add the according test case coverage to the positive transactional dataset; otherwise we add it to the negative one.

We implemented our approach as a tool named F-CPminer. The implementation was carried out in Gecode5 , an open and efficient constraint programming solver. Our implementation includes Algorithm 1 for mining top-k suspicious patterns (see Sect. 4) and Algorithm 2 that processes the top-k patterns and returns at the end an accuracy fault localization with a ranking on statements according to their suspiciousness. All experiments were conducted on a 3.10 GHz processor Intel Core i5-2400 with 8 GB of memory, running Ubuntu 12.04 LTS.

To make a fair comparison between our F-CPminer and Tarantula, we have implemented Tarantula and we evaluated the statement suspiciousness as presented in [START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF]. We used a common metric in fault localization, the exam score [START_REF] Wong | A survey of software fault localization[END_REF] that measures the effectiveness of a given fault localization technique. The exam score gives the percentage of statements that a developer will check before the one containing the fault. It is clear that the best approach is the one that has the lowest percentage of exam score.

Tarantula and F-CPminer can return a set of equivalent statements in terms of suspiciousness (i.e., with the same suspiciousness degree). In this case, the effectiveness depends on which statement is to check first. For that reason, we report two exam scores, the optimistic and the pessimistic one, denoted respectively o-exam and p-exam. We talk about o-exam (resp. p-exam) when the first (resp. last) statement to be checked in the set of equivalent statements is the faulty one. We also define a third metric, ∆-exam = oexamp-exam, representing the margin of the exam score. In other words, ∆-exam represents the distance between the optimistic and the pessimistic score.

k = n/5 k = n/2 k = n k = n × 2 k = n ×

Comparing various values of k for mining top-k suspicious patterns

Our first experiment aims to select the best value of k leading to the best ∆-exam score. Here, a small value for ∆-exam means that the method gives small sets of statements with equivalent suspiciousness, which leads to an accurate localization. We recall that Algorithm 1 for mining top-k suspicious patterns takes as input parameter k. For this experiment, we selected three programs from Replace, Schedule2 and Tcas of different sizes (resp. 245, 128 and 65 executable statements). We have varied the value of k from (n/5) to (5×n), where n represents the program size in terms of executable statements.

Figure 2 shows the impact of k on the precision of the returned ∆-exam score (i.e., the distance between the optimistic and pessimistic exam values). The first observation that we can make is that by increasing the value of k, we reduce the distance between o-exam and p-exam (i.e., ∆-exam). The second observation is that starting from a value of k greater or equal to the program size measured in executable statements (i.e., n), ∆-exam becomes stable. Throughout the rest of this section and according to our tests, k is fixed to the program size.

(0,0) (0,0) (4,4) (1,1) (2,0) (6,4) (29,24) tie game (1,2) (3,3) (2,2) (0,0) (0,0) (3,5) (12,13) (21,25) (p-exam, o-exam)
6.4 F-CPminer versus Tarantula Table 7 gives an exam score based comparison between F-CPminer and Tarantula on the 111 faulty programs of Siemens suite. For each class of program (e.g., Tcas includes 37 faulty versions), we report the averaged number of positive test cases |T + |, the averaged number of negative test cases |T -|, the value of k used in F-CPminer for mining top-k suspicious patterns (corresponding to the program size), the averaged p-exam, o-exam and ∆-exam values ± the standard deviation for each class of programs.

The first observation is that F-CPminer wins on 5 program classes out of 7 according to the reported o-exam and p-exam values. For instance, if we take PrintTokens2 class with its 9 faulty versions (see Table 6), using F-CPminer the fault is localized after examining only 2.66% of the code in the pessimistic case (p-exam) with a standard deviation less than 4% and only 1.72% in the optimistic case (o-exam) with a standard deviation less than 3%. On the other side, Tarantula needs to examine 16.44% in the pessimistic case and 15.5% of the code in the optimistic case, in both cases, the standard deviation is greater than 17%.

Our second observation relates to ∆-exam values that are better in Tarantula than F-CPminer in most cases. It is worth mentioning that ∆-exam relies heavily on the value of k (see Figure 2). Using k equal to the program size leads us to a difference not exceeding 4% for one program class and 1% for the other classes comparing to Tarantula.

In order to complement the results given by Table 7, we report in Table 8 a pairwise comparison between Tarantula and F-CPminer on the 111 faulty programs. For each category of program (e.g., Tcas), we give the number of programs where Tarantula or F-CPminer wins by comparing the exam scores (o-exam and p-exam). We also report the case when the exam values are equal (i.e., tie game).

The results of Table 8 match with the results given in Table 7, where F-CPminer is the winner on 5 out of the 7 classes. For instance, if we take the Tcas program class and by comparing the o-exam (resp. p-exam) value, F-CPminer is more accurate on 19 (resp. 20) faulty programs, while Tarantula is better only on 6 (resp. 4) faulty programs. Finally, both approaches obtain the same o-exam (resp. p-exam) on 12 (resp. 13) programs.

From Tables 7 and8, a general observation is that F-CPminer is highly competitive as compared to Tarantula. Indeed, F-CPminer wins on almost 60% of benchmark programs (61 to 62 out of 111 programs), while Tarantula is doing better on 25% of benchmark programs ( 24 to 29 out of 111 programs). Finally, the two approaches behave similarly on 15% of the benchmark programs ( 21 to 25 out of 111 programs).

Figure 3 shows the effectiveness comparison, based on p-exam and oexam, between Tarantula and F-CPminer. The x-axis reports the exam values while y-axis reports the cumulative percentages of located faults on the 111 programs. Let us start with the pessimistic case (i.e. p-exam). Until 10% of exam, both approaches behave similarly by locating more than 40% of faults. From 10% to 20%, we observe a difference of 10% of located faults in favour of F-CPminer (i.e., 60% for F-CPminer instead of 50% of located faults for Tarantula). For values of exam in the interval [20%, 40%], this difference is reduced up to 6%. Between 60% and 70%, Tarantula arrives to catch up F-CPminer by locating more than 84% of the faults. However, F-CPminer enables to detect 100% of faults with a exam value equal to 80%, while Tarantula needs to reach 100% of exam to locate all the faults. For the optimistic case (i.e. o-exam) showed with dashed curves, F-CPminer is acting quickly from the beginning by locating more faults than Tarantula. It is important to stress that the two curves do not intersect and the one of F-CPminer is always above the one of Tarantula. Let us note that after 30% of exam, F-CPminer in the pessimistic case detects the same percentage of faults (i.e., 74%) than Tarantula in the optimistic case.

In order to strengthen the previous results, we carried out a statistical test using the Wilcoxon Signed-Rank Test, where we are not able to assume a normal distribution of the population (Lyman [START_REF] Ott | An introduction to statistical methods and data analysis[END_REF]. Here, the one-tailed alternative hypothesis is used with the following null hypothesis (Tarantula more efficient than F-CPminer): H 0 : The exam score using F-CPminer ≥ The exam score using Tarantula. That is, the alternative hypothesis H 1 states that our F-CPminer examined a fewer number of statements before locating the fault than Tarantula (F-CPminer more efficient than Tarantula). With this test, we are able to conclude that the use of F-CPminer is more efficient than Tarantula (H 1 accepted) with 79.99% confidence in the pessimistic case (p-exam) and 97.99% confidence in the optimistic case (o-exam). This test strengthens the previous results on the comparison of F-CPminer with Tarantula.

Overall, we can conclude that F-CPminer enables to locate most of the faults more quickly than Tarantula in terms of effectiveness (i.e., exam metric). Moreover, by considering the cumulative exam over all the 111 faulty programs (i.e. the total number of executable statements that have to be checked for the whole of the 111 faulty programs), F-CPminer has to check 2861 of statements against 3482 for Tarantula using p-exam (gain of 17%). With o-exam, F-CPminer needs to check 1808 statements against 2496 for Tarantula (gain of 27%).

Finally, we point out the following observations over the tested 111 programs :

-For 110 programs, the faulty statement is in the first pattern (i.e., SM 1 ).

-For 96 programs, the faulty statement is in the first statement disappearing from SM 1 . -For 14 programs, the faulty statement is located in Ω 2 .

-For only one program, the faulty statement is located in Ω 3 .

These observations show clearly the effectiveness of the ranking strategy adopted by Algorithm 2.

Impact of the size of test cases on fault localization accuracy

In this section, we study the impact of test cases on F-CPminer and Tarantula. For that, we have varied the number of test cases given as inputs by varying the sizes of datasets T + and T -. For each program, we reduced the size of its datasets from 100% to 10% by removing randomly at each time 10% of test cases. At each time, we report the exam score (o-exam and p-exam) for F-CPminer and Tarantula. In this section, we selected four programs to present (PrintTokens2-v3, Tcas-v28, TotInfo-v18 and Replace-v22) that we consider quite representative in terms of program and datasets sizes. Figure 4 reports the results on the selected programs.

PrintToken2-v3. For this instance of program, when considering the full datasets (100% of test cases), F-CPminer and Tarantula obtain approximately the same results in terms of exam score. Once we start reducing the number of test cases, we observe that F-CPminer keeps approximately the same accuracy until 10% of test cases. In the other side, Tarantula exhibits a chaotic behavior where the accuracy decreases (i.e (exam score) increases) significantly after a reduction of 60% of test cases by reaching a score of 30% of code to exam. These observations are made on 6 versions where only one is from the same PrintToken2 class.

Tcas-v28. For this instance, F-CPminer and Tarantula obtain approximately the same results in terms of exam score with a slight gain in favour of our tool. We can observe that datasets reduction have no substancial impact on the accuracy of the two tools. This result is observed on 32 versions where 14 versions belong to the same Tcas class.

TotInfo-v18. Here, F-CPminer shows a stable behavior and more or less the same accuracy during the datasets reduction, while Tarantula is greatly impacted by this reduction (i.e., the exam score changes in a strange manner). The same observations are made on 22 versions where 8 of them are from TotInfo class.

Replace-v22. For this instance, we can see that the two approaches are more or less affected by reducing the size of the datasets, Tarantula is still more stable until a reduction of 50% of test cases. F-CPminer is stable from 75% to 35% of datasets reduction. For this case, the same observations (i.e., the two approaches are affected by the reduction) are made on 34 other versions from which 7 programs belong to the Replace class.

These observations are particularly informative and highlight the fact that our approach can be less sensitive to the number of considered test cases. This is especially true when we have 60 programs out of 111 where F-CPminer behaves in the same manner (stable behavior) as the three first programs (see fig. 4). And this is especially true where on 62 programs Tarantula follows a chaotic behavior and only on 17 programs, the behavior of Tarantula is more stable than F-CPminer. Thanks to the quality of our top-k patterns extraction and to the processing step, which enables to analyze the dependencies between the extracted patterns. This is not the case for Tarantula, where the adding or removing of test cases can lead to less accurate results.

Analyzing CPU times for F-CPminer

This section analyzes the CPU times of our approach. It is important to recall that Tarantula approach evaluates the suspiciousness degree of each statement using the formula 6 without taking in account any statement dependencies. Therefore, the combinatorial explosion due to the possible combinations of statements is not tackled by Tarantula. Consequently, the CPU times obtained by Tarantula are very shorts (in milliseconds).

Table 9 reports for each class of program both the average CPU times and the standard deviation for the two steps performed by F-CPminer (i.e. the top-k extraction corresponding to Algorithm 1 and the post-processing step corresponding to Algorithm 2).

The first observation that can be made is that extracting top-k patterns is the most costly step for F-CPminer. This is in part explained by the very high number of candidate patterns (i.e. |L I |). However, in our experiments, the CPU times are not exceeding 235 seconds for the worst case (see Replace programs) and it is less than 0.2 seconds in the best case (see Tcas programs). The second observation that can be made is that the CPU times spent by the post-process step for a final localization of the faulty statement is negligible (of order of milliseconds for the best and worst cases).

Figure 5 shows the CPU times variation for extracting top-k patterns for the 29 faulty versions (from v1 to v29) of the Replace program. For this class, the standard deviation of CPU time is quite large (i.e., 86.02s) according to an average CPU time equal to 147.69s. The deviation is not negligible due to the fact that all program versions take the same test cases but different T + and T -according to the introduced fault. For instance, for v5, is especially true when the extracted patterns are the more frequent ones in the negative dataset and less frequent ones in the positive dataset.

Impact of the size of test cases and parameter k on CPU times. Figure 6 shows the impact of varying the size of the dataset (i.e. test case) and the value of k on the CPU times for the top-k extraction. For this experiment, we have selected the program Replace(v5), which represents the worst case for extracting top-k. We can observe that CPU time decreases considerably when reducing the number of test cases (see Figure 6(a)). The same trend is observed when varying the value of k (see Figure 6(b)).

In fact, the extraction of top-k patterns is closely related to the size of the transactional datasets (test cases coverage) and the value of k. The behavior of Figure 6(b) can be explained by the fact that with a large value of k, a large number of patterns are extracted and the least patterns of the top-k list have low values of frequency. This increases drastically the number of candidate patterns to be explored by Algorithm 1 and thus the CPU times. In contrast, with small values of k, the least patterns in the top-k list tend to have high values of frequency. Consequently, constraints added dynamically during the mining process (line 7 of Algorithm 1) will refine the pruning condition leading to more and more powerful pruning of the search space. This explains the decreasing of CPU times.

Related works

In this section, we present some related works in the area of fault localization.

Using failing/passing executions. As introduced previously, Tarantula [START_REF] Jones | Empirical evaluation of the tarantula automatic fault-localization technique[END_REF] is one of the most popular fault localization technique that records information linked to failing/passing executions in terms of how each test case covers statements. Another technique rather close to Tarantula is the one proposed by Cleve and Zeller [START_REF] Cleve | Locating causes of program failures[END_REF] and based on program states. The technique compares states of passing and failing test cases. In a previous work [START_REF] Zeller | Isolating cause-effect chains from computer programs[END_REF], Zeller shows that locating a fault just by considering the search space (variables, values) is not sufficient in general. Indeed, it is possible to return statements as faulty ones by comparing the different states in passing/failing test cases, but a fault in a program can produce a difference on all next states in the program. In [START_REF] Cleve | Locating causes of program failures[END_REF], a search during execution is performed to locate the first transition that leads to a fail. But these techniques do not differentiate between two fail-ing (resp. passing) test cases, and consequently they ignore the dependencies between statements that can help us to locate the fault.

Using dependence analysis. Other techniques take in account the cause effect chains with a dependence analysis by using, for instance, program slicing [START_REF] Agrawal | Debugging with dynamic slicing and backtracking[END_REF]. The disadvantage here is that fault can be located in a quite large slice (static slicing) and/or can be time/space consuming (dynamic slicing). In the same perspective, Renieris and Reiss [START_REF] Renieres | Fault localization with nearest neighbor queries[END_REF] use the notion of nearest neighbor, where they confront a failing trace with the nearest passing trace. Here the distance between two traces is expressed with the difference between the set of executed statements. In the case where no nearest passing trace can be obtained, the technique builds the program dependence graph and checks the adjacent nodes of the failing trace one by one with the hope of finding the location of the fault.

Using data mining. In [START_REF] Cellier | Dellis: A data mining process for fault localization[END_REF], Cellier et al. propose a data mining process DeLLIS which computes program element clusters and shows dependencies between program elements. They compute all differences between execution traces and, at the same time, gives a partial ordering of those differences. In [START_REF] Cellier | Formal concept analysis enhances fault localization in software[END_REF], Cellier et al. propose a methodology that combines between association rules to search for possible causes of failure and formal concept analysis to assist in fault localization. They try to identify rules between statement execution and corresponding test case failure, and then measure the frequency of each rule. Nessa et al. [START_REF] Nessa | Software fault localization using n-gram analysis[END_REF] generate statement subsequences of length N , referred to as N-grams, from the trace data. The failed execution traces are then examined to find the N-grams with a rate of occurrence higher than a certain threshold in the failed executions. A statistical analysis is conducted to determine the conditional probability that an execution fails given that a certain N-gram appears in its trace.

Our approach aims at exploiting recent progress on cross-fertilization between data mining and constraint programming in order to model the fault localization problem as finding the top-k patterns of statements occurring more frequently in failing executions and less frequently in the passing executions. These top-k patterns are then processed by analysing their dependencies, in order to infer a ranking on the suspiciousness degree of the statements.

Conclusion

In this paper we have proposed a new approach based on itemset mining and constraint programming to deal with fault localization problem. Our approach proceeds in two steps. In the first step, we have formally defined the problem of locating faults in programs as a mining task using CP for modelling and solving the arising constraints. Solving the underlying CP model enables us to get the top-k most suspicious set of statements . The second step aims at ranking in a more accurate way the whole top-k statements by taking benefit of two main observations i) where faults are introduced in a program can be seen as a pattern (set of statements), which is more frequent in failing executions than passing ones; ii) the difference between a more suspicious pattern and a less suspicious one is a set of statements that appears/disappears in one or the other; this difference helps us to know more about the location of the fault. We have shown how these two properties can be exploited in an ad-hoc ranking algorithm producing accurate localization. Finally we have compared experimentally our approach implemented in F-CPminer with Tarantula on a set of faulty programs. The results we obtained show that our approach enables to propose a more precise localization as compared to Tarantula.

As future works, we plan to experiment our approach on programs with complex faults (more than one faulty statement). We also plan to explore other observations on the behavior of a faulty program and adding them as constraints for mining the location of faults.

  Definition 8 (top-k patterns) Let m be a measure, and k an integer. top-k is the set of k best patterns according to m : {x ∈ L I | freq T (x) ≥ 1 ∧ ∃y 1 , . . . , y k ∈ L I : ∀1 ≤ j ≤ k, m(y j ) > m(x)} Example 5 In our running example, the top-4 closed frequent patterns (i.e., m = f req) are: B 6 , E 6 , BE 5 , BCD 5 .

  and M have the same frequency in T -, M should have a less positive frequency to dominate M . We also define the indistinct relation = R between patterns. Definition 10 (Indistinct relation) Two patterns M and M are indistinct (denoted by M

Algorithm 1 :

 1 Extraction of top-k most suspicious patterns S 1 , ..., S k Input T + , T -, k Output S: top-k most suspicious patterns 1 C ← elementary(M ); S ← ∅; i ← 1; 2 repeat 3 P ← SolveN ext(C); 4 if P = ∅ then S.add(P ); i ← i + 1 until (i > k) or (P = ∅); 5 Sort S according to decreasing order R ; 6 while P = ∅ do 7 C.add(M R S k ); 8 P ← SolveN ext(C); 9 if P = ∅ then 10 S.remove(S k ); 11 Insert P in S according to decreasing order R ; end end return S;

Fig. 2 :

 2 Fig. 2: Comparing various values of k.

Fig. 3 :

 3 Fig. 3: Tarantula and F-CPminer effectiveness comparison.

Fig. 4 :

 4 Fig. 4: Impact of the size of the datasets on the performances of F-CPminer and Tarantula.

  (|T -|, |T + |) = (271, 5271) and F-CPminer requires about 400 s. to extract the top-k patterns. For comparison, for v24, (|T -|, |T + |) = (3, 5539) and F-CPminer needs less CPU times to complete the extraction of top-k pattens (22 s). In fact, extracting top-k patterns from a small negative dataset (i.e.T -) and a considerable positive dataset (i.e. T + ) is not time consuming. This

Fig. 5 :

 5 Fig. 5: CPU times variation for top-k extraction (Algorithm 1) for Replace program.

Fig. 6 :

 6 Fig. 6: Impact of varying the size of the dataset and the value of k on CPU times.

  Definition 1 (Test case) A test case tc i is a tuple D i , O i , where D i is a collection of input settings for determining whether a program P works as expected or not, and O i is the expected output.

Definition 2 (Passing and failing test case) Let D i , O i a test case and A i be the current output returned by a program P after the execution of its input D i . If A i = O i , tc i is considered as a passing (i.e. positive), failing (i.e. negative) otherwise.

Definition 3 (Test suite) A test suite T = {tc 1 , tc 2 , ..., tc m } is a set of m test cases that are intended to test whether the program P follows the specified set of requirements.

  : {e 1 , e 2 , e 3 , e 4 , e 10 } 0 5 4 S 5 : {e 1 , e 2 , e 4 , e 10 } 1 5 5 S 6 : {e 1 , e 2 , e 3 , e 4 , e 6 , e 10 } 0 4 6 S 7 : {e 1 , e 2 , e 4 , e 6 , e 10 } 1 4 7 S 8 : {e 1 , e 2 , e 3 , e 4 , e 5 , e 10 } 0 3 8 S 9 : {e 1 , e 2 , e 3 , e 4 , e 6 , e 7 , e 10 } 0 2 9 S 10 : {e 1 , e 2 , e 3 , e 4 , e 6 , e 8 , e 9 , e 10 }

	Patterns of statements	f req +	f req -	Rank
	S 1 : {e 1 , e 2 , e 3 , e 10 }	0	6	1
	S 2 : {e 1 , e 2 , e 10 }	1	6	2
	S 3 : {e 1 , e 10 }	2	6	3
	S 4 0	2	9

Table 3 :

 3 The top-k suspicious patterns returned by Algorithm 1.

	Statements	Rank	List
	e 3	1	Ω 1
	e 2	2	Ω 1
	{e 1 , e 10 }	4	Ω 2
	e 4	5	Ω 3
	e 6	6	Ω 3
	e 5	7	Ω 3
	e 7	10	Ω 3
	{e 8 , e 9 }	10	Ω 3

Table 4 :

 4 The most suspicious statements with their ranks returned by Algorithm 2.

	Statements	Tarantula	Growth rate	Rank
	e 3	1	∞	1
	e 5	1	∞	1
	e 7	1	∞	1
	e 2	0.66	2	4
	e 4	0.62	1.67	5
	e 6	0.57	1.33	6
	e 1	0.5	1	8
	e 10	0.5	1	8
	e 8	0.4	0.67	10
	e 9	0.4	0.67	10

Table 6 :

 6 Faulty programs from the Siemens suite (111 programs).

	Program	Description	Faulty versions	LOC	LEC Test cases
	Replace	Pattern replacement	29	514	245	5542
	PrintTokens2	Lexical analyzer	9	358	200	4056
	PrintTokens	Lexical analyzer	4	348	195	4071
	Schedule	Priority scheduler	5	294	152	2650
	Schedule2	Priority scheduler	8	265	128	2680
	TotInfo	Information measure	19	272	123	1052
	Tcas	Altitude separation	37	135	65	1578
	LOC: lines of code in the correct version -LEC: lines of executable code

Table 7 :

 7 F-CPminer versus Tarantula (exam score comparison).

					p-exam	o-exam	∆-exam
	Program	|T + | |T -| k	F-CPminer	Tarantula	F-CPminer	Tarantula	F-CPminer	Tarantula
	Replace	5450	92	245	12.34 ± 10.75	11.04 ± 13.95	10.47 ± 9.37	9.31 ± 12.12	1.87 ± 3.11	1.73 ± 3.14
	PrintTk2	3827	229	200	2.66 ± 3.26	16.44 ± 17.86	1.72 ± 2.18	15.50 ± 17.37	0.94 ± 1,48	0.94 ± 1.48
	PrintTk	4016	55	195	5.89 ± 3.70	13.58 ± 11.21	3.97 ± 3.46	11.66 ± 11.40	1.92 ± 1.71	1.92 ± 1.71
	Sched	2506	144	152	25.26 ± 26.51	6.71 ± 3.84	22.63 ± 25.39	5 ± 3.54	2.63 ± 1.44	1.71 ± 0.67
	Sched2	2646	34	128	44.82 ± 21.23	62.01 ± 28.73	30.46 ± 22.50	51.17 ± 26.08	14.36 ± 7.18	10.84 ± 5.53
	TotInfo	1015	37	123	11.51 ± 7.85	23.62 ± 16.30	5.56 ± 7.06	17.88 ± 14.62	5.95 ± 5.05	5.74 ± 5.14
	Tcas	1542	36	65	40.33 ± 27.50	43.95 ± 31.91	16.13 ± 10.84	20.12 ± 14.61	24.19 ± 17.30	23.83 ± 17.62

Table 8 :

 8 Pairwise comparisons between Tarantula and F-CPminer.

	Winner	Replace	Print2	Print	Sched Sched2	Totinfo	Tcas	Total
	F-CPminer	(12,12)	(6,6)	(2,2)	(1,1)	(7,7)	(14,14) (19,20) (61,62)
	Tarantula	(16,15)

Table 9 :

 9 Comparing the CPU times for the two steps of F-CPminer on different programs (in seconds).

	Program	k	Top-k extraction (Algo.1) Top-k processing (Algo.2)
	Replace	245	147.69 ± 86.02	0.051 ± 0.015
	PrintTokens2 200	146.25 ± 65.89	0.045 ± 0.013
	PrintTokens	195	61.49 ± 34.35	0.030 ± 0.003
	Schedule	152	27.41 ± 15.28	0.012 ± 0.004
	Schedule2	128	12.66 ± 05.37	0.016 ± 0.004
	TotInfo	123	2.53 ± 01.01	0.014 ± 0.006
	Tcas	65	0.16 ± 00.02	0.001 ± 0.000
		400
		300
	time(s)	200
		100
		0
			Versions

A complete description of Siemens suite can be found in[START_REF] Hutchins | Experiments of the effectiveness of dataflow-and controlflow-based test adequacy criteria[END_REF]

In our approach, it is possible to tackle versions with segmentation faults by considering them as negative examples. In these cases, the last instruction which is mostly at the origin of the error, is not seen by the execution trace. Moreover, Gcov is not able to generate properly the covering instructions. Thus, these faulty versions pose many out-of-scope issues which disturb a safe discussion of the experimental results.

https://gcc.gnu.org/onlinedocs/gcc/Gcov.html

www.gecode.org

Conflict of Interest: The authors declare that they have no conflict of interest.


Fanny Dufossé
email: fanny.dufosse@inria.fr

Kamer Kaya
email: kaya@sabanciuniv.edu

Ioannis Panagiotas
email: ioannis.panagiotas@ens-lyon.fr

Bora Uçar
email: bora.ucar@ens-lyon.fr

Further notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices

Keywords: doubly stochastic matrix, Birkhoff-von Neumann decomposition 2010 MSC: 15B51, 05C50, 05C70

The well-known Birkhoff-von Neumann (BvN) decomposition expresses a doubly stochastic matrix as a convex combination of a number of permutation matrices. For a given doubly stochastic matrix, there are many BvN decompositions, and finding the one with the minimum number of permutation matrices is NP-hard. There are heuristics to obtain BvN decompositions for a given doubly stochastic matrix. A family of heuristics is based on the original proof of Birkhoff and proceeds step by step by subtracting a scalar multiple of a permutation matrix at each step from the current matrix, starting from the given matrix. At every step, the subtracted matrix contains nonzeros at the positions of some nonzero entries of the current matrix and annihilates at least one entry, while keeping the current matrix nonnegative. Our first result, which supports a claim of Brualdi [Canad. Math. Bull. 25 (1982), pp. 191-199], shows that this family of heuristics can miss optimal decompositions. We also investigate the performance of two heuristics from this family theoretically.

Introduction

A square matrix A = [a ij ] ∈ R n×n is doubly stochastic if a ij ≥ 0 for all i, j and Ae = A T e = e, where e is the column vector of all ones. In other words, the entries of a doubly stochastic matrix are nonnegative, and the sum of entries in any row or column is equal to one. By Birkhoff's Theorem, there exist α 1 , α 2 , . . . , α k ∈ (0, 1] with k i=1 α i = 1 and k different permutation matrices P 1 , P 2 , . . . , P k such that

A = α 1 P 1 + α 2 P 2 + • • • + α k P k . ( 1
)
This representation is also called Birkhoff-von Neumann (BvN) decomposition. For a given matrix there are many different BvN decompositions. Marcus-Ree Theorem [START_REF] Marcus | Diagonals of doubly stochastic matrices[END_REF] states that for a dense matrix, there are decompositions with k ≤ n 2 -2n+2 permutation matrices. For a sparse, fully indecomposable matrix with τ nonzeros, the same result holds [START_REF] Brualdi | Convex polyhedra of doubly stochastic matrices: I. Applications of the permanent function[END_REF][START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] with k ≤ τ -2n+2. Brualdi [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF] gives lower bounds on the number of permutation matrices in any BvN decomposition of a given matrix. A recent work [START_REF] Dufossé | Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices[END_REF] shows that the problem of finding a BvN decomposition with the smallest number k of permutation matrices is strongly NP-complete.

There are heuristics to compute a BvN decomposition for a given matrix A. In particular, the following family of heuristics is based on the constructive proof of Birkhoff. Let A (0) = A. At every step j ≥ 1, find a permutation matrix P j having its ones at the positions of the nonzero elements of A (j-1) , use the minimum nonzero element of A (j-1) at the positions identified by P j as α j , set A (j) = A (j-1) -α j P j , and repeat the computations in the next step j + 1 until A (j) becomes void. Any heuristic of this type is called generalized Birkhoff heuristic. The original Birkhoff heuristic chooses a permutation matrix P j at step j which contains a one at the position of the minimum nonzero element of A (j-1) . Dufossé and Uçar [START_REF] Dufossé | Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices[END_REF] propose choosing a permutation matrix P j where the minimum nonzero element of A (j-1) identified by P j is maximum. We refer to this last heuristic as the Greedy heuristic.

Doubly stochastic matrices and their associated BvN decompositions have been used in several operations research problems and applications. Classical examples are concerned with allocating communication resources, where an input traffic is routed to an output traffic in stages [START_REF] Chang | On service guarantees for inputbuffered crossbar switches: A capacity decomposition approach by Birkhoff and von Neumann[END_REF]. Each routing stage is a (sub-)permutation matrix and is used for handling a disjoint set of communications. The number of stages correspond to the number of (sub-)permutation matrices. A recent variation of this problem appears in routing in data centers [START_REF] Kulkarni | Proc. 19th International Conference Integer Programming and Combinatorial Optimization (IPCO 2017)[END_REF]. BvN decompositions are also used to build preconditioners for solving sparse linear systems [START_REF] Benzi | Preconditioning techniques based on the Birkhoff-von Neumann decomposition[END_REF]. Here, the number k of the permutation matrices is related to the cost of applying the preconditioner.

Our contribution in this note is threefold. First, we answer an open question asked by Brualdi [3,. Brualdi predicts that there are BvN decompositions which cannot be written as a convex combination of other decompositions and which cannot be obtained by a generalized Birkhoff heuristic. We give an example affirming the prediction of Brualdi. Second, we investigate the worst case performance guarantee of the original Birkhoff heuristic and show that it can be very bad. To do so, we describe a family of matrices where the smallest number of permutation matrices required is 3 while the Birkhoff heuristic finds decompositions with n permutation matrices, for a matrix of size n × n. Kulkarni et al. [START_REF] Kulkarni | Proc. 19th International Conference Integer Programming and Combinatorial Optimization (IPCO 2017)[END_REF]Theorem 7] also present examples to show that the worst case performance of Birkhoff heuristic can be very bad. Our construction is more explicit and seems simpler. Our third contribution is to show an approximation guarantee of the Greedy heuristic.

A result on the polytope of BvN decompositions

Let S(A) be the polytope of all BvN decompositions for a given, doubly stochastic matrix A. The extreme points of S(A) are the ones that cannot be represented as a convex combination of the other decompositions. Brualdi [3, pp. 197-198] observes that any generalized Birkhoff heuristic obtains an extreme point of S(A), and predicts that there are extreme points of S(A) which cannot be obtained by a generalized Birkhoff heuristic. In this section, we substantiate this claim by showing an example.

Any BvN decomposition of a given matrix A with the smallest number of permutation matrices is an extreme point of S(A); otherwise the other BvN decompositions expressing the said point would have smaller number of permutation matrices.

Lemma 1. There are doubly stochastic matrices whose polytopes of BvN decompositions contain extreme points that cannot be obtained by a generalized Birkhoff heuristic.

We will prove the lemma by giving an example. We use computational tools based on a mixed integer linear programming (MILP) formulation of the problem of finding a BvN decomposition with the smallest number of permutation matrices. We first describe the MILP formulation.

Let A be a given n × n doubly stochastic matrix, and Ω n be the set of all n × n permutation matrices. There are n! matrices in Ω n ; for brevity, let us refer to these permutations by P 1 , . . . , P n! . We associate an incidence matrix M of size n 2 × n! with Ω n . We fix an ordering of the entries of A so that each row of M corresponds to a unique entry in A. Each column of M corresponds to a unique permutation matrix in Ω n . We set m ij = 1 if the ith entry of A appears in the permutation matrix P j , and set m ij = 0 otherwise. Let a be the n 2 -vector containing the values of the entries of A in the fixed order. Let x be a vector of n! elements, x j corresponding to the permutation matrix P j . With these definitions, MILP formulation for finding a BvN decomposition of A with the smallest number of permutation matrices can be written as:

minimize n! j=1 s j (2) subject to M x = a , (3)
1 ≥ x j ≥ 0, for j = 1, . . . , n! , (4)
x j = 1, for j = 1, . . . , n! , (5)
s j ≥ x j , for j = 1, . . . , n! , (6)
s j ∈ {0, 1}, for j = 1, . . . , n! . (7)
In this MILP, s j is a binary variable which is 1 only if x j > 0, otherwise 0. The equality (3), the inequalities (4), and the equality ( 5) guarantee that we have a BvN decomposition of A. This MILP can be used only for small problems. In this MILP, we can exclude any permutation matrix P j from a decomposition by setting s j = 0.

Proof of Lemma 1. Let the following 10 letters correspond to the numbers underneath a b c d e f g h i j 1 2 4 8 16 32 64 128 256 512 .

Consider the following matrix whose row sums and column sums are 1023, hence can be considered as doubly stochastic

A =       a + b d + i c + h e + j f + g e + g a + c b + i d + f h + j f + j e + h d + g b + c a + i d + h b + f a + j g + i c + e c + i g + j e + f a + h b + d       . ( 8
)
Observe that the entries containing the term 2 form a permutation matrix, for = 0, . . . , 9. Therefore, this matrix has a BvN decomposition with 10 permutation matrices. We created the MILP above and found 10 as the smallest number of permutation matrices by calling the CPLEX solver [START_REF]IBM ILOG CPLEX Optimizer[END_REF] via the NEOS Server [START_REF] Dolan | The NEOS Server 4.0 administrative guide[END_REF][START_REF] Czyzyk | The NEOS Server[END_REF][START_REF] Gropp | Optimization environments and the NEOS Server[END_REF]. Hence the described decomposition is an extreme point of S(A). None of the said permutation matrices annihilate any entry of the matrix. Therefore, at the first step no entry of A gets reduced to zero, regardless of the order of permutation matrices. Thus, this decomposition cannot be obtained by a generalized Birkhoff heuristic.

One can create a family of matrices with arbitrary sizes by embedding the same matrix into a larger one of the form Let us investigate the matrix A and its BvN decomposition given in the proof of Lemma 1. Let P a , P b , . . . , P j be the 10 permutation matrices of the decomposition, corresponding to a, b, . . . , j. We solve 10 MILPs in which we set s t = 0 for one t ∈ {a, b, . . . , j}. This way, we try to find a BvN decomposition of A without P a , without P b and so on, always with the smallest number of permutation matrices. The smallest number of permutation matrices in these 10 MILPs were 11. This certifies that the only BvN decomposition with 10 permutation matrices necessarily contains P a , P b , . . . , P j . It is easy to see that there is a unique solution to the equality M x = a of the MILP with x t = 0 for t / ∈ {a, b, . . . , j}, as the submatrix M containing only the corresponding 10 columns has full column rank.

B = 1023 • I O O A , A =       3
Any generalized Birkhoff heuristic obtains at least 11 permutation matrices for the matrix A of the proof of Lemma 1. One such decomposition is shown in a tabular format in Fig. 1. In Fig. 1a, we write the matrix of the lemma explicitly for convenience. Then in Fig. 1b, we give a BvN decomposition. The column headers (the first line) in the table contains the coefficients of the permutation matrices. The nonzero column indices of the permutation matrices are stored by rows. For example, the first permutation has the coefficient 129 and columns 3, 5, 2, 1, and 4 are matched to the rows 1-5. The bold indices signify entries whose values are equivalent to the coefficients of the corresponding permutation matrices, at the time where the permutation matrices are found. Therefore, the

A = 513 •    0 0 0 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0   + 257 •    0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0   + 127 •    0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0   + 63 •    0 0 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0   + 31 •    0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0   + 15 •    0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0   + 7 •    0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1   + 3 •    0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 1 0 0 0 0   + 2 •    0 0 0 0 1 0 0 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0 0   + 2 •    0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1   + 2 •    1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 1 0   + 1 •    1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 1    Figure 2:
The output of the Greedy heuristic for the matrix given in the proof of Lemma 1.

corresponding entries become zero after the corresponding step. For example, in the first permutation, a 5,4 = 129.

The output of the Greedy heuristic for the matrix A is given in Fig. 2 for reference. It contains 12 permutation matrices.

Analysis of two known heuristics for computing BvN decompositions

Empirical results in the literature demonstrated that the Greedy heuristic [4] can obtain a much smaller number of permutation matrices compared to the Birkhoff heuristic. Here, we compare these two heuristics theoretically. First we show that the original Birkhoff heuristic does not have any constant ratio approximation guarantee. Furthermore, for an n × n matrix, its worst-case approximation ratio is Ω(n).

We begin with a small example shown in Fig. 3. We decompose a 6 × 6 matrix A (0) which has an optimal BvN decomposition with three permutation matrices; the main diagonal, the one containing the entries equal to 4, and the one containing the remaining entries. For simplicity, we used integer values in our example. However, since the row and column sums of A (0) is equal to 6, it can be converted to a doubly stochastic matrix by dividing all the entries to 6. Instead of the optimal decomposition, in the figure, we obtain the permutation matrices as the original Birkhoff heuristic does. Each red-colored entry set is a permutation and contains the minimum possible value, 1.

In the following, we show how to generalize the idea for having matrices of arbitrarily large size, with three permutation matrices in an optimal decomposition, while the original Birkhoff heuristic obtains n permutation matrices.

Lemma 2. The worst-case approximation ratio of the Birkhoff heuristic is Ω(n).

Proof. For any given integer n ≥ 3, we show that there is a matrix of size n × n whose optimal BvN decomposition has 3 permutations, whereas the Birkhoff heuristic obtains a BvN decomposition with exactly n permutation matrices. The example in Fig. 3 is a special case for n = 6 for the following construction process.

A (0) =         1 4 1 0 0 0 0 1 4 1 0 0 0 0 1 4 1 0 0 0 0 1 4 1 1 0 0 0 1 4 4 1 0 0 0 1         A (1) =         0 4 1 0 0 0 0 1 3 1 0 0 0 0 1 3 1 0 0 0 0 1 3 1 1 0 0 0 1 3 4 0 0 0 0 1         A (2) =         0 4 0 0 0 0 0 0 3 1 0 0 0 0 1 2 1 0 0 0 0 1 2 1 1 0 0 0 1 2 3 0 0 0 0 1         A (3) =         0 3 0 0 0 0 0 0 3 0 0 0 0 0 0 2 1 0 0 0 0 1 1 1 1 0 0 0 1 1 2 0 0 0 0 1         A (4) =         0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 1 1 1 0 0 0 1 0 1 0 0 0 0 1         A (5) =         0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 1        
Let f : {1, 2, . . . , n} → {1, 2, . . . , n} be the function f (x) = (x mod n) + 1. Given a matrix M, let M = F(M) be another matrix containing the same set of entries where the function is f (•) is used on the coordinate indices to redistribute the entries of M on M . That is m i,j = m f (i),f (j) . Since f (•) is oneto-one and onto, if M is a permutation matrix then F(M) is also a permutation matrix. We will start with a permutation matrix, and run it through F for n -1 times to obtain n permutation matrices, which are all different. By adding these permutation matrices, we will obtain a matrix A whose optimal BvN decomposition has three permutation matrices, while the n permutation matrices used to create A correspond to a decomposition that can be obtained by the Birkhoff heuristic.

Let P 1 be the permutation matrix whose ones, which are partitioned into three sets, are at the positions

1st set (1, 1), 2nd set (n, 2) , 3rd set
(2, 3), [START_REF] Brualdi | Notes on the Birkhoff algorithm for doubly stochastic matrices[END_REF][START_REF] Dufossé | Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices[END_REF]

, . . . , (n -1, n) . ( 9
)
Let us use F(•) to generate a matrix sequence P i = F(P i-1 ) for 2 ≤ i ≤ n. For example, P 2 's nonzeros are at the positions

1st set (2, 2), 2nd set (1, 3) , 3rd set (3, 4), (4, 5)
, . . . , (n, 1) .

We then add the P i s to build the matrix

A = P 1 + P 2 + • • • + P n .
We have the following observations about the nonzero elements of A:

1. a i,i = 1 for all i = 1, . . . , n, and only P i has a one at the position (i, i).

These elements are from the first set of positions of the permutation matrices, as identified in [START_REF] Dolan | The NEOS Server 4.0 administrative guide[END_REF]. When put together, these n entries form a permutation matrix P (1) . 2. a i,j = 1 for all i = 1, . . . , n and j = ((i+1) mod n)+1, and only P h , where h = (i mod n) + 1, has a one at the position (i, j). These elements are from the second set of positions of the permutation matrices, as identified in [START_REF] Dolan | The NEOS Server 4.0 administrative guide[END_REF]. When put together, these n entries form a permutation matrix P (2) . 3. a i,j = n -2 for all i = 1, . . . , n and j = (i mod n) + 1, where all P for ∈ {1, . . . , n} \ {i, j} have a one at the position a i,j . These elements are from the third set of positions of the permutation matrices, as identified in [START_REF] Dolan | The NEOS Server 4.0 administrative guide[END_REF]. When put together, these n entries form a permutation matrix P (3) multiplied by the scalar (n -2).

In other words, we can write A = P (1) + P (2) + (n -2) • P (3) , and see that A has a BvN decomposition with three permutation matrices. We note that each row and column of A contains three nonzeros; and hence three is the smallest number of permutation matrices in a BvN decomposition of A.

Since the minimum element in A is 1, and each P i contains one such element, the Birkhoff heuristic can obtain a decomposition using P i for i = 1, . . . , n. Therefore, the Birkhoff heuristic's approximation is no better than n 3 , which can be made arbitrarily large.

We note that Greedy will optimally decompose the matrix A used in the proof above.

We now analyze the performance of the Greedy heuristic. We start with a few experiments which indicate that there is a connection between the performance of Greedy and the values in the matrix.

We create a set of n × n matrices. To do that, we first fix a set of z permutation matrices {C 1 , . . . , C z } of size n × n. These permutation matrices with varying values of α will be used to generate the matrices. The matrices are parametrized by the subscript i and each A i is created as follows:

A i = α 1 • C 1 + α 2 • C 2 + • • • + α z • C z
where each α j for j = 1, . . . , z is a Table 1: Experiments showing the dependence of the performance of Greedy on the values of the matrix elements. n is the matrix size; i ∈ {10, 20, 30, 40, 50} is the parameter for creating matrices using α j ∈ [1, 2 i ]; z is the number of permutation matrices used in creating A i . Five experiments for a given pair of n and i. Greedy obtains k i permutation matrices. The average and the maximum number of permutation matrices obtained by Greedy for five random instances are given. k i/z is a lower bound to the performance of Greedy, as z ≥ Opt. ], and we also set a randomly chosen α j equivalent to 2 i to guarantee the existence of at least one large value even in the unlikely case that all other values are not large enough. As can be seen, each A i has the same structure and differs from the rest only in the values of α j 's that are chosen. As a consequence, they all can be decomposed by the same set of permutation matrices. We present our results in two sets of experiments shown in Table 1, for two different n. In both cases we create five random A i for i ∈ {10, 20, 30, 40, 50}, that is we have five matrices with the parameter i, and there are five different i. We have n = 30 and z = 20 in the first set, and n = 200 and z = 100 in the second set. Let k i be the number of permutation matrices Greedy obtains for a given A i . The table gives the average and the maximum k i of five different A i , for given n and i. By construction, each A i has a BvN with z permutation matrices. Since z is no smaller than the optimal value, the ratio ki z gives a lower bound on the performance of Greedy. As seen from the experiments, as i increases, the performance of Greedy gets increasingly worse. This shows that a constant ratio worst case approximation of Greedy is unlikely. While the performance depend on z (for example, for small z, Greedy is likely to obtain near optimal decompositions), it seems that the size of the matrix does not largely affect the relative performance of Greedy. Now we attempt to explain the above results theoretically. Lemma 3. Let α 1 P 1 + • • • + α k P k be a BvN decomposition of a given doubly stochastic matrix A with the smallest number k of permutation matrices. Then, for any BvN decomposition of A with ≥ k permutation matrices, we have ≤ k • maxi α i mini αi . If the coefficients are integers (e.g., when A is a matrix with constant row and column sums of integral values), we have ≤ k • max i α i .

Proof. Consider a BvN decomposition α

1 P 1 + • • • + α P with ≥ k . Assume without loss of generality that α 1 ≥ • • • ≥ α k and α 1 ≥ • • • ≥ α .
We know that the coefficients of these two decompositions sum up to the same value. That is

i=1 α i = k i=1 α i .
Since α is the smallest of α, and α 1 is the largest of α , we have

• α ≤ k • α 1 , and hence k ≤ α 1 α .
By assuming integer values, we see that α ≥ 1 and thus

≤ k • max i α i .
This lemma evaluates the approximation guarantee of a given BvN decomposition. It does not seem very useful, because of the fact that even if we have min i α i , we do do not have max i α i . Luckily, we can say more in the case of Greedy.

Corollary 1. Let k be the smallest number of permutation matrices in a BvN decomposition of a given doubly stochastic matrix A. Let α 1 and α be the first and last coefficients obtained by the Greedy heuristic for decomposing A. Then, ≤ k • α1 α . Proof. This is easy to see as Greedy obtains the coefficients in a non-increasing order [4, Lemma 3], and α 1 ≥ α j for all 1 ≤ j ≤ k for any BvN decomposition containing α j . Lemma 3 and Corollary 1 give a posteriori estimates of the performance of the Greedy heuristic, in that one looks at the decomposition and tells how good it is. This potentially can reveal a good performance. For example, when Greedy obtains a BvN decomposition with all coefficients equivalent, then we know that it is an optimal BvN. The same cannot be told for the Birkhoff heuristic though (consider the example proceeding Lemma 2). We also note that the ratio given in Corollary 1 should usually be much larger than the practical performance.

Conclusion

We have investigated heuristics for obtaining Birkhoff-von Neumann decomposition of doubly stochastic matrices and presented three results. First, there are matrices whose decompositions with the smallest number of permutation matrices cannot be obtained by any Birkhoff-like heuristic. Second, the worstcase approximation ratio of the original Birkhoff heuristic is Ω(n). Third, the performance of the Greedy heuristic depends on the values of matrix elements, and a bound using the first and the last coefficients found by Greedy is obtained.

The shown bound for the performance of Greedy is expected to be much larger than what one observes in practice, as the bound can even be larger than the upper bound on the number of permutation matrices. A tighter analysis should be possible to explain the practical performance of the Greedy heuristic (which was demonstrated earlier [START_REF] Dufossé | Notes on Birkhoff-von Neumann decomposition of doubly stochastic matrices[END_REF]).

Figure 1 :

 1 Figure 1: The matrix A from Lemma 1, and a BvN decomposition with 11 permutation matrices which can be obtained by a generalized Birkhoff heuristic. Each column in (b) corresponds to a permutation, where the first line gives the associated coefficient, and the following lines give the column indices matched to the rows 1 to 5 of A.

Figure 3 :

 3 Figure 3: A sample matrix to show that the original Birkhoff heuristic can obtain BvN decomposition with n permutation matrices while the optimum one has 3.