Anne Benoit

Loris Marchal

Yves Robert

Oliver Sinnen
email: o.sinnen@auckland.ac.nz

Mapping pipelined applications with replication to increase throughput and reliability

Mapping and scheduling an application onto the processors of a parallel system is a difficult problem. This is true when performance is the only objective, but becomes worse when a second optimization criterion like reliability is involved. In this paper we investigate the problem of mapping an application consisting of several consecutive stages, i.e., a pipeline, onto heterogeneous processors, while considering both the performance, measured as throughput, and the reliability. The mechanism of replication, which refers to the mapping of an application stage onto more than one processor, can be used to increase throughput but also to increase reliability. Finding the right replication trade-off plays a pivotal role for this bi-criteria optimization problem. Our formal model includes heterogeneous processors, both in terms of execution speed as well as in terms of reliability. We study the complexity of the various subproblems and show how a solution can be obtained for the polynomial cases. For the general NP-hard problem, heuristics are presented and experimentally evaluated. We further propose the design of an exact algorithm based on A* state space search which allows us to evaluate the performance of our heuristics for small problem instances.

platform found in industry and academia, well suited for applications with negligible or low inter-task communication volume.

We focus in this paper on the pipeline skeleton, i.e., an application with linear dependencies between stages, to be executed on such clusters. The pipeline skeleton is typical of streaming applications like video and audio encoding and decoding applications, etc. [6, 18, 10, 19, 20]. Each stage has its own communication and computation requirements: it reads an input file from the previous stage, processes the data and outputs a result to the next stage. For each data set, initial data is input to the first stage, and final results are output from the last stage. The pipeline workflow operates in synchronous mode: after some latency due to the initialization delay, a new task is completed every period.

Key metrics for a given workflow are the throughput and the reliability. The throughput measures the aggregate rate of processing data, and it is the rate at which data sets can enter the system. Equivalently, the inverse of the throughput, defined as the period, is the time interval required between the beginning of the execution of two consecutive data sets. The reliability of the application is the probability that all computations will be successful. The period can be minimized by dealing different data sets to different processors in a round-robin fashion, while reliability is increased by replicating computations on a set of processors (i.e., computing each data set several times). The application fails to be executed only if all processors involved in a same redundant computation fail during execution. In this paper, we focus on bi-criteria approaches, i.e., minimizing the failure probability under period constraints, or the converse. Thus we must investigate the trade-off of using processors to deal data sets and reduce the period, or to replicate computations and increase the reliability. The optimization problem can be stated informally as follows: which stage to assign to which processor?

We target applications and platforms on which the cost of communications is negligible in comparison to the computation cost. SpeedHom platforms are made of identical speed processors, while SpeedHet platforms consist of different speed processors. In terms of reliability (defined as the inverse of the overall failure probability), we consider either processors with identical failure probabilities, denoted FailureHom, or with different failure probabilities FailureHet. All four combinations of processor speed and reliability are addressed in this paper.

.

, the replication is such than a given number of processor failures is supported. In [2], we defined the reliability of a pipelined application similarly than in the present

Introduction

Mapping applications onto parallel platforms is a difficult challenge. Several scheduling and load-balancing techniques have been developed for homogeneous architectures (see [START_REF] Shirazi | Scheduling and load balancing in parallel and distributed systems[END_REF] for a survey) but the advent of heterogeneous clusters has rendered the mapping problem even more difficult. Typically, such clusters are composed of different-speed processors interconnected by some communication network. They constitute a common paper, and used replication to decrease the failure probability of an application.

To the best of our knowledge, this paper presents the first study of problems in which both replication for performance and replication for reliability are used simultaneously, and thus considering the two antagonist criteria, throughput and reliability. First we formally detail the model and the optimization problems in Section 2. Then in Section 3 we establish the complexity of all problem instances, and derive some new approximation results. Bi-criteria problems are NP-hard in most cases because of the antagonistic criteria. In particular, it is NP-hard to minimize the failure probability, given a bound on the period of the application, when targeting heterogeneous platforms. Thus, we introduce polynomial time heuristics in Section 4 to propose practical solutions to this problem. One contribution is the design of an exact algorithm based on A* (see Section 5), which allows us to evaluate the absolute performance of our heuristics on small problem instances in Section 6. Finally, we state some concluding remarks in Section 7.

Framework

We outline in this section the characteristics of the applicative framework, as well as the model for the target platform. Next we detail the bi-criteria optimization problem.

Applicative framework

We consider a pipeline of n stages S i , 1 ≤ i ≤ n, as illustrated on Fig. 1. Input data is fed into the pipeline and processed from stage to stage, until output data exits the pipeline after the last stage.

The i-th stage S i receives an input from the previous stage, performs a number of w i computations, and outputs data to the next stage. Communication costs are assumed to be negligible in comparison with the computation costs, thus the time required to compute one data set for stage S i is proportional to w i . S n ... ...

w 1 w 2 w n S 1 S 2 S i w i
Figure 1: The application pipeline.

Target platform

We target a platform with p processors P u , 1 ≤ u ≤ p, fully interconnected as a (virtual) clique. The speed of processor P u is denoted as s u , and it takes X/s u time-units for P u to execute X floating point operations. We consider either SpeedHom platforms made of identical processors: s u = s for 1 ≤ u ≤ p, or general SpeedHet platforms. Moreover, we assume that communications are negligible compared to computation costs, thus we do not characterize here the network bandwidths and communication times.

In order to evaluate and optimize the reliability of the mappings produced, we associate a failure probability 0 < f u < 1 with each processor P u , 1 ≤ u ≤ p. It is the probability that P u fails during the execution of the application. We consider constant failure probabilities as we are dealing with pipelined applications. These applications are meant to run during a very long time, and therefore we address the question of whether the processor will fail or not at any time during execution. It might seem odd that this probability is completely independent of the execution time, since one may believe that the longer a processor executes, the larger the chance that it fails. However, we target a steady-state execution of the application, for instance in a scenario in which we would loan/rent resources. Computers could be suddenly reclaimed by their owners, as during an episode of cycle-stealing [START_REF] Awerbuch | Making commitments in the face of uncertainty: how to pick a winner almost every time[END_REF][START_REF] Bhatt | On optimal strategies for cycle-stealing in networks of workstations[END_REF][START_REF] Rosenberg | Optimal schedules for cycle-stealing in a network of workstations with a bag-of-tasks workload[END_REF]. The failure probability should thus be seen as a global indicator of the reliability of a processor. Also note that we consider in this work only fail-silent (a faulty processor does not produce any output) and fail-stop (no processor recovery) processor failures. We do not consider link failures since in a grid framework, a different path can be found to overcome such failures (and we consider communication costs as negligible).

A platform composed of processors with identical failure probabilities is denoted Fail-ureHom and otherwise FailureHet. Note that it seems natural to consider FailureHet platforms when processors have different speeds, thus for SpeedHet platforms, while Speed-Hom platforms are more likely to be FailureHom. However, we consider all combinations in the complexity study in Section 3.

Mapping problem

The general mapping problem consists of assigning application stages to platform processors. Instead of only considering the simple, but restrictive mapping of one stage per processor, we consider so called interval mapping, where several consecutive stages S i , S i + 1, . . . , S j -1, S j , i.e., an interval, are mapped together. Such mappings have been extensively studied, see [START_REF] Subhlok | Optimal mapping of sequences of data parallel tasks[END_REF][START_REF] Subhlok | Optimal latency-throughput tradeoffs for data parallel pipelines[END_REF][START_REF] Benoit | Mapping pipeline skeletons onto heterogeneous platforms[END_REF][START_REF] Benoit | Complexity results for throughput and latency optimization of replicated and data-parallel workflows[END_REF].

The cost model associated with interval mappings is the following. We search for a partition of [1..n] into m ≤ p intervals I j = [d j , e j ] such that d j ≤ e j for 1 ≤ j ≤ m, d 1 = 1, d j+1 = e j +1 for 1 ≤ j ≤ m-1 and e m = n. Interval I j is mapped onto a set of processors, A j . These processors are organized into l j teams (l j ≤ |A j |), and each team is in charge of one round of the deal. We stress that processors within a team perform redundant computations, while different teams assigned to the same interval execute distinct data sets in a round-robin fashion. Note that we could envision that faster teams would perform more computations, i.e., operate on more data sets, than slower ones. However, enforcing a round-robin execution is the key to preventing out-of-order execution, thereby allowing for a reasonably simple orchestration of all operations. We assume that a processor cannot participate in two different teams, thus for all 1 ≤ j, j ′ ≤ m, A j ∩ A j ′ = ∅. The period is then expressed as:

P = max 1≤j≤m e j i=d j w i l j × min Pu∈A j s u (1)
Indeed, for each interval, the computation is slowed down by the slowest processor enrolled for this interval (because data sets are distributed in round robin to each team), but the period of the slowest processor is divided by l j , since this processor gets only one data set every over l j ones.

Given a partition of processors into l teams T k , where l = m j=1 l j and 1 ≤ k ≤ l, the failure probability of the application is computed as follows, since the computation is successful if at least one processor per team does not break down during execution:

F = 1 - 1≤k≤l (1 - u∈T k f u ) (2)
The optimization problem is to determine the best interval mapping, over all possible partitions into intervals, and over all processor assignments. The objective can be to minimize either the period, or the failure probability, or a combination: given a threshold period, what is the minimum failure probability that can be achieved? and the counterpart: given a threshold failure probability, what is the minimum period that can be achieved?

3 Complexity Results

Mono-criterion problems

All mono-criterion problems have been studied previously, but we recall the results in this section since they introduce the bi-criteria complexity results. Failure probability [START_REF] Benoit | Optimizing latency and reliability of pipeline workflow applications[END_REF]. Finding the interval mapping which minimizes the failure probability F is easy on any kind of platforms: the value of Equation ( 2) is minimized with l = 1 and T 1 = {P 1 , . . . , P p }. Thus, all stages are grouped as a single interval, and one single team consisting of all p processors is processing this interval. Period [START_REF] Benoit | Complexity results for throughput and latency optimization of replicated and data-parallel workflows[END_REF]. In order to minimize the period, we first note that there is no need to add several processors in a team since these extra processors only aim at increasing the reliability.

On a SpeedHom platform, P is minimized by grouping all stages as a single interval and having p teams working on this interval, with one single processor per team. Since all speeds are identical and equal to s, in Equation (1), we have min s u = s thus no computational power is lost due to differences in processor speeds.

However, on SpeedHet platforms, if two processors of different speeds process the same interval, then the slowest one becomes the bottleneck and the other one is not fully used. In this case, we derived a sophisticated polynomial time algorithm for homogeneous applications with w i = w for 1 ≤ i ≤ n. The idea consists in first proving that there exists an optimal solution which replicates intervals of stages onto intervals of consecutive speed processors. Then, a binary search and a dynamic programming approach are combined to compute which intervals of processors should be used onto which stages (see [START_REF] Benoit | Complexity results for throughput and latency optimization of replicated and data-parallel workflows[END_REF]). The problem turns out NP-hard in the general case of heterogeneous pipelines, since it is no longer possible to consider only intervals of processors. To establish the completeness, we use a reduction from NUMERICAL 3D MATCHING [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF], which is NP-complete in the strong sense (see [START_REF] Benoit | Complexity results for throughput and latency optimization of replicated and data-parallel workflows[END_REF]).

Bi-criteria problems

To the best of our knowledge, it is the first time that the bi-criteria problem considering period and failure probability minimization is studied. Since period minimization is NP-hard on a SpeedHet platform, all bi-criteria problems are NP-hard on such platforms. First we derive complexity results for SpeedHom platforms, and we provide a polynomial time algorithm for FailureHom platforms. Then we prove that the problem becomes NP-hard when considering FailureHet platforms. Approximation results are provided in Section 3.3. Preliminary result for SpeedHom platforms. We prove that there is an optimal solution which consists of a single interval of stages if we restrict to SpeedHom platforms with identical speed processors.

Lemma 1. For SpeedHom platforms, there exists a mapping of the pipeline as a single interval of stages which minimizes the failure probability (resp. period) under a fixed period (resp. failure probability) threshold.

Proof. The proof works by contradiction. Suppose that there is an optimal mapping which consists of several interval of stages I 1 , . . . , I m . Processors are partitioned into l teams T 1 , . . . , T l , and we use the notations introduced in Section 2.3. The period obtained by this solution is denoted P opt , and for 1 ≤ j ≤ m, we have e j i=d j w i ≤ P opt × l j × s. We build a new solution consisting of a single interval of stages which are processed by the same l teams as in the optimal solution. According to Equation (2), the value of the failure probability remains the same, since the teams are identical. The period of the new solution is P = P n i=1 w i l×s . If P > P opt , then for 1 ≤ j ≤ m, e j i=d j w i < P × l j × s, and if we sum this equation for all j, we obtain n i=1 w i < P × l × s = n i=1 w i . Thus, P ≤ P opt , and the new solution consisting of a single interval of stages is optimal.

Note that this result does not hold on SpeedHet platforms. Consider for instance the following problem with two stages and two processors such that w 1 = s 1 = 1 and w 2 = s 2 = 10. The threshold on the failure probability is fixed to 1, thus it cannot be violated. In this case, we can obtain a period 1 with two intervals while we cannot do better than 11/10 with a single interval. With identical failure probabilities. Building upon Lemma 1, we are able to derive a polynomial algorithm to solve the bi-criteria problem when further restricting to FailureHom platforms.

Theorem 1. For SpeedHom-FailureHom platforms, the optimal mapping which minimizes the failure probability under a fixed period threshold can be determined in polynomial time O(p); the optimal mapping which minimizes the period under a fixed failure probability threshold can be determined in polynomial time O(p log(p)).

Proof. Consider first the case with a fixed value of the period, P * . According to Lemma 1, we restrict the search to solutions consisting in a single interval of stages. We can thus compute the minimum number of teams required to achieve the threshold period: l min = P n i=1 w i P * ×s . According to Equation (2), a higher value of l can only increase the failure probability: if l > l min teams, we can obtain a lower failure probability be removing ll min teams. Moreover, the minimum of the function is reached by greedily assigning all processors to teams in order to have balanced teams:

∀1 ≤ k, k ′ ≤ l, -1 ≤ |T k | -|T k ′ | ≤ 1.
To see this, assume that the optimal solution involves two teams with x and y ≥ x + 2 processors respectively. Using two teams of x + 1 and y -1 processors leads to a smaller failure probability. Indeed, (1

-f x+1 )(1 -f y-1 ) > (1 -f x )(1 -f y ) ⇔ f x+1 + f y-1 < f x + f y ⇔ f y-x-1 (1 -f ) < 1 -f .
The processor assignment can thus be done in O(p), where p is the number of processors.

For the converse problem, we cannot easily derive the optimal number of teams from the failure probability threshold F * . Rather, we try all p possible number of teams. With 1 ≤ l ≤ p teams, the period is P = P n i=1 w i l×s , which is increasing when l decreases. Starting with l = p, we use the previous algorithm to compute the corresponding failure probability, and check if this failure probability is lower than F * . We stop as soon as we get a failure probability under the threshold. The overall complexity is thus O(p 2 ) since we may run p times the previous algorithm which works itself in O(p). We could also perform a binary search on the number of teams, thus reducing this complexity to O(p log(p)).

With different failure probabilities. In the previous algorithm for FailureHom platforms, we exploit the fact that balancing the number of processors into teams minimizes the failure probability. This property does not hold with different failure probabilities. Indeed, we prove that the problem becomes NP-had in this case, even though Lemma 1 still holds.

Theorem 2. For SpeedHom-FailureHet platforms, the problem of finding a mapping which respects both a fixed threshold period and a fixed threshold failure probability is NP-complete.

Proof. The problem clearly belongs to the class NP: given a solution, it is easy to verify that it is an interval mapping, that the teams are partitioning the set of processors, and to compute its period and reliability with Equations ( 1) and ( 2) in polynomial time.

To establish the completeness, we use a reduction from 3-PARTITION, which is NPcomplete in the strong sense [START_REF] Garey | Computers and Intractability, a Guide to the Theory of NP-Completeness[END_REF]. We consider an arbitrary instance I 1 of SP: given 3n positive integer numbers {a 1 , . . . , a 3n } and a bound B, assuming that B 4 < a i < B 2 for all i and that 3n i=1 a i = nB, are there n subsets I 1 , I 2 , . . . , I n such that

I 1 ∪ I 2 . . . ∪ I n = {1, 2, . . . , 3n}, I j ∩ I j ′ = ∅ if j = j ′ ,

and

i∈I j a i = B for all j (and |I j | = 3 for all j). Because 3-PARTITION is NP-complete in the strong sense, we can encode the 3n numbers a i in unary and assume that the size of I 1 is O(n + M ), where M = max i {a i } (or equivalently O(n + B)).

We build the following instance I 2 of our problem: the application has only one pipeline stage with w = n, and the SpeedHom-FailureHet platform consists of 3n processors with speeds s u = 1 and failure probabilities f u = 2 -au , for 1 ≤ u ≤ 3n. The threshold period is fixed to P = 1, and the threshold reliability is fixed to R = (1 -2 -B ) n . We ask whether there exists a mapping which respects both thresholds. The size of I 2 is polynomial in the size of I 1 : the f u can be encoded in binary with O(M ) bits, and the bound R can be encoded in binary with O(nB) bits.

We now show that I 1 has a solution if and only if I 2 has a solution. Suppose first that I 1 has a solution, I 1 , . . . , I n . We build the following solution for I 2 : n teams are working on the unique stage, with processors in team T j being the three processors P u such that u ∈ I j , for 1 ≤ j ≤ n. Since there are n teams, the period of this mapping is n/n = 1 = P, and the reliability is:

n j=1 (1 - u∈T j f u ) = n j=1 (1 -2 - P u∈T j au ) = n j=1 (1 -2 -B ) = R.
Thus we have a solution to I 2 .

Conversely, if I 2 has a solution with l teams T 1 , . . . , T l , then the period is n/l ≤ P = 1, and so l ≥ n. If l > n, we can build a solution with a period 1 and a higher reliability by removing ln teams, so we consider in the following that l = n. The reliability is

n j=1 (1 -2 - P u∈T j
au ). This quantity is always strictly larger than n j=1 (1 -2 -B ) unless when u∈T j a u = B for 1 ≤ j ≤ n (see Lemma 2 below, where U j = u∈T j a u and C = 2 -nB , hence (1-C 1/n ) n = R). Thus, the processor indices of the mapping correspond to a solution of I 1 , which concludes the proof.

Lemma 2. Given n real numbers U j , 0 < U j < 1, whose product C = n j=1 U j is given, the quantity R = n j=1 (1 -U j ) has a unique maximum which is obtained when all the U j have the same value C 1/n . Proof. The proof is by induction. For n = 2, R = (1 -U 1 )(1 -C U 1 ) = 1 -U 1 -C U 1 + C, which has a unique maximum when U 1 = C U 1
, hence the result. Assume the result is true for n -1 variables and consider the case with n variables. We have

R = n-1 j=1 (1 -U j ) (1 -U n ). By induction the quantity n-1 j=1 (1 -U j ) is maximum when the first n -1 variables U j , whose product is C Un , have the same value C Un 1 n-1 . We derive that R = 1 -C Un 1 n-1 n-1 (1 -U n ).
Considering R as a function R(U n ) and differentiating, we obtain

R ′ (U n ) =   (n -1) 1 - C U n 1 n-1 n-2 1 n -1 C 1 n-1 U -( 1 n-1 )-1 n (1 -U n )   -1 - C U n 1 n-1 n-1 .
We simplify

R ′ (U n ) = 1 - C U n 1 n-1 n-2 C 1 n-1 U -( 1 n-1 )-1 n (1 -U n ) -1 + C 1 n-1 U -( 1 n-1 ) n , which has a unique zero when C 1 n-1 U -( 1 n-1 )-1 n = 1, hence when U n = C 1 n . The value of U j for j < n is thus C Un 1 n-1 = C 1 n
, which concludes the proof.

Approximation results

According to Lemma 1, it is always efficient to use only one single interval for SpeedHom platforms. In this section, we revisit the (mono-criterion) period minimization problem for SpeedHet platforms. This problem is NP-hard but the following theorem allows us to compare the single interval solution with the optimal one. Unfortunately, there are cases for which the period of the optimal solution with one single interval is m times greater than the period of the optimal solution with m intervals.

Lemma 3. The optimal single-interval mapping for period minimization on a SpeedHet platform with p processors can be found in time O(p log p).

Proof. Simply sort the processors by non-increasing speed: s 1 ≥ s 2 ≥ . . . s p . According to Equation (1), if we use the fastest i processors, we obtain a period P = W i×s i , where W = n i=1 is the total work. Therefore the optimal solution is to select the value of i, 1 ≤ i ≤ p, such that i × s i is maximal and to use the fastest i processors.

Theorem 3. For a workflow with n stages, single interval mapping is a n-approximation algorithm for the period minimization problem on SpeedHet platforms. Moreover, this approximation factor cannot be improved.

Proof. Consider that the optimal solution consists in m ≤ n intervals, with l j processors in charge of interval I j and s j being the speed of the slowest processor for this interval. The period P m of this mapping is such that e j i=d j w i ≤ l j s j P m . We can build a solution with one single interval which is an m-approximation of the solution: let k be one of the interval for which the product l k s k is maximal, i.e., l j s j ≤ l k s k for 1 ≤ j ≤ m. Then we use only processors in charge of this interval to compute the whole pipeline. The achieved period is:

P 1 = n i=1 w i l k s k = m j=1 ej i=dj w i l k s k ≤ m j=1 ej i=dj w i l j s j ≤ mP m .
Thus, P 1 ≤ P m , which proves the approximation result.

To prove that we cannot obtain a better approximation result, we consider the following problem instance with n pipeline stages. Let K be a constant, which can be arbitrarily large For 1 ≤ i ≤ n, we have w i = 1, and l i = K i-1 processors of speed s i = 1/K i-1 . The optimal solution has n intervals of size 1, and allocates stage S i to the l i processors of speed s i . The optimal period is P = 1. There is no idle time in the previous mapping, which shows its optimality. Furthermore, there is no other mapping without idle time, which shows the uniqueness of the optimal solution.

With a single interval mapping, we need to decide which processors to use. If we use the first j sets of processors, which are the fastest ones, we obtain a period

P j = n i=1 w i ( j k=1 l k ) min j k=1 s k = n ( j k=1 l k )s j = nK j-1 1 + K + • • • + K j-1 = nK j-1 (K -1) K j -1 = n(1 + o(K)).
Therefore, the optimal single-interval mapping, whose period is the minimum of the P j , has a period which is arbitrarily close to n for K large enough. This concludes the proof.

The proof of Theorem 3 actually shows that single-interval mapping is a m-approximation for the period, where m is the number of intervals in the optimal solution (but of course m is unknown a priori, and we can only bound it by n, the number of stages). However, for bi-criteria problems, we stress that is always good to reduce the number of intervals in order to reduce the failure probability. Hence it looks very challenging to derive approximation algorithms for such problems on SpeedHet-FailureHet platforms.

Heuristics

In this section we propose heuristics for the general mapping problem on processors that are heterogeneous both in terms of computation speed, SpeedHet, and in terms of reliability, FailureHet. We concentrate on the optimization problem of minimizing the failure probability under a fixed period, because this is the more frequently encountered scenario. The opposite problem of minimizing the period under a fixed failure probability can be tackled by doing a simple binary search over the prescribed period using the proposed algorithms. We will see this in the evaluation of Section 6.

From Lemma 1 we know that for SpeedHom platforms, an optimal mapping can be found using a single-interval mapping. For this reason, we start with the design of a heuristic that employs only one interval. We expect this approach to be reasonable even for slightly speed-heterogeneous processors. Further, it will serve as a building block for our multiple-interval heuristic which we discuss subsequently.

Single-interval mapping

The heuristic proposed here, called OneInterval, minimizes the failure probability F under a fixed upper period bound P on a SpeedHet-FailureHet platform. As our approach is to group all stages into a single interval, the heuristic has to determine the number of teams and how the processors are assigned to these teams. The idea is to test all possible numbers of teams, i.e., from 1 to p. For each fixed number of teams l, we discard all processors that would not allow to achieve the given period P. This procedure is motivated by Lemma 3 for finding the optimal mapping for period minimization. The remaining processors are assigned to the l teams, so that the total failure probability is minimized.

From Lemma 2 we know that to minimize the total failure probability, we need to minimize the maximum team failure probability among all teams. Note that this problem, similar to a partition problem, remains NP-hard. However, we use a simple greedy heuristic which is known to be efficient: we consider processors by non-decreasing failure probability and assign the next processor to the team with the highest failure probability. This heuristic, detailed in Algorithm 1, has a time complexity of O(p 2 log p).

Multiple-interval mapping

Theorem 3 demonstrates that even an optimal single-interval solution can be far from the absolute optimal on a SpeedHet platform. Hence we propose in this section a multipleinterval mapping heuristic, called MultiInterval. Again the heuristic minimizes the failure probability F under a fixed upper period bound P on a SpeedHet-FailureHet platform.

The proof of Theorem 3 shows that it can be beneficial for the minimal period to have multiple intervals. Our proposed heuristic therefore starts with min(n, p) intervals. Either we create one interval per stage, or we create at most p intervals, starting from the first stage and adding stages to the current interval while its computation cost is lower than

P n i=1 w i p
(step 1). The processors are then distributed over these intervals in such a way that the maximum ratio of interval computation weight to accumulated processor speed is Algorithm 1: Heuristic OneInterval optimizing F with fixed P using single interval for l = 1 to p do Discard processors whose speed is less than minimum speed necessary to achieve period with l teams: procs = {P u , 1 ≤ u ≤ p : s u ≥ P n i=1 w i l×P } Order processors procs by non-decreasing failure probability f u into list L Create l teams T k , 1 ≤ k ≤ l, set failure probability of all teams to f T k = 1 foreach P u in L do Find team T max with highest failure probability Assign P u to T max ; f Tmax ← f Tmax × f u Compute the total failure probability F l = 1 -l k=1 (1 -T k ) Choose among the p previous solutions the one with the lowest total failure probability F l minimized: we greedily add processors, sorted by non-increasing speed, to intervals with the highest current ratio (step 2). Teams are formed and processors are allocated in each interval using the single-interval Algorithm 1, including the utilization of processors that were unused in the allocation of previous intervals (step 3). If no solution can be found for an interval, we increase the bound on the period for this interval until Algorithm 1 returns a valid allocation. Then, if the period bound is not achieved for at least one interval, we merge the interval with the largest period, either with the previous or with the next interval, until the bound on the period is respected (step 4). The actual optimization of the failure probability happens in the last step (step 5), in which we merge intervals with highest failure probability as long as it is beneficial. Indeed, we know from Section 3.1 that the fewer intervals and teams the (potentially) better is the failure probability. In other words, Equation (2) for the failure probability F is minimized for l = 1. Note that Algorithm 1 is called each time we try to merge intervals (steps 4 and 5).

This heuristic, described in has a time complexity of O(p 3 log p).

Optimal algorithm using A*

In order to find optimal solutions for small instances of our problem, we employ an A* best-first state space search algorithm [START_REF] Dechter | Generalized best-first search strategies and the optimality of A*[END_REF][START_REF] Russell | Artificial Intelligence: A Modern Approach[END_REF]. A* has successfully been employed before for optimal scheduling, e.g. [START_REF] Shahul | Optimal scheduling of task graphs on parallel systems[END_REF]. For the given problem it is especially interesting, as the non-linear nature of the failure probability rules out other exact approaches, e.g. Let w(I) be the computation amount of interval I let speed(I) be the sum of all processor allocated to interval I (zero for the moment) ✄ 2. Assign processors to intervals foreach each P u , 1 ≤ u ≤ p, in non-decreasing speed s u order do Assign P u to interval I max with highest ratio w(I max )/s(max ) ✄ 3. Build Team in each interval for each interval I i in non-decreasing w(I i ) s(I i ) order do Apply Algorithm 1 for I i and its allocated processors plus unused processors of previous intervals If Algorithm 1 is unable to find a solution within the period bound, increase this bound and re-run Algorithm 1, until we have a valid solution Let P(I) = w(I)/(n I × s minI ), where n I is the number of teams chosen for interval I and s minI the minimum speed of processors allocated to interval I ✄ 4. Merge intervals to reach period bound while the period bound is not reached, and there are at least two intervals do Find interval I i with highest period P(I i ) Merge I i with I i-1 if it exists; apply Algorithm 1 on this interval Merge I i with I i+1 if it exists; apply Algorithm 1 on this interval Accept the merging which results in smallest period for the current interval ✄ 5. Merge intervals to decrease failure probability while we can decrease the failure probability, and there are at least two intervals do Find interval I i with highest failure probability F I i Merge I i with I i-1 if it exists; apply Algorithm 1 on this interval Merge I i with I i+1 if it exists; apply Algorithm 1 on this interval Accept the merging which results in better remaining highest failure probability, provided that the period bound is satisfied Return the current solution, provided it satisfies the period bound the OPEN list, so that the list remains sorted by non-decreasing c value. The algorithm terminates when the state s that is taken out of the OPEN list is a goal state, i.e., a complete solution. With a good cost estimation function c, A* can find a complete mapping which is guaranteed to be optimal (with certain conditions on the function c(s)) without exploring the entire solution space. The determining factor in the effectiveness of A* is the quality of the underestimated c(s) values, i.e., how tight they are to the real costs.

A CLOSED set of the states that have already been visited, which is used in classical A*, is not necessary to maintain, because our state space construction does not create duplicates.

Solution space

The first step in the design of A* is to describe the solution space of the problem and the successive construction of larger partial solutions. In order to find the optimal solution for our general mapping problem, the grouping into intervals and the mapping of the processors is integrated into one solution space. The proposed approach consists in sequentially assigning stages to intervals, and processors to teams for these intervals. We start with an initial state, which is an empty state, where no stage nor processor has been assigned yet. A state can be expanded by considering the first stage S j which is not yet allocated. We have two choices when adding this stage: (i) Stage S j is included in the previous interval, and inherits from the team structure: only one new state is added to the solution space. (ii) Stage S j starts a new interval, without any processor assigned to it at the moment, which forms the new state. Then, the state expansion process continues by assigning the not yet allocated processors, i.e., the free processors, to the current interval. For each free processor P u , we have three choices, hence (potential) new states:

1. P u is not allocated to the current interval; 2. P u is allocated to an existing team for the current interval; 3. P u is allocated to a new team for the current interval. Some invalid states may be created this way, but are immediately discarded, such as a state where no processor at all has been allocated to an interval. All valid states are added in the sorted OPEN list of states. Finally, a goal state is a state where all stages have been included into intervals, and the process of allocating processors to the last interval is complete.

Fig. 2 shows the state tree after expanding all states, for a small example of an application with two stages, on two processors. For this example, the expansion reaches six goal states and four invalid states.

Underestimate functions

The c(s) function is the central part of an A* algorithm. It is an underestimate of the exact minimum cost c * (s) of any goal state that is based on the state s. If the function c(s) fulfills c(s) ≤ c * (s) for any state s, it is called admissible. With an admissible c(s) function, A* is guaranteed to find an optimal solution. The number of examined states

Legend [S1] [S1] P1 [S1] P1, P2 [S1] (P1) [S1, S2] (P1) [S1] [S2] (P1) [S1] (P1, P2) [S2] [S1] (P1) (P2) [S2] [S1, S2] (P1) (P2) [S1] (P1) (P2) (P1, P2) [S1, S2] [S1, S2] [S2] [S1] [S2] [S1] P1 [S2] [S1] [S1] (P1, P2) [S1] P2 (P1) [S2] [S1] P2 [S2] [S1] empty state [S1] (P2) P1 (P1) (P1) (P2) (P2) (P2) (P1) (P2) (P2)
Figure 2: Complete state tree for two-stage application mapped on two processors. E.g., state at bottom right corner is made of a single interval [S 1 , S 2 ] and two teams with one processor each: (P 1 ) and (P 2 ). depends on how close c(s) is to c * (s). In general, the more accurate, the fewer states have to be examined.

In our problem, we have to design two such underestimate functions, both for F and P. Like in the previous section, we concentrate on the problem of minimizing the failure probability F with a given period bound P. Thus, we first need an underestimate function f for the failure probability: for a partial solution s, f (s) is an underestimate of the minimum failure probability of any complete solution based on s. We also have to check that a solution never exceeds the period bound P. Therefore, we also use an underestimate function p for the period: p(s) is an underestimate of the minimum period of any complete solution based on the partial state s.

Estimate of the failure probability

To understand the underestimation of the failure probability of a state s, recall that the failure probability does not depend on the structure of intervals, but only on how processors are grouped into teams, Equation [START_REF] Benoit | Optimizing latency and reliability of pipeline workflow applications[END_REF]. Adding a new team to a partial mapping can never reduce the failure probability, only increase it.

The underestimation must therefore assume that all remaining stages are added to the currently open interval (having more than this last interval would at least require one additional team). It is straight forward to calculate the failure probabilities of the already closed teams of the closed intervals. The difficulty arises for underestimating the failure probabilities of the teams of the last interval. It must be assumed that all free processors are assigned to the already existing teams of the last interval (if there is no team yet, all free processors are assigned to a single team). Assigning the free processors to the existing teams in an optimal way is an NP-hard problem though.

From Lemma 2 we know that the best failure probability would be achieved if we perfectly balanced the failure probability among the teams. Thus, we make the ideal-izing assumption that processors are perfectly divisible, and that we have an "amount of reliability" (corresponding to the failure probability of a single team made of all free processors). This amount is then distributed to the existing teams, by starting with the teams having the worst failure probability, to minimize the maximum failure probability. The result of this procedure gives us an underestimate of the best achievable failure probability for the open teams. Together with the failure probabilities of the closed teams we obtain an underestimate for the considered partial state. The complete process is described in Algorithm 3.

Algorithm 3: Estimate of the failure probability for incomplete teams Input: K teams with T 1 , . . . , T K , with failure probabilities f team T 1 , . . . , f team T K , and p available processors with failure probabilities f 1 , . . . , f p Compute the aggregated available failure probability F ← p u=1 f u ; we consider now that we have a divisible processor with failure probability F . while F < 1 do Select the X teams T max 1 , . . . , T max X with maximum failure probability f team max Compute the second maximum failure probability f team max2 (the maximum among all other teams) ρ ← max

f team max2 f team max , F 1/X foreach team T i in T max 1 , . . . , T max X do
Decrease the failure probaility of T i by ρ:

f team i ← f team i × ρ Update the available failure probability: F ← F/ρ X

Estimate of the period

Our estimation of the period takes into account the stages already structured into intervals, and the processors already allocated to teams for these intervals. We simply use Equation (1) on those closed intervals. In order to better estimate the period of a final solution expanded from the current state, we also compute the minimum period that can be achieved with the other stages and processors. As outlined in Section 3.1, the mono-criterion problem of period minimization is already NP-hard in the heterogeneous case. However, we can derive a simple lower bound on the period, by assuming a perfect load-balancing of processor speeds: P ≤ P S i w i P Pu su . We use this bound, with all stages that do not belong to the closed intervals, and all processors that are not allocated to a closed interval, to derive a second underestimate of the achievable period. The final underestimate of the period is then the maximum of the previous two estimations.

Optimization

A classical optimization for A* is employed. We first compute the solution given by a heuristic (in practice, the single-interval heuristic OneInterval). During A*, all states that have an underestimated failure probability larger than the one obtained from the heuristic can be immediately discarded. This does not create less states (as those states would not have been expanded) but can save queue memory which is crucial for A*.

Evaluation results

In this section, we present the simulations we have undertaken to evaluate the performance of the two proposed heuristics and the A* algorithm.

Simulation settings

We have randomly generated several workload scenarios. Stage computational weights are uniformly distributed in the interval [1; 10]. Processor speeds are uniformly distributed in the interval [1; 10], while processor failure probabilities are uniformly distributed in the interval [0.1; 0.9]. The number of stages randomly varies from 5 to 20, just as the number of processors. However, for A*, we had to limit both numbers to 10: for larger scenarios, the memory demand of A* is too high. Finally, all tests were performed on a quad-processor machine (64-bit AMD Operon at 2.2GHz) with 32Gbytes of RAM.

Results

We first give the results of both heuristics compared to A*. Figure 3(a) shows the behavior of the different algorithms on a given scenario (with 8 processors and 8 stages). We see that the heuristics approach the optimal (minimal) failure probability computed by A*. For a period bound between 3 and 7, MultiInterval achieves a better performance than OneInterval: the use of several intervals is then crucial to get the optimal performance. On the contrary, for a period bound around 2, OneInterval is closer to the optimal: because of its initial construction of the interval, and of the choice of the intervals to merge, MultiInterval sometimes fails to find the best one-interval solution. Figure 3(b) shows the running times of the algorithms: the heuristics are extremely fast, whereas A* need a few seconds to run on this example: when the period is neither too small nor too large (say in the interval [1.5; 4.5]), A* is confronted with a large solution space, hence it performs many operations. However, the main limitation of A* is the memory needed to store all partial states generated during the execution; this is why we we had to limit our search with A* to smaller scenarios.

Figures 4(a) and 4(b) show the distribution of the ratio between the failure obtained by a heuristic, and the optimal failure computed with A*: when this ratio is close to one, the performance of the heuristic is good. The box at 10 accounts for all scenarios where the heuristic does not find a solution, whereas A* finds one. We can see that both heuristics perform similarly. On average in the cases where the heuristics find a solution, their failure probability is 20% above the optimal failure probability. In very few cases (a couple of scenarios, representing less than 0.1% of all cases), the failure probability can be up to 10 times the optimal one. However, the heuristic sometimes fail to find a solution within the period bound (whereas A* exhibits a solution): this happens in 9.6% of the cases for OneInterval, and 11.5% for MultiInterval. In the complete set of experiments (with larger scenarios), OneInterval is better than MultiInterval in 61% of the cases, and MultiInterval is better than OneInterval in 20% of the cases (in the remaining 19% cases, either they both fail, or they find the same solution). On average, OneInterval gives a failure probability which is 2% higher than MultiInterval.

We also compared OneInterval with the optimal solution restricted to one interval: it is easy to limit the search space of A* to solution using only one interval. This shows that OneInterval is very good, with an average failure probability only 0.05% larger than the optimal one (and 5% larger in the worst case). This proves that the problem restricted to one interval is not very difficult, but the difficult part is to find how to split stages into intervals, and how to allocate processors to these intervals.

Conclusion

In this paper, we have investigated the problem of mapping pipelined applications onto heterogeneous platforms, while maximizing both the throughput and the reliability of the application. The main challenge was to perform a good trade-off between reliability and performance in the use of processors. Given that the complexity of monocriterion problems had already been established, we presented new results concerning bi-criteria optimization. In particular, we derived a polynomial-time algorithm for Speed-Hom-FailureHom platforms, while we proved that the problem is NP-hard as soon as one level of heterogeneity is introduced. Moreover, we provided approximation results to compare the optimal mapping consisting of a single interval with any other solution.

To address the practical solution of these difficult problems, we proposed polynomialtime heuristics for the most general problem instance SpeedHet-FailureHet. Strongly inspired by the theoretical study, the main heuristic is built upon a simpler heuristic which proposes a single interval mapping. Because of the strong non-linearity of the failure probability function, we were not able to formulate an integer linear program to compute the optimal solution for small problem instances. Rather, we investigated the use of A*, and proposed non-trivial underestimate functions leading to an efficient execution of this (exponential-time) algorithm. The experimental results demonstrated the very good behavior of the heuristics, which failed only in 10% of the tests, and which otherwise returned a failure probability only 20% worse than the optimal one.

Future work includes the investigation of further approximation results, enhanced multi-interval heuristics, and improved cost estimation and pruning techniques for the exact A* based algorithm.

Figure 3 :

 3 Figure 3: Behavior on an example scenario.

Figure 4 :

 4 Figure 4: Performance of the heuristics compared to the optimal solution

  Integer Linear Programming. A* is a generic algorithm that performs a best-first search over all possible solutions. Every node in the search space is a state s that represents a partial solution and has an underestimated cost value c(s) associated with it. This value represents the minimum cost of any complete solution based on the partial solution represented by s. For our mapping problem a state s is a partial mapping and the cost value c(s) is either the failure probability or the period. A* maintains an OPEN list of states, populated with the initial state in the beginning. At each step, A* expands

	Algorithm 2: Heuristic MultiInterval optimizing F with fixed P using multiple
	intervals
	✄ 1. Make initial intervals
	if n ≤ p then
	Create n intervals, one for each stage
	else
	Create p intervals, each with approximate computation weight of	P n i=1 w i p

the state with the lowest c(s) value, i.e., it creates new states that are larger partial solutions. A c value is calculated for each of these new states and they are inserted into


Sabine Barrat
email: barrat@loria.fr

Salvatore Tabbone
email: tabbone@loria.fr

A Bayesian network for combining descriptors: application to symbol recognition

Keywords: Symbol recognition, Descriptor combination, Variable selection, Probabilistic graphical models, Bayesian networks

In this paper, we propose a descriptor combination method, which enables to improve significantly the recognition rate compared to the recognition rates obtained by each descriptor. This approach is based on a probabilistic graphical model. This model also enables to handle both discrete and continuous-valued variables. In fact, in order to improve the recognition rate, we have combined two kinds of features: discrete features (corresponding to shape measures) and continuous features (corresponding to shape descriptors). In order to solve the dimensionality problem due to the large dimension of visual features, we have adapted a variable selection method. Experimental results, obtained in a supervised learning context, on noisy and occluded symbols, show the feasibility of the approach.

Introduction

Pattern recognition applications have to face the problem of describing a large number of different objects for recognition. A recognition system should be robust to variability (geometric transformations, noise, occlusions, …) and to scalability, when a large number of classes and images should be recognized. Symbol recognition is a field within pattern recognition for which a lot of efforts have already been made [START_REF] Lladós | Symbol recognition by error-tolerant subgraph matching between region adjacency graphs[END_REF][START_REF] Valveny | Graphics recognition. Recent advances and new opportunities[END_REF][START_REF] Yang | Symbol recognition via statistical integration of pixellevel constraint histograms: a new descriptor[END_REF][START_REF] Zhang | Symbol recognition with kernel density matching[END_REF]. Symbol recognition is usually decomposed and learning (factorization of parameters fitting, computation of probability distributions …). We propose, in this paper, an original method of descriptor combination applied to symbol classification. Thus, we have adapted the probabilistic graphical model theory to the symbol recognition problem. In this model, continuous and discrete variables are combined. Continuous variables correspond to shape descriptors, and the discrete ones correspond to shape measures. Thanks to this combination, the proposed classifier is more robust to deformations and to the size of database, when the number of symbols increases. The originality of our approach also relies on the use of a variable selection method [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], to overcome the dimensionality problem related to the size of feature vectors and the inherent network complexity.

The organization of the paper is as follows. In Sect. 2, the main properties of a Bayesian network-based classifier are introduced and lead to the presentation of our probabilistic model for symbol recognition. The visual features used to represent the symbols are described in Sect. 3. The feature selection algorithm which allows us to increase the recognition rate by focusing only on the main features while reducing the dimensionality problem is also explained in Sect. 4. Our method is evaluated on a database of noisy and occluded symbols (Sect. 5). Finally, Sect. 6 brings conclusions and opens new perspectives to our work.

Representation and classification of images

Context and objectives

Our work is focused on symbol recognition by combining descriptors. Given an image database, where each image contains one symbol, we try to recognize the "perfect" symbol (the model) represented in each image. In fact, the symbols contained in the images are not perfect: they can be noisy, deformed and can have occlusions. This recognition problem can be view as a classification problem: our aim is to assign each image to the class corresponding to the perfect symbol (the model) of this image. However, no perfect symbol is available. Therefore, we cannot just minimize a distance between each image of the database and each perfect symbol. On the other hand, this classification task can be resolved by using a supervised learning method, from a subset of the database where the class label (the perfect symbol) is known for each image.

Moreover, in order to describe all kinds of shapes properly, even deformed or noisy shapes, and thus increase the recognition rate, our proposition is to combine several shape descriptors. Now, shape descriptors can provide vectors of continuous or discrete values: let f j be a query image characterized by a set of features F composed of:

-m continuous visual features, denoted v 1 , . . . , v m , -n discrete visual features, denoted D F 1, . . . , D F n.

The chosen visual features are issued from 3 shape descriptors and 3 shape measures. Shape descriptors provide vectors of continuous values, and each shape measure provides a single discrete value.

Consequently, it seems appropriate to propose a classifier that enables to manage both discrete and continuous features. Although most classification methods handle only discrete data and thus require a pre-processing step of discretization in order to transform each continuous-valued variable into a discrete one, few classification methods can handle both discrete and continuous-valued variables. It is the case of Support Vector Machines [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF], Random Forests [START_REF] Breiman | Random forests[END_REF], and Bayesian classifiers [START_REF] Friedman | Bayesian network classifiers[END_REF]. Support Vector Machines (SVM) and Random Forests (RF) are well known for their ability to handle highdimensional data. On the contrary, Bayesian classifiers are sensitive to the dimensionality of the data, but they often perform well in many domains. Therefore, we have chosen to construct a Bayesian classifier for its ability to combine discrete and continuous-valued variables. Moreover, we show that this Bayesian classifier, associated with a variable selection method, is competitive with SVM, even on high-dimensional data.

Bayesian classifiers

Let I be a new image designed by a particular instance f = { f 1 , . . . , f n } of the feature vector F = {F 1 , . . . , F n }. Our aim is to assign I to a class c i among k classes. Each c i is a particular instance of the variable C. The Naïve Bayes (N B) is a simple probabilistic classification algorithm that often performs well in many domains. This classifier encodes a distribution P N B (F 1 , . . . , F n , C), from a given training set (composed of labeled data). The resulting probabilistic model can be used to classify the new instance I . In fact, the Bayes rule is applied to compute the probability of c i given the particular instance f . Then the classifier based on N B returns the label c i , i ∈ {1, . . . , k}, that maximizes the posterior probability P i = P N B (c i | f 1 , . . . , f n ), where:

P i = P N B ( f 1 , . . . , f n |c i ) × P N B (c i ) P N B ( f 1 , . . . , f n ) and P N B ( f 1 , . . . , f n ) = k j=1 P N B ( f 1 , . . . , f n |c j ) × P N B (c j ).
However, we are interested in the probability distributions of discrete and continuous features and their conditional dependence relations. Let us consider each component of continuous vectors (issued from shape descriptors) as a continuous random variable and the discrete values (provided by shape measures) as discrete variables. This model is too big to be represented as a unique joint probability distribu-tion. Therefore, it is required to introduce some sparse and structural a priori knowledge: the Naïve Bayes has to be extended to take into account continuous and discrete variables. In this perspective, the probabilistic graphical models, and especially Bayesian networks, are a good way to solve this kind of problem. In fact, within Bayesian networks, the joint probability distribution is replaced by a sparse representation only among the variables directly influencing one another. Interactions among indirectly related variables are then computed by propagating inference through a graph of these direct connections. Consequently, Bayesian networks are a simple way to represent a joint probability distribution over a set of random variables, to visualize the conditional properties and to compute complex operations like probability learning and inference, according to graph-based computations.

Bayesian networks

Definitions

Formally, a Bayesian network for a set of random variables V (continuous or/and discrete) is a pair B = G, . The first component, G, is a directed acyclic graph whose vertices correspond to random variables V 1 , . . . , V n , and whose edges represent direct dependencies between variables. The graph G encodes independence assumptions: each variable V i is independent of its non-descendants given its parents in G. The second component of the pair, , represents the set of parameters that quantifies the network. It contains a parameter θ v i |Pa(v i ) = P B (v i |Pa(v i )) for each possible value v i of V i , and Pa(v i ) of Pa(V i ), where Pa(V i ) denotes the set of parents of V i in G. That is, the Bayesian network, in its initial state, contains the initial a priori probabilities of each node of the network: P B (v i |Pa(v i )). Thanks to the conditional independence assumption of each variable given its parents, the joint probability distribution P B (V 1 , . . . , V n ) can be reduced to this formula:

P B (V 1 , . . . , V n ) = n i=1 P B (V i |Pa(V i )) = n i=1 θ v i |Pa(v i )
The framework of Bayesian networks comes with many algorithms for performing inference (i.e., the computation of posteriors probabilities) and learning (factorization of parameters fitting, computation of probability distributions,…). The algorithms we used in this work are briefly described below.

Parameter learning

Only one has a description of a model, knowing the structure of the graph and probabilistic forms for each variables, one wants to estimate the numerical values of each parameter.

Let assume we have either discrete or continuous variables (or a mix of them), and, for the simple case, a set of data describing many possible cases for each variables. The data set can either be complete or have missing data. In each case, a different solution will be used. In the case, where the data set has no missing values, an approach is to consider the parameters having the highest probabilities to generate the most similar data set if the Bayesian network was used to draw random values according to the probability distribution it describes (hence the name "generative model"). This method is known as the Maximum Likelihood. Let call D the data set, then P(d|M) is the probability of a data d ∈ D to be generated by the model M and is called the likelihood of M given d. Therefore, the likelihood of M given the full data set D is:

L(M|D) = P(D|M) = d∈D P(d|M)
For the sake of computational simplicity (or to help deriving an analytic form), the log-likelihood is often used:

L(M|D) = d∈D log 2 P(d|M)
Therefore, the principle of maximum likelihood prefers to choose parameters with the highest likelihood:

θ = argmax θ L(M θ |D)
In general, a maximum likelihood is obtained by counting the frequencies over the total count. Whenever you do not have enough data for each case (missing data) one of the most popular algorithm is the Expectation-Maximization (EM) algorithm. The general purpose of this algorithm is explained in detail in [START_REF] Dempster | Maximum likelihood from incomplete data via the em algorithm[END_REF]. During the Expectation phase, the data set is locally completed, then a Maximization step is performed to find the current maximum likelihood estimate (as seen above) using the completed data. In a Bayesian network, the first step of the EM algorithm can be easily done using a map a posteriori algorithm, i.e., computing the most probable values of the missing data variables given other known variables. The second step is then executed and can either be done with an optimization algorithm if no analytical form of the maximum likelihood is known, or with the previous approach. These two steps are repeated until convergence. The algorithm is initialized with random probability distribution parameters.

The EM algorithm is also used to learn the parameters of Gaussian distributions, which are considered as missing data.

Inference

An inference algorithm is necessary to compute the posterior probability distributions of unobserved nodes. According to the Bayesian network topology, the inference process propagates the values from the leaf level to the inferred node. Many algorithms can be used [START_REF]Learning in Graphical Models[END_REF]. The most popular is the message passing algorithm [START_REF] Kim | A computational model for combined causal and diagnostic reasoning in inference systems[END_REF]. In this technique, each node is associated with a processor, which can send some messages to its neighbors, in an asynchronous way, until it reaches stability.

Bayesian network classifiers

Bayesian networks can be used as classifiers. For example, the Naïve Bayes can be represented by the structure in Fig. 1, where:

-C refers to the class variable, -F 1 , . . . , F n are the feature variables.

The Naïve Bayes is a simple and efficient model, but it requires discrete variables. Since, we have to manage continuous (values provided by shape descriptors) and discrete (values provided by shape measures) variables, this model has to be extended to take into account continuous and discrete variables.

A Gaussian-Mixtures and Bernoulli Mixture model

A Bayesian network classifier, which handles both discrete and continuous-valued variables, is proposed. We present a hierarchical probabilistic model, the Gaussian-Mixtures and Bernoulli Mixture model, in order to classify large databases of symbols. In fact, the observation of some peaks on the different histograms of the vector components provided by shape descriptors has led us to consider that the continuous visual features can be estimated by mixtures of Gaussian densities. The discrete variables have a Bernoulli distribution. In fact, these variables can take two values: 1 if the corresponding shape measure provides a value smaller than 0.5, or else 2. Finally, the proposed model is inspired by the Naïve Bayes. Indeed, the class variable is connected to each other. Now let F be the training set composed of m instances f j = { f j 1 , . . . , f j n }, ∀ j ∈ {1, . . . , m}, ∀i ∈ {1, . . . , n}, where n is the dimension of the signatures provided by the concatenation of the feature vectors issued from the computation of all the descriptors for each image on the training Gaussian Class Component GMM Fig. 2 A Probabilistic graphical model as GMMs set. Each instance f j , ∀ j ∈ {1, . . . , m} is then characterized by n continuous variables. A supervised classification is considered; then, F instances are divided into k classes c 1 , . . . , c k . Let G 1 , . . . , G g be g groups whose each has a Gaussian density with a mean µ l , ∀l ∈ {1, . . . , g} and a covariance matrix l . Besides, let π 1 , . . . , π g be the proportions of the different groups, θ l = (µ l , l ) be the parameter of each Gaussian and = (π 1 , . . . , π g , θ 1 , . . . , θ g ) the global mixture parameter. The probability density of F conditionally to the class c i , ∀i ∈ {1, . . . , k} can be defined by

P i ( f, ) = g l=1 π l p( f, θ l )
where p( f, θ l ) is the multivariate Gaussian defined by the parameter θ l .

We have one Gaussian Mixture model per class, which can be represented by the probabilistic graphical model in Fig. 2, where:

-The "Class" node is a discrete node, which can take k values corresponding to the pre-defined classes c 1 , . . . , c k . -The "Component" node is a discrete node which corresponds to the components (i.e., the groups G 1 , . . . , G g ) of the mixtures. This variable can take g values, i.e., the number of Gaussians used to compute the mixtures. It is an hidden variable which represents the weight of each group (i.e., the π l , ∀l ∈ {1, . . . , g}). -The "Gaussian" node is a continuous variable which represents each Gaussian G l , ∀ ∈ {l = 1, . . . , g} with its own parameter (θ l = (µ l , l )). It corresponds to the set of feature vectors in each class. The model can be completed by the discrete variables, denoted D F 1, . . . , D F n, where n is the number of shape measures, and D F i represents the value of each shape measure. Dirichlet priors [START_REF] Robert | A Decision-Theoretic Motivation[END_REF] have been used for the probability estimation of the variables D F 1, . . . , D F n. That is we introduce additional pseudo counts at every instance in order to ensure that they are all "virtually" represented in the training set. Therefore every instance, even if it is not represented in the training set, will have a not null probability. Like the continuous variables, the discrete variables corresponding to the discrete measures are included in the graphical model by connecting them to the class variable. Now our classifier can be depicted by the Fig. 3. The hidden variable "α" shows that a Dirichlet prior is used. The box around the variable D F denotes n repetitions of D F for each shape measure.

This Bayesian classifier means that continuous and shape features, representing images, are assumed to have been generated conditional on the same class. Therefore, the resulting Bernoulli and Gaussian mixture parameters should correspond: concretely if an image, represented by continuous visual descriptors, has an high probability under a certain class, then its discrete shape measures should have an high probability under the same class.

In order to classify a query image f j , the class node C is inferred thanks to the massage passing algorithm. This image, characterized by its continuous shape features v j 1 , . . . , v j m and its discrete shape features D F 1 j , . . . , D F k j is considered as an "evidence" represented by:

P( f j ) = P(v j 1 , . . . , v j m , D F 1 j , . . . , D F n j ) = 1
when the network is evaluated. Thanks to the inference algorithm, the probabilities of each node are updated in function of this evidence. After the belief propagation, we know, ∀i ∈ {1, . . . , k}, the posterior probability:

P(c i | f j ) = P(c i |v j 1 , . . . , v j m , D F 1 j , . . . , D F n j )
The query f j is assigned to the class c i which maximizes this probability.

Symbol description

This section explains how we have adapted the theoretical method before mentioned to a symbol recognition problem. We present the visual features we used. The set of chosen features is composed of 3 different off-the-shelf shape descriptors and 3 shape measures. The choice of these features is not really important because, the aim of this paper is to show that combining shape features improves the symbol classification, whatever the used features. The distinction between shape descriptors and shape measures is determined by the size of the features: we consider single value features as shape measures and feature vectors as shape descriptors. Moreover, shape measures are discretized with a discretization threshold fixed at 0.5. This discretization has sense with this kind of features, because each of them is composed of a single value normalized between 0 and 1. Thus, a shape measure has an intrinsic meaning. Finally, this discretization enables to consider shape measures as discrete variables and thus to show the interest of discrete and continuous features combination for shape recognition.

In this perspective, we have chosen three pixel shape descriptors: the Generic Fourier Descriptor (GFD), the Zernike descriptor, and the R-signature 1D and three shape measures: compactness, rectangularity and ellipticity. We briefly present each descriptor and measures below.

Shape descriptors

Generic Fourier descriptor: Generic Fourier descriptor is based on Fourier transform [START_REF] Zhang | Shape-based image retrieval using general fourier descriptor[END_REF]. The rotation invariance is achieved by using the modified polar Fourier transform (MPFT), and the scaling invariance is achieved after normalization.

Zernike descriptor: Zernike descriptor [START_REF] Kim | A modified zernike moment shape descriptor invariant to translation, rotation and scale for similarity-based image retrieval[END_REF] is a descriptor based on Zernike moments. Zernike moments of a given shape are calculated as correlation values of the shape with Zernike basis functions, in that all the pixels of the shape, independently of their position, contribute with the same weight to the Zernike moments. These moments are rotation invariant. To make the Zernike moments of the shape descriptor invariant also to translation and scaling, a given shape is normalized, by obtaining the smallest circle centered at the center of mass, covering all the shape pixels. Then the obtained circle is adjusted to match the radius of Zernike moment basis functions. The Zernike shape descriptor consists of low-order magnitudes of Zernike moments. R-signature 1D: The R-signature 1D [START_REF] Tabbone | Technical symbols recognition using the two-dimensional radon transform[END_REF] uses Radon transform to represent an image. The Radon transform is the projection of an image in a particular plan. This projection has interesting properties. According to these geometrical properties, a 1D signature of the transform is created. This signature checks the properties of invariance to some geometrical transformations, such as the translation and the scaling (after normalization). The rotation invariance is achieved by a cyclic permutation of the signature, or directly from its Fourier transform.

Shape measures

Compactness:

The compactness measure C represents the ratio of the shape area to the area of a circle (the most compact shape) having the same perimeter:

C = 4π A P 2
where P is the perimeter and A, the area. This measure is invariant to translation, rotation, and scaling.

Rectangularity:

The rectangularity degree [START_REF] Rosin | Measuring rectangularity[END_REF] R is equal to the ratio of the shape area to the area of its minimal bounding box:

R = A L * l
where A is the shape area and L (respectively l) is the length (respectively the width) of the minimal bounding box.

Ellipticity:

The ellipticity degree is obtained from the ratio of the major axis to the minor axis [START_REF] Teague | Image analysis via the general theory of moments[END_REF]:

= 1 - b a
where a is the major axis and b the minor axis. This measure is invariant to rotation, translation, and homothety.

Dimensionality reduction

Only the n, n ∈ {1, 2, 3} continuous descriptors we want to combine are computed on each image, we dispose of n signatures per image. The concatenation of these signatures provides us a new feature vector per image. The large dimensions of initial visual signatures and their concatenation imply a dimensionality problem. In fact, a too large feature dimension increases the computation time and causes a wrong Gaussian mixture learning, because of the Small Sample Size (SSS) problem: there is a disproportion between the training set size and the feature vector dimension. To overcome this problem, we have used a dimensionality reduction method. A lot of methods have been proposed in the literature to reduce the dimension of vectors [START_REF] Denoeux | Dimensionality reduction and visualization of interval and fuzzy data: a survey[END_REF][START_REF] Piccardi | Maximum-likelihood dimensionality reduction in gaussian mixture models with an application to object classification[END_REF]. Among dimensionality reduction methods, we consider especially feature selection methods, because they enable to reduce dimension by selecting a subset of initial features, on the contrary to the methods which reduce dimension by providing new variables issued from initial variable combinations [START_REF] Duda | Pattern Classification. 2nd edn[END_REF]. Methods of variable selection are then more suitable to our problem, because our aim is to reduce the number of features, in order to reduce the size of our Bayesian network and our method complexity.

The most popular methods of variable selection are heuristics based on sequential runs, which consist in iteratively adding or removing variables [START_REF] Devijver | Pattern Recognition: A Statistical Approach[END_REF]. In these approaches, it is possible to begin with an empty set of variables and to add variables to the variables which are already selected (it is the Sequential Forward Selection (SFS) [START_REF] Pudil | Floating search methods in feature selection[END_REF]), or to begin with the set of all variables et to remove variables in this set (it is the Sequential Backward Selection (SBS)). These methods are known for their simplicity and their rapidity. However, they are known for their instability too. Moreover, since they do not explore all possible subsets of variables and they do not enable to come back during the process, they are not optimal. Thus, we have chosen a feature selection method since it enables to extract from the feature vectors, just the most relevant and discriminating features, with a minimal information loss. The regression method LASSO (Least Absolute Shrinkage and Selection Operator) [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF] has been used for its stability and implementation efficiency. Moreover, the LASSO method especially enables to select variables and takes into account the class variable values to select a subset of variables. In Sect. 5, we compare the results obtained on our database, by the LASSO and the SFS method which stay one of the most popular.

The principle of LASSO is to shrink the regression coefficients by imposing a penalty on their size (we speak of shrinkage methods too). These coefficients minimize a penalized residual sum of squares:

β lasso = arg min β N i=1 (y i -β 0 - p j=1 x i j β j ) 2 subject to p j=1 |β j | ≤ s.
The linear form of the LASSO has been applied in a pre-processing stage, totally independent of our Bayesian classifier, on our visual features. For each training set, y i represents the sum of the mean vector features of the class c i , and x j = {x j 1 , x j 2 , . . . , x j p } the p features of the instance j. Then, just the subset of the selected variables is used in our model.

The LASSO uses a L1 penalty: p j=1 |β j |. This constraint implies that for small values of s (s ≥ 0), some of the coefficients β will be null. So choosing s is like choosing the number of predictors in a regression model. Therefore, the variables corresponding to the coefficients different from zero are selected.

The LASSO solutions have been computed by the Least Angle Regression (LAR) procedure [START_REF] Efron | Least angle regression[END_REF]. This algorithm exploits the special structure of the LASSO problem and provides an efficient way to compute the solutions simultaneously for all values of s.

Thus, we can distinguish two independent phases in our method: first, the LASSO is used to select the most relevant visual features. Secondly, our Gaussian-Mixtures and Bernoulli Mixture model is used to classify new images represented by the pre-selected continuous visual features and the three discrete shape measures.

Experimental results

We have used the symbols of GREC database [START_REF] Valveny | Symbol recognition contest : a synthesis[END_REF] for our tests (see Fig. 4), especially created for the symbol recognition contest GREC'2005.

This database is mainly defined from two application domains, architecture and electronic, because these symbols are most largely used by graphic recognition teams and represent a great number of different forms. We have 50 different symbol models for which we have applied some noises based on Kanungo [START_REF] Kanungo | A statistical, nonparametric methodology for document degradation model validation[END_REF] model. These noises are similar to noise obtained when a document is scanned, printed, or photocopied. Moreover, we have applied to these symbols some rotations of different degrees and different zooms, in order to obtain a database of 3, 600 images, constituted of 72 different images per model. We have evaluated our method by performing a crossvalidation by using 75% of the database for the training and the remaining 25% for the tests. The tests are repeated four times in order to use each database instance for the training and the tests. The recognition rate is obtained by taking the mean recognition rate of the 4 tests.

Since, we want to improve the recognition rate by combining descriptors and selecting variables, we limit ourselves to the experiments comparing:

-the classification after variable selection with the LASSO vs. the classification without automatic variable selection and with a well-known variable selection method, -the classification by combining 2 or 3 continuous descriptors vs. the classification with only one continuous descriptor, -the classification by combining discrete and continuous features vs. the classification with continuous descriptors only, -the classification by combining 3 continuous descriptors with our method vs. two state-of-art classifiers.

First of all, we can remark that the reducing the size of the descriptors improves the recognition rate for all the classifiers. Moreover, the variable selection with the LASSO method has enabled us to significantly reduce the number of variables. Table 1 shows the mean number of variables selected for each descriptor with the LASSO method compared to the well-known SFS method [START_REF] Pudil | Floating search methods in feature selection[END_REF]. We can see that the LASSO method enables to select fewer variables than the SFS method (see Table 1). Table 2 shows the mean recognition rates in function of different variable selection methods, by combining the 3 available descriptors with our Gaussian-Mixtures and Bernoulli Mixture model (GM-B) and two state-of-art classifiers: a classical SVM classifier [START_REF] Chang | LIBSVM: a library for support vector machines[END_REF] and the fuzzy k-nearest neighbor (FKNN) [START_REF] Keller | A fuzzy k-nearest neighbor algorithm[END_REF]. The recognition rates for these three classifiers without variable selection and after a variable subset selection with the SFS and LASSO methods, or with some random selections, are compared. For the random selection, the number of variables chosen randomly is set to the one obtained with the LASSO. The FKNN has been computed with k = 1 and with k = m where m is the mean number of images per class in the training set. The results in Table 2 show that the variable selection with 123 Let us consider Table 3. The notation G (respectively Z and R) means that the GFD descriptor (respectively the Zernike descriptor and the R-signature 1D) has been used. Finally, the "+" operator indicates that we have combined the descriptors. The recognition rates confirm that combining 2 or 3 descriptors performs always better than any of them alone. In fact, we observe that the combination of 2 descriptors increases the recognition rate by 18% on average compared to the use of only one descriptor. Besides we can notice that the combination of the 3 descriptors is better, by 18.3% on average, to use just one of them. Moreover, even if we obtain a high recognition rate with Zernike descriptor, this rate will not decrease when we combine this descriptor with one or two other descriptors, whatever these descriptors and even if the added descriptors have a low rate (it is the case with the R-signature 1D). The bad behavior of a descriptor does not impede the other descriptor behaviors.

Finally, the last line of the Table 2 shows the effectiveness of our approach compared to the SVM classifier and the FKNN. The results have been obtained by combining the 3 descriptors and after using the variable selection method LASSO. It appears that the proposed GM-B results are always better than the ones of SVM and FKNN.

The initial database of 3,600 instances has been extended to a database of 5,400 instances by randomly generating occlusions on the half image set of each class from the initial database. In fact, we can meet this kind of degradation when we have to segment a graphical document where symbols are embedded into the graphic or are partially occluded by text for example. The generated occlusions are from different sizes, and their locations in the images have been chosen randomly. Now we have a larger and more distorted database, composed of 108 instances per class. For example, the Fig. 5 presents 4 symbol models (first column) and 5 occluded images disturbed by some noises derived from these models. Our classifier has been applied after variable subset selection with the LASSO. This time, our method has been evaluated by performing 3 cross-validations whose each proportion of the training set is 25,50 and 75% of the database, the remaining, respectively, 75,50, and 25% are hold for testing set. In each case, the tests are repeated 10 times in order that each database instance would be used for the training and the test. For each training set size, the recognition rate is obtained by taking the mean recognition rate of the 10 tests. On this database the LASSO method has enabled to select quite the same number of variables than with the initial database: 12 variables on average from GFD features, 13 from Zernike features and 13 from R-signature 1D features. Let us consider Table 4. The used notations are the same as the ones previously used in Table 3. Moreover, the notation DF means that the three discrete shape measures have been used. The recognition rates show the descriptor combination interest. Indeed, even if the classification is less efficient on this database than on the initial one, the results show the combination of continuous descriptors improves the recognition rate. Moreover, the addition of the 3 discrete shape measures outperforms these results. In fact, the integration of the discrete measures improves the recognition rate by 3.8% on average compared to the recognition rate obtained by the combination of the 3 continuous descriptors. Finally, Table 5 shows that the proposed GM-B classifier performs always better than the SVM and the FKNN classifiers. In the same way, Table 6 shows the maximal and minimal values and the mean and the standard deviation of the recognition rates obtained by the 3 compared classifiers, during the 10 tests for a training on 50% of the database. The standard deviation is small whatever the classifier and shows a low variability of recognition rate following the different training and testing sets.

Table 7 shows CPU times of the SVM classifier, the FKNN, and the proposed GM-B classifier, for training and for the SVM and the GM-B classifier), the SVM classifier is faster than the two others. The CPU time is higher for the GM-B model because it depends on the number of Gaussians (in this case, 2) and the pre-defined precision of the EM algorithm. However, the processing time remains weak since it takes less than 0, 03 s per image. Following the discussion in Sect. 4, we can remark that without variable selection the CPU time raises drastically for the GM-B classifier.

Conclusion and future works

In this paper, we have proposed an original adaptation of the Bayesian theory to combine descriptors. We have shown that a Bayesian network has good properties for symbol recognition. In the proposed model, the bad behavior of a descriptor does not impede the behavior of the others. Moreover, we can take into account different types of descriptors. Indeed, we have combined discrete shape measures with continuous shape descriptors. This combination provides a classifier more robust to variability and scalability. Moreover, we have adapted the LASSO method, which solves our dimensionality problem and thus decreases our method complexity and especially increases the recognition rate. The experimental results are very promising and show the efficiency of our method.

In our future works, we want to use our approach in the case of very complex symbols like electrical wiring diagrams. In this case, to recover the maximum amount of information, it is useful to add more descriptors in our combination framework. Even if the LASSO has shown good results, it does not take into account correlation between variables and after selection, variables are always correlated. However, it is well

Fig. 1

 1 Fig. 1 Naïve Bayes

Fig. 3

 3 Fig.3The Gaussian-mixtures and Bernoulli mixture model

Fig. 4

 4 Fig. 4 GREC symbol database

Fig. 5

 5 Fig. 5 Examples of noisy and occluded symbols for different models

Table 1

 1 Mean numbers of variables in function of variable selection method

	Number of variables from	GFD	Zernike	R-signature 1D
	Without selection	225	34	180
	SFS	141	34	48
	LASSO	13	15	13

Table 2

 2 8% on average compared to the SFS selection. Thus the LASSO method is more robust, experimentally, on this database, than the SFS method. In fact, the shrinkage methods like the LASSO are well known to be more stable than iterative methods, to select variables among a large set of variables but with few examples. Thus, the variables selected with the LASSO method have been used for our following experiments.

	Mean recognition rates (in %) for SVM classifier,	Variable selection method	SVM classifier	FKNN k = 1	FKNN k = m	GM-B
	FKNN and GM-B in function of variable selection method	Without selection Random selection		87.6 90.8	8 9 .9 9 3 .7	8 8 .6 9 1 .5	8 9 .8 9 3 .3
					SFS			94.1	9 7 .2	9 5 .3	9 6 .7
					LASSO			95.7	9 8 .8	9 6 .2	100
	Table 3 Mean recognition rates (in %) of GM-B after variable selection
	with the LASSO
	G	Z	R	G + Z	G + R	Z + R	G + Z+ R
	99	100	46.1	100	99.3	100	100
	the LASSO method improves the recognition rate by 8.7%
	on average compared to the classification without variable
	selection, by 5.3% on average with random selection, and
	by 1.

Table 4

 4

	Mean recognition rates (in %) of GM-B after variable	Training part (%)	G	Z	R	G + Z	G + R	Z + R	G + Z + R	G + Z + R + DF
	selection with the LASSO-database including occluded symbols	25 50	70.4 71	79 80.7	39.3 40.2	85.5 87.6	75.5 76.3	82.2 83.4	93.3 93.7	96.8 98.6
		75	75.7	85.1	41.2	89.4	79.1	87.6	96.2	99.2

Table 5

 5 Mean recognition rates (in %), by combining continous and discrete features (G + Z + R + DF), for SVM classifier, FKNN and GM-B after variable selection with the LASSO-database including occluded symbols

	Training part (%)	SVM classifier	FKNN k = 1	FKNN k = m	GM-B
	25	89.2	91.9	91.7	96.8
	50	91	95.2	93	98.6
	75	92.5	97.1	94.7	99.2

Table 6

 6 Statistical measures (in %) for SVM classifier, FKNN and GM-B recognition rates, after variable selection with the LASSO, by combining continous and discrete features (G + Z + R + DF)-database including occluded symbols (training set = 50% of the database)

	Measure	SVM classifier	FKNN k = 1	FKNN k = m	GM-B
	Min	90.4	9 4 .8	9 2 .9	9 8 .5
	Max	91.7	9 5 .4	9 3 .03	98.65
	Mean	91	95.2	9 3	9 8 .6
	Standard deviation	0.4	0 .17	0.04	0.07

The CPU times are given for all the test images known that the descriptors are often partially redundant since they address the same task. In this case, a more appropriate reduction method should be investigated. Moreover, it can be interesting to annotate some symbols and add the information given by possible keywords associated with a subset of training data. Our motivation is based on the property of Bayesian networks to enable to manage, in a same network, different kinds of information (in this case different media), and to their ability to handle missing values.