Stanislav Harizanov

Jean-Christophe Pesquet

Gabriele Steidl

Epigraphical Projection for Solving Least Squares Anscombe Transformed Constrained Optimization Problems

This papers deals with the restoration of images corrupted by a non-invertible or ill-conditioned linear transform and Poisson noise. Poisson data typically occur in imaging processes where the images are obtained by counting particles, e.g., photons, that hit the image support. By using the Anscombe transform, the Poisson noise can be approximated by an additive Gaussian noise with zero mean and unit variance. Then, the least squares difference between the Anscombe transformed corrupted image and the original image can be estimated by the number of observations. We use this information by considering an Anscombe transformed constrained model to restore the image. The advantage with respect to corresponding penalized approaches lies in the existence of a simple model for parameter estimation. We solve the constrained minimization problem by applying a primal-dual algorithm together with a projection onto the epigraph of a convex function related to the Anscombe transform. We show that this epigraphical projection can be efficiently computed by Newton's methods with an appropriate initialization. Numerical examples demonstrate the good performance of our approach, in particular, its close behaviour with respect to the I-divergence constrained model.

Introduction

The Poisson distribution exhibits a mean/variance relationship. This mean/variance dependence can be reduced by using variance-stabilizing transformations (VST), one of which is the Anscombe transform [START_REF] Anscombe | The transformation of Poisson, binomial and negative-binomial data[END_REF] defined as

T : [0, +∞) n → (0, +∞) n : v = (v i ) 1≤i≤n → 2 v i + 3 8 1≤i≤n .
It transforms Poisson noise to approximately Gaussian noise with zero-mean and unit variance. The Anscombe transform has been employed in order to solve inverse problems where one wants to recover an original signal u ∈ [0, +∞) m from observations f = P(Hu),

where P denotes an independent Poisson noise corruption process and H ∈ [0, +∞) n×m is a linear degradation operator, e.g. a blur. Note that we consider images of size M × N columnwise reshaped as vectors of length m = M N .

In this context, one of the possible uses of the Anscombe transform is i) to transform the degraded observations f , ii) to apply a data recovery technique which is valid for an additive white zero-mean Gaussian model and iii) to apply an inverse transform to the so-recovered signal [START_REF] Chaux | Wavelet-based restoration methods: Application in 3d confocal microscopy images[END_REF] (see also [START_REF] Mikkitalo | Optimal inversion of the Anscombe transformation in low-count Poisson image denoising[END_REF] for more recent developments). Note that this method appears mainly to be well-founded for denoising problems. When a linear degradation operator H is present, a better approach consists of adopting a variational framework [START_REF] Chaux | Nested iterative algorithms for convex constrained image recovery problems[END_REF][START_REF] Dupé | A proximal iteration for deconvolving Poisson noisy images using sparse representations[END_REF] where one minimizes a data fidelity term

u → T (Hu) -T (f ) 2 2 (1)
penalized by a (sum of) regularization term(s) serving to incorporate prior information about the sought signal u. The approach is also closely related to a Maximum A Posteriori (MAP) estimate, where the function in [START_REF] Anscombe | The transformation of Poisson, binomial and negative-binomial data[END_REF] is substituted for the neg-log-likelihood of the Poisson noise, i.e., the I-divergence (generalized Kullback-Leibler divergence)

u → D(f, Hu) := 1 n , f log f Hu -f + Hu if Hu > 0, +∞ otherwise,
where • , • denotes the standard Euclidean inner product and 1 n denotes the vector consisting of n entries equal to 1 (see [START_REF] Figueiredo | Restoration of Poissonian images using alternating direction optimization[END_REF][START_REF] Setzer | Deblurring Poissonian images by split Bregman techniques[END_REF]). One of the drawbacks of these penalized methods is that multiplicative constants weighting the regularization terms (the so-called regularization parameters) need to be set carefully, which may be a difficult task.

A way of circumventing this problem consists of adopting a constrained approach instead of a regularized one, by imposing that

T (Hu) -T (f ) 2 2 ≤ τ (2)
where τ ∈ [0, +∞). Based on the statistical properties of the Anscombe transform and the law of large numbers, a consistent choice for the above bound is τ = n, when the number of observations n is large. In this work, we will investigate such an approach by solving the following problem:

minimize u∈C Φ(Lu) subject to T (Hu) -T (f ) 2 2 ≤ τ, (3)
where C is a nonempty closed convex subset of [0, +∞) m , L ∈ R q×m , and Φ : R q → (-∞, +∞] is a proper, lower-semicontinuous, convex function. A typical choice for C is the nonnegative orthant of R m . The classical Total Variation objective function [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] is obtained, as a special case, when Φ is an ℓ 2,1 norm and L corresponds to a discrete gradient operator. Constrained models based on the I-divergence have been considered in [START_REF] Carlavan | Sparse Poisson noisy image deblurring[END_REF][START_REF] Teuber | Minimization and parameter estimation for seminorm regularization models with I-divergence constraints[END_REF], where in the second paper special attention was paid to the relation between the parameters of the constrained and the penalized problem via discrepancy principles. Note that recently penalized versus constrained problems in a rather general form were handled in [START_REF] Aravkin | Variational properties of value functions[END_REF]. In [START_REF] Cherchia | A proximal approach for constrained cosparse modelling[END_REF], the I-divergence constraint was replaced through a polyhedral approximation technique and an epigraphical projection method was applied to solve the problem. In this work, we will also take advantage of an epigraphical projection approach to solve the Anscombe constrained model (3) and we will show that the required epigraphical projections can be easily determined in this context. The structure of this paper is as follows: Section 2 recalls the notation. In Section 3 we determine the epigraphical projection for a function related to constraint ( 2) which plays a central role in the primal dual algorithms established in Section 4. In particular, we provide a good starting point for the involved Newton method. Numerical examples are presented in Section 5 emphasizing the good approximation of the I-divergence constrained approach achieved by our Anscombe constrained model. Finally, a summary of our contribution and some conclusions are given in Section 6.

Notation

Let Γ 0 (R n ) denote the set of proper, lower-semicontinuous, convex functions mapping from

R n to (-∞, +∞]. The epigraph of ϕ ∈ Γ 0 (R n ) is the nonempty, closed, convex subset of R n+1 defined as epi ϕ := {(v, ζ) ∈ R n × R : ϕ(v) ≤ ζ}.
For a nonempty, closed, convex set C ⊂ R m we denote by ι C ∈ Γ 0 (R m ) its indicator function

ι C (u) := 0 if u ∈ C, +∞ otherwise,
and by P C the orthogonal projector onto C. Beyond epigraphs of functions from Γ 0 (R n ) we will consider the half-space V τ := {ζ ∈ R n : 1 n , ζ ≤ τ }. Using this notation and defining, for every i ∈ {1, . . . , n},

ϕ i : [0, +∞) → [0, +∞) : s → 2 √ s -T (f ) i 2 , problem (3)
can be rewritten as

minimize (u,ζ)∈R m ×R n ι C (u) + Φ(Lu) + n i=1 ι epi ϕi (Hu) i + 3 8 , ζ i + ι Vτ (ζ). (4)
Now one can choose a primal-dual splitting algorithm as those proposed in [START_REF] Bot | Convergence analysis for a primal-dual monotone + skew splitting algorithm with application to total variation minimization[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Vu | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF] to solve this problem. One step in all these algorithms consists of the orthogonal projections onto the epigraphs of ϕ i for all i ∈ {1, . . . , n} which is the topic of the next section.

In this section, we deal with the projection onto the epigraph of the function ϕ ∈ Γ 0 (R) defined as

ϕ(s) := (2 √ s -z) 2 if s ≥ 0, +∞ otherwise, (5)
where z > 0, see Fig. 1. Proposition 1. Let ϕ be defined by (5) with z > 0. Then the epigraphical projection of (x, ζ) ∈ R 2 is given by

z 2 4 z 2 s epi ϕ (x, ζ) Pepi ϕ(x, ζ)
P epi ϕ (x, ζ) =              (max{x, 0}, ζ) if (max{x, 0}, ζ) ∈ epi ϕ, t++z 2
2 , t 2 + if 4x ≥ z 2 , t-+z 2 2 , t 2 - if 4x < z 2 ,
where t + , resp. t -is the unique root in [0, +∞), resp. in (-z, 0) of the cubic polynomial

p : t → 17t 3 + 3zt 2 + (3z 2 -16ζ -4x)t + z(z 2 -4x). ( 6
)
Proof. The function ϕ fulfills ϕ(0) = z 2 and

ϕ ′ (s) = 4 - 2z √ s      < 0 if 0 < s < z 2 4 , = 0 if s = z 2 4 , > 0 if s > z 2 and therefore lim s→0 s>0 ϕ ′ (s) = -∞. Thus, if x ≤ 0 and ζ ≥ z 2 , then P epi ϕ (x, ζ) = (0, ζ). In addition, if (x, ζ) ∈ epi ϕ, then P epi ϕ (x, ζ) = (x, ζ).
We consider the remaining cases when (max{x, 0}, ζ) ∈ epi ϕ. The tangent vector of the curve associated with the graph of ϕ reads (1, ϕ ′ (s)), s > 0. The uniquely determined orthogonal projection (x, ζ) := P epi ϕ (x, ζ) has to satisfy

x ζ -x ζ ⊥ 1 ϕ ′ (x) and ζ = ϕ(x), x > 0 (7)
which leads to

0 = (x -x) √ x + 2 ζ -(2 √ x -z) 2 (2 √ x -z), x > 0.
Substituting t := 2 √ x -z > -z, this can be rewritten as

0 = 17 t3 + 3z t2 + (3z 2 -16ζ -4x) t + z(z 2 -4x) = p( t), t > -z.
Conversely, if t > -z is a root of the polynomial p in ( 6), then x = t+z 2 2 fulfills (7). When x ≥ z 2 /4 then also x ≥ z 2 /4 (see Fig. 1), thus we are interested in the restriction of ϕ to [z 2 /4, +∞) i.e., the nonnegative roots of p. The restriction of ϕ to [z 2 /4, +∞) is convex, monotonically increasing, and (x, ζ) / ∈ epi ϕ. Hence, there is a unique point (x, ζ) on its graph that satisfies [START_REF] Chaux | Wavelet-based restoration methods: Application in 3d confocal microscopy images[END_REF], i.e., p has a unique root in [0, +∞). Analogously, when x < z 2 /4 then also x < z 2 /4, thus we are interested in the restriction of ϕ to [0, z 2 /4) i.e., the roots of p in the interval (-z, 0). The restriction of ϕ to [0, z 2 /4) is convex and monotonically decreasing, and the uniqueness of the root follows by the same arguments. Finally, it can be noticed that ζ = ϕ t+z 2 2 = t2 since t > -z, which completes the proof.

The next proposition states that the root t + , resp. t -, of polynomial p can be computed efficiently by Newton's method with initial value t 0 := 2 max{x, 0}z. Indeed, we have seen in our numerical examples that t 0 is a very good starting point.

Proposition 2. Let (max{x, 0}, ζ) ∈ epi ϕ and t 0 := 2 max{x, 0} -z. Let the polynomial p be defined by [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. Then the Newton method for finding a zero of p with initial value t 0 converges (after a finite number of steps) monotonically to the root

t + if 4x ≥ z 2 , resp., t -if 4x < z 2 . Proof. 1. First we show that i) p(t 0 )p(0) ≤ 0, ii) p ′ (t 0 ) > 0,
where equality in i) holds true iff 4x = z 2 . If x < 0, then t 0 = -z and consequently, since (0, ζ) ∈ epi ϕ, i.e., z 2 > ζ, we obtain:

p(0) = z(z 2 -4x) > 0 and p(t 0 ) = p(-z) = 16z(ζ -z 2 ) < 0, which proves i). Further, since p ′ : t → 51t 2 + 6zt + 3z 2 -4x -16ζ, we obtain p ′ (t 0 ) = p ′ (-z) = 48z 2 -16ζ -4x ≥ 16(4z 2 -ζ) > 0. If x ≥ 0, then t 0 = 2 √ x -z. Consequently, we have p(0) = -z(z + 2 √ x)t 0 and p(t 0 ) = 17t 3 0 + 2zt 2 0 + (3z 2 -16ζ -4x)t 0 + (t 2 0 + z 2 -4x)z = 16t 0 (t 2 0 -ζ).
Since (x, ζ) ∈ epi ϕ, i.e., t 2 0 > ζ we conclude that i) holds true. Finally ii) follows by

p ′ (t 0 ) = 51t 2 0 +6zt 0 +3z 2 -16ζ -4x = 8(6t 2 0 -2ζ +x) = 16(t 2 0 -ζ)+8(4t 2 0 +x) > 0.
Since, in both cases, p(t 0 ) = 0, equality arises in i) iff p(0) = 0 i.e. 4x = z 2 .

2. The following result is well-known, see, e.g., [START_REF] Hanke-Bourgeois | Grundlagen der Numerischen Mathematik und des Wissenschaftlichen Rechnens[END_REF]Theorem 18.3]: Newton's method for finding the unique root of a differentiable, convex, strictly increasing function on an interval converges monotonically if we start at the right endpoint of the interval. There is an analogue result for concave functions.

3. Since p ′′ : t → 6(17t+z), p is convex on [-z 17 , +∞) and concave on (-∞, -z 17 ]. 3.1 Let t 0 > 0 which implies 4x -z 2 > 0. Then p(0) < 0 and, according to Part 1i), p(t 0 ) > 0. Hence, by Proposition 1, t + is the unique root of p in (0, t 0 ). Since p is continuous, p(t + ) = 0 and p(t 0 ) > 0, we necessarily have p ′ (t + ) ≥ 0 (otherwise there would exist another root of p on (t + , t 0 )). Thus, since p is strictly convex, it is strictly monotone increasing on [t + , t 0 ] and we can invoke the argument in Part 2 of the proof. 3.2 Let t 0 < 0 which implies 4x -z 2 < 0. Then, p(0) > 0 and p(t 0 ) < 0 and by Proposition 1 we know that t -∈ (t 0 , 0) is the unique root of p in (-z, 0). If t -≤ -z 17 , then we are done by similar arguments as in 3.1 for concave functions. It remains to study the case when t -> -z 17 . If t 0 > -z 17 , then we know by the strict convexity of p on [t 0 , +∞) that p ′ is strictly increasing on this interval and by Part 1ii) we further have p ′ (t) > p ′ (t 0 ) > 0 for every t > t 0 . Thus, p itself is strictly increasing on [t 0 , +∞). Consequently, one Newton step with initialization t 0 generates t 1 = t 0 -p(t 0 )/p ′ (t 0 ) and the convexity inequality p(t 1 ) ≥ p(t 0 ) + p ′ (t 0 )(t 1 -t 0 ) shows that p(t 1 ) ≥ 0. We thus are in the setting of Part 2 and the method converges monotonically. Finally, let t 0 ≤ -z 17 so that 4x ≤ 16 17

2 z 2 . (8)
Since we must have p(-z 17 ) < 0, a straightforward calculation yields

p - z 17 = 2z 17
z 2 17 + 7z 2 + 8ζ -32x < 0 ⇔ 8ζ < 32x - z 2 17 -7z 2 . ( 9
)
We shall now prove that p ′ (t) > 0 for all t ∈ R, so that p is strictly increasing. Since p ′ is a quadratic polynomial with minimimum at -z 17 , we have only to show that p ′ -z 17 = 2 24 17 z 2 -8ζ -2x > 0. Plugging in ( 8) and ( 9), we indeed obtain

1 2 p ′ - z 17 > 24 17 z 2 -34x + z 2 17 + 7z 2 > 25 17 z 2 - 128 17 z 2 + 7z 2 = 16 17 z 2 > 0.
The sequence (t k ) k∈N generated by Newton's algorithm is such that there exists k 0 ∈ N \ {0} such that t k0 ≥ t -. Otherwise, (t k ) k∈N would be an increasing sequence which would necessarily converge to t -and there would exist k 1 ∈ N such that p is convex over [t k1 , +∞[. Then, we would have p(t k1+1 ) ≥ 0. Thus, after a finite number of steps, the algorithm arrives at t k0 ≥ t -and we can apply Part 2 of the proof. 3.3 Let t 0 = 0 which implies 4z -t 2 = 0. Then, we get t + = t 0 .

Primal-Dual Algorithms

We can apply the projection onto the epigraph of ϕ in combination with any primal-dual algorithm proposed in [START_REF] Bot | Convergence analysis for a primal-dual monotone + skew splitting algorithm with application to total variation minimization[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF][START_REF] Combettes | Primal-dual splitting algorithm for solving inclusions with mixtures of composite, Lipschitzian, and parallel-sum type monotone operators[END_REF][START_REF] Vu | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF] or an alternating direction method of multipliers. For example, we use here the primal-dual hybrid gradient algorithm from [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Pock | A convex relaxation approach for computing minimal partitions[END_REF] with an extrapolation (modification) of the dual variable which will be designated by PDHGMp. Based on the following reformulation of (4), minimize

(u,ζ),(v1,v2,η) ι C (u) + ι Vτ (ζ) + Φ(v 2 ) + n i=1 ι epi ϕi (v 1,i , η i ) subject to   H 0 L 0 0 I   u ζ +   3/8 0 0   =   v 1 v 2 η   , (10)
this algorithm reads:

Algorithm 1 (PDHGMp for solving the Anscombe constrained problem) Initialization: u (0) , ζ (0) , (p

(0) j ) 1≤j≤3 = (p (0) j ) 1≤j≤3 , θ ∈ (0, 1], (ρ, σ) ∈ (0, +∞) 2 with ρσ < 1/ max{1, H * H + L * L 2 } For k = 0, 1, . . . repeat until a stopping criterion is reached 1. u (k+1) = PC u (k) -σρ H * p(k) 1 + L * p(k) 2 2. ζ (k+1) = PV τ ζ (k) -σρp (k) 3 3. (v (k+1) 1,i , η (k+1) i ) = Pepi ϕ i p (k) 1,i + (Hu (k+1) )i + 3/8 , p (k) 3,i + ζ (k+1) i , i = 1, . . . , n 4. v (k+1) 2 = prox σ -1 Φ (p (k) 2 + Lu (k+1) ) 5. p (k+1) 1 = p (k) 1 + Hu (k+1) + 3/8 -v (k+1) 1 6. p (k+1) 2 = p (k) 2 + Lu (k+1) -v (k+1) 2 7. p (k+1) 3 = p (k) 3 + ζ (k+1) -η (k+1) 8. p(k+1) j = p (k+1) j + θ(p (k+1) j -p (k) j ), j = 1, 2, 3.
The projection in step 1 is quite simple if C is the nonnegative orthant of R m .

Step 3 requires the epigraphical projections discussed in the previous section,

Step 4 can be performed by coupled soft shrinkage with threshold σ -1 if we use the ℓ 2,1 -norm. The other steps can be computed in a straightforward way.

We will compare this algorithm with PDHGMp applied to the I-divergence constrained problem minimize u∈C Φ(Lu) subject to D(f, Hu) ≤ τ I (11) using a similar splitting to [START_REF] Chierchia | Epigraphical projection and proximal tools for solving constrained convex optimization problems -part I[END_REF] but without the extra-variables ζ and η:

Algorithm 2 (PDHGMp for solving the I-divergence constrained problem) Initialization: u (0) , (p

(0) j ) 1≤j≤2 = (p (0) j ) 1≤j≤2 , θ ∈ (0, 1], (ρ, σ) ∈ (0, +∞) 2 with ρσ < 1/ H * H + L * L 2 For k = 0, 1, . . . repeat until a stopping criteria is reached 1. Step 1 of Algorithm 1 2. minimize v 1 v1 -(p (k) 1 + Hu (k+1) ) 2 2
subject to D(f, v1) ≤ τI as in [START_REF] Teuber | Minimization and parameter estimation for seminorm regularization models with I-divergence constraints[END_REF]. p(k+1)

j = p (k+1) j + θ(p (k+1) j -p (k) j ), j = 1, 2.
In [START_REF] Zanella | Efficient gradient projection methods for edge-preserving removal of Poisson noise[END_REF] (see also [START_REF] Bardsley | Regularization parameter selection methods for ill-posed Poisson maximum likelihood estimation[END_REF]) statistical arguments were used to show that τ I = 1 2 n is a good estimate in case of moderate Poisson noise. In [START_REF] Carlavan | Sparse Poisson noisy image deblurring[END_REF] this estimate was improved in case f has many zero components.

Numerical Examples

In this section, we demonstrate the performance of our algorithm by numerical examples implemented in MATLAB (Intel Core i7-870 Processor with 8M Cache, 2.93 GHz, 8 GB physical memory). We have tested the two original images u, namely 'cameraman' (256 × 256) and 'brain' (184 × 140), depicted in Fig. 2 and denoted by B1, resp. B2 in the following. ν is the intensity of the image. We tested ν = 100, 600, 1200, 2000, 3000 and denoted the blurred, noisy images by B1 ν and B2 ν . In order to synthetically add Poisson noise to the noise-free image, we applied the MATLAB routine imnoise(X,'poisson'). For a quantitative comparison of the images, we computed the peak signal to noise ratio (PSNR) and the (MAE) defined by PSNR = 10 log 10

| max u-min u| 2 1 n u-u 2 2
, and MAE = 1 nν u -u 1 . The 'true' constraints between the blurred, noisy image f and the original image u are given in Table 1. As can be seen, the estimates τ I = n/2 and τ = n are good approximations of the true constraints D(f, Hu) and T (Hu) -T (f ) 2 2 . We computed a minimizer of our functional (3) with C = [0, +∞) n , the ℓ 2,1norm for Φ, and the discrete gradient operator L (q = 2n) by using Algorithm 1 with τ = n. We compared the result with the I-divergence constrained approach (11) and Algorithm 2 with τ I = n/2. The parameters σ and ρ appearing in PDHGMp (in this setting convergence is theoretically guaranteed for σρ < 1/9) are fitted such that the algorithms give (up to two digits after the comma) the same PSNR, MAE and TV-norm (times 10 5 or 10 6 ) after 1000 iterations as after 100000 iterations. In the I-divergence constrained approach with the brain data, we stopped after 5000 iterations. Furthermore, we have chosen θ = 1. Fig. 3 shows the restoration results for B1 1200 and Fig. 4 for B2 1200 . Table 2 summarizes the results for the different intensities. As expected, we observe that the outcomes of the two algorithms are very similar. More precisely, if u A , resp. u I denotes the output of the restoration procedure with Anscombe, resp. I-divergence constraints, then we get for image B1 that u A -u I 2 /(ν √ n) ranges from 0.004 to 1.74e-4 and max |u A -u I |/ν from 0.0612 to 0.0031 for the different noise levels. Finally, Table 3 compares Algorithms 1 and 2 for different constraints τ and τ I and a central part of the cameraman of size 130 × 130 with ν = 3000 after 1000 iterations.

Summary and Conclusions

We have considered a constrained restoration model for images corrupted by a linear transform and Poisson noise by making use of the Anscombe transform. 1.

In contrast with penalized approaches, a main advantage of the proposed one is that it makes it possible to employ a simple estimate for the model parameter. We have provided proximal algorithms to find a minimizer of the model, which are based on epigraphical projections, and we have shown that the performance is similar to a recently introduced I-divergence constrained model. Future research directions include the following: i) replacing or combining the discrete gradient operator L with other ones (discrete higher order operators, nonlocal means, wavelet-like transforms) and handling other problems than deblurring ones, ii) considering convex optimization problems involving multiple constraints for which the epigraphical projection approach may be quite efficient, see [START_REF] Chierchia | Epigraphical projection and proximal tools for solving constrained convex optimization problems -part I[END_REF], iii) restoring images with Poisson+Gauss noise, see [START_REF] Jezierska | A primal-dual proximal splitting approach for restoring data corrupted with Poisson-Gaussian noise[END_REF][START_REF] Li | A reweighted ℓ2 method for image restoration with Poisson and mixed Poisson-Gaussian noise[END_REF], and iv) finding numerically efficient methods to map the constraint bound to the parameter of the corresponding penalized functional.

Fig. 1 .

 1 Fig. 1. The epigraph of ϕ for z = 3 and the epigraphical projection Pepi ϕ(x, ζ) of some point (x, ζ).

3. - 5 .

 5 Steps 4. -6. of Algorithm 1 6.

Fig. 2 .

 2 Fig. 2. Original images 'cameraman' (left) and phantom of a brain image (right).

  The images were blurred by a matrix H corresponding to a Gaussian kernel with standard deviation 1.3 and mirrored boundary (we have then m = n). Their gray values are interpreted as photon counts in the range [0, ν], where D(f, Hu)/n T (Hu) -T (f )

Fig. 3 .

 3 Fig. 3. Result for the 'cameraman' image B11200 corresponding to Table 3. Corrupted image (left), restoration result by the Anscombe constrained model and Alg. 1 (middle), difference image between the middle image and the one recovered by I-divergence constrained model and Alg. 2 (right). The gray values in the difference image are between -10 and 10, while the image values were scaled up to 1200.

  ) -T (f )

Fig. 4 .

 4 Fig. 4. Result for the 'brain' image B21200 corresponding to Table 3. Corrupted image (left), restoration result by the Anscombe constrained model and Alg. 1 (middle), difference image between the middle image and the one recovered by the I-divergence constrained model and Alg. 2 (right). The gray values in the difference image are between -15 and 15, while the image values were scaled up to 1200.

  ) -T (f )

Table 1 .

 1 The original values of D(f, Hu)/n, T (Hu) -T (f ) 2 2 /n and PSNR, MAE of f .

				2 2 /n PSNR MAE
	B1100	0.5075	1.0086	20.58 66.41e-3
	B1600	0.5020	1.0034	23.38 41.59e-3
	B11200	0.5018	1.0039	23.79 37.17e-3
	B12000	0.4979	0.9960	23.97 34.87e-3
	B13000	0.4994	0.9989	24.06 33.58e-3
	B2100	0.4954	0.8866	18.34 82.45e-3
	B2600	0.5131	1.0004	20.17 61.01e-3
	B21200	0.5122	1.0178	20.37 57.90e-3
	B22000	0.5063	1.0085	20.49 56.02e-3
	B23000	0.4956	0.9899	20.52 55.16e-3

Table 2 .

 2 Results of Algorithms 1 and 2 with τ = n and τI = n 2 .

					2 2 -n D(f, Hu) -n/2 TV-norm PSNR MAE
	B1100	0.09 1.225 0.12	2.2837 -	-2.5643	1.2217e+5 24.28 30.76e-3 1.2646e+5 24.39 30.36e-3
	B1600 0.0599 2.8	-1.5323 -	-0.0049	8.9230e+5 25.58 25.69e-3 8.9402e+5 25.59 25.67e-3
	B11200 0.042	3	9.0207 -	-0.045	1.9190e+6 26.08 24.16e-3 1.9198e+6 26.08 24.15e-3
	B12000 0.027 3.03	0.2634 -	--0.1911	3.3007e+6 26.35 23.30e-3 3.3017e+6 26.36 23.30e-3
	B13000 0.0329 4.001	-0.3559 -	-0.2767	5.1667e+6 26.64 22.50e-3 5.1673e+6 26.64 22.49e-3
	B2100 0.55 0.25	0.9730 -	-0.001	0.9349e+5 19.91 59.78e-3 1.0459e+5 20.39 53.35e-3
	B2600	0.040 3.04 0.050	-2.7434 -	-9.0129	7.4158e+5 21.81 40.05e-3 7.5693e+5 21.91 39.28e-3
	B21200	0.034 3.97 0.042	0.0079 -	-4.2683	1.5598e+6 22.33 36.26e-3 1.5673e+6 22.36 36.11e-3
	B22000	0.021 4.108 0.041	0.7789 -	-0.8585	2.6780e+6 22.72 33.53e-3 2.6885e+6 22.73 33.48e-3
	B23000	0.0182 5.413 0.0282	-0.5284 -	-1.5536	4.0565e+6 22.93 32.08e-3 4.0603e+6 22.94 32.03e-3

Table 3 .

 3 2 2 -n D(f, Hu) -n/2 TV-norm PSNR MAE Results of Algotithms 1 and 2 on a part of B13000 for different constraining parameters τ = scale•n and τI = scale• n 2 . The optimal is computed as in Table

	0.8	0.2	0.975	-0.0813 -	--0.0585	2.2400e+6 26.51 22.03e-3 2.2416e+6 26.52 22.02e-3
	1.0	0.0269	4	-2.4868 -	--0.0425	1.7071e+6 25.60 21.20e-3 1.7073e+6 25.60 21.19e-3
	1.0095 (optimal) 0.0239 4.002	-0.0616 -	--1e-4	1.6964e+6 25.56 21.29e-3 1.6914e+6 25.54 21.33e-3
	1.2	0.011	8	1.1405 -	-0.0611	1.5563e+6 24.93 22.96e-3 1.5565e+6 24.93 22.96e-3
	2	0.004 30.065	-1.7501 -	--0.0257	1.3388e+6 23.70 27.02e-3 1.3391e+6 23.71 27.02e-3


Henri Casanova
email: henric@hawaii.edu

Fanny Dufossé
email: dufosse@i3s.unice.fr

Yves Robert

Frédéric Vivien

Scheduling Tightly-Coupled Applications on Heterogeneous Desktop Grids

Platforms that comprise volatile processors, such as desktop grids, have been traditionally used for executing independent-task applications. In this work we study the scheduling of tightly-coupled iterative masterworker applications onto volatile processors. The main challenge is that workers must be simultaneously available for the application to make progress. We consider three additional complications: (i) one should take into account that workers can become temporarily reclaimed and, for data-intensive applications; (ii) one should account for the limited bandwidth between the master and the workers; and (iii) workers are strongly heterogeneous, with different computing speeds and availability probability distributions.

In this context, our first contribution is a theoretical study of the scheduling problem in its off-line version, i.e., when processor availability is known in advance. Even in this case the problem is NP-hard. Our second contribution is an analytical approximation of the expectation of the time needed by a set of workers to complete a set of tasks and of the probability of success of this computation. This approximation relies on a Markovian assumption for the temporal availability of processors. Our third contribution is a set of heuristics, some of which use the above approximation to favor reliable processors in a sensible manner. We evaluate these heuristics in simulation. We identify some heuristics that significantly outperform their competitors and derive heuristic design guidelines.

I. INTRODUCTION

In this paper we study the problem of scheduling parallel applications onto volatile processors. We target typical scientific iterative applications in which a master process parallelizes the execution of each iteration across worker processes. Each iteration requires the execution of a fixed number of tasks, with a global synchronization at the end of each iteration. In [START_REF] Casanova | Scheduling parallel iterative applications on volatile resources[END_REF] we have studied this problem when these tasks are independent. In this work instead we consider tightly-coupled tasks that exchange data throughout each iteration, thus requiring that workers be simultaneously available. This work and that in [START_REF] Casanova | Scheduling parallel iterative applications on volatile resources[END_REF] cover the two extremes of the parallelization spectrum, and are together representative of An extended abstract of this paper appears in Euro-PDP'2013.This work has been supported in part by the French ANR Rescue project. a large class of scientific applications. However, the new requirement that workers be all simultaneously available for tightly-coupled applications introduces another level of difficulty in all resource selection and scheduling decisions.

We consider a platform that consists of processors that alternate between periods of availability and periods of unavailability. When available each processor runs a worker process, and a master process can choose to enroll a subset of these workers to participate in the application execution. Worker unavailability can be due to software faults, in which case unavailability may last only the time of a reboot. A hardware failure can lead to a longer unavailability period, until a repair is completed and followed by a reboot. We consider a third source of processor unavailability, which comes from cyclestealing scenarios: when a processor is contributed to the platform by an individual owner, this owner can reclaim it at any time without notice for some unknown length of time. A difference here is that the processor is merely preempted (as opposed to being terminated) until the processor is no longer reclaimed. A worker process on this processor can later resume its computation. Accordingly, we use a 3-state availability model: UP (available), DOWN (crashed, computation is lost) and RECLAIMED (preempted, but computation can resume later). Our platform model also accounts for the fact that, due to bandwidth limitation, the master is only able to communicate simultaneously with a limited number of workers (to send them the application program as well as task data). This limitation corresponds to the bounded multi-port model [START_REF] Hong | Adaptive allocation of independent tasks to maximize throughput[END_REF]. It turns out that limiting the communication capacity of the master dramatically complicates the design of scheduling strategies. But without this limitation, it would be in principle possible to enroll thousands of new processors at each iteration, which is simply not feasible in practice even if this many processors are available.

Given the above application and platform models, and given a deadline (typically expressed in hours or days), the scheduling problem under study is that of maxi-mizing the expected number of application iterations successfully completed before the deadline. Informally, during each iteration, one must use the "best" processors among those that are simultaneously UP ; these could be the fastest ones, or those expected to remain UP for the longest time. In addition, with processors failing, becoming reclaimed, and becoming UP again later, one has to decide when and how to change the set of currently enrolled processors. Each such change comes at a price: first, the application program needs to be sent to newly enrolled processors, thereby consuming some of the master's bandwidth; second, and more importantly, iteration computation that was only partially completed is lost due to the tight coupling of tasks.

Our contribution in this work is threefold. First, we determine the complexity of the off-line scheduling problem, i.e., when processor availability is known in advance. Even with such knowledge the problem is NP-hard. Second, we compute approximations of the expectation of the time needed by a set of processors to complete a set of tasks and of the probability that this computation succeeds. These approximations provide a sound basis for making sensible scheduling decisions. Third, we design several on-line heuristics that we evaluate in simulation. Some of these contributions assume Markovian processor availability, which is not representative of real-world platforms but provides a tractable framework for obtaining theoretical and experimental results in laboratory conditions. This paper is organized as follows. In Section II, we discuss related work. We give a formal definition of the application and platform models in Section III. In Section IV we define the scheduling problem and establish off-line complexity results. In Section V, we introduce a 3-state Markovian model of processor availability, and use this model to compute approximations of relevant probabilistic quantities. In Section VI, we describe several heuristics for solving the on-line scheduling problem, which are evaluated in simulation in Section VII. Finally, we summarize our results and provide perspectives on future work in Section VII-B.

II. RELATED WORK

Iterative applications that can be implemented in master-worker fashion are widely used in computational linear algebra for sparse linear systems [START_REF] Chazan | Chaotic relaxation[END_REF], [START_REF] Strikwerda | A probabilistic analysis of asynchronous iteration[END_REF], [START_REF] Uc | Partitioning sparse matrices for parallel preconditioned iterative methods[END_REF] or eigenvalue problems [START_REF] Sleijpen | A jacobidavidson iteration method for linear eigenvalue problems[END_REF]), image processing [START_REF] Zhang | A fast parallel algorithm for thinning digital patterns[END_REF], [START_REF] Censor | Component averaging: An efficient iterative parallel algorithm for large and sparse unstructured problems[END_REF], signal processing [START_REF] Byrne | A unified treatment of some iterative algorithms in signal processing and image reconstruction[END_REF], [START_REF] Crouse | Waveletbased statistical signal processing using hidden markov models[END_REF], [START_REF] Katsaggelos | A regularized iterative image restoration algorithm[END_REF], etc. While iterative applications can be asynchronous [START_REF] Strikwerda | A probabilistic analysis of asynchronous iteration[END_REF], [START_REF] Charr | Jacep2pv2: A fully decentralized and fault tolerant environment for executing parallel iterative asynchronous applications on volatile distributed architectures[END_REF], [START_REF] Bahi | Java and asynchronous iterative applications: large scale experiments[END_REF], in this work we only consider the synchronous case.

Several authors have proposed scheduling approaches for such applications [START_REF] Bahi | Parallel Iterative Algorithms: From Sequential to Grid Computing[END_REF], [START_REF] Jagadish | Array architectures for iterative algorithms[END_REF], [START_REF] Heddaya | Mapping parallel iterative algorithms onto workstation networks[END_REF], [START_REF] Legrand | Mapping and load-balancing iterative computations on heterogeneous clusters with shared links[END_REF]. In this work we consider volatile compute resources, such as those found in desktop grids, whose volatility has been studied in [START_REF] Nurmi | Modeling Machine Availability in Enterprise and Wide-area Distributed Computing Environments[END_REF], [START_REF] Wolski | An Analysis of Availability Distributions in Condor[END_REF], [START_REF] Javadi | Mining for Statistical Models of Availability in Large-Scale Distributed Systems: An Empirical Study of SETI@home[END_REF], [START_REF] Bouteiller | MPICH-V2: a Fault Tolerant MPI for Volatile Nodes based on Pessimistic Sender Based Message Logging[END_REF]. Several authors have studied the "bag-of-tasks" scheduling problem on these platforms, either at an Internet-wide scale or within an Enterprise [START_REF] Kondo | Resource Management for Rapid Application Turnaround on Enterprise Desktop Grids[END_REF], [START_REF] Zhou | Wave Scheduler: Scheduling for Faster Turnaround Time in Peer-based Desktop Grid Systems[END_REF], [START_REF] Estrada | The Effectiveness of Threshold-Based Scheduling Policies in BOINC Projects[END_REF], [START_REF] Anglano | Fault-aware scheduling for Bag-of-Tasks applications on Desktop Grids[END_REF], [START_REF] Byun | MJSA: Markov job scheduler based on availability in desktop grid computing environment[END_REF], [START_REF] Estrada | A distributed evolutionary method to design scheduling policies for volunteer computing[END_REF], [START_REF] Wingstrom | Probabilistic Allocation of Tasks on Desktop Grids[END_REF], [START_REF] Heien | Computing Low Latency Batches with Unreliable Workers in Volunteer Computing Environments[END_REF]. Most of these works propose simple greedy scheduling algorithms that rely on mechanisms to select processors based on static criteria (e.g., processor clock-rates, benchmark results, time zone), on simple statistics of past availability [START_REF] Kondo | Resource Management for Rapid Application Turnaround on Enterprise Desktop Grids[END_REF], [START_REF] Estrada | The Effectiveness of Threshold-Based Scheduling Policies in BOINC Projects[END_REF], [START_REF] Estrada | A distributed evolutionary method to design scheduling policies for volunteer computing[END_REF], [START_REF] Zhou | Wave Scheduler: Scheduling for Faster Turnaround Time in Peer-based Desktop Grid Systems[END_REF], and on predictive statistical models of availability [START_REF] Wingstrom | Probabilistic Allocation of Tasks on Desktop Grids[END_REF], [START_REF] Heien | Computing Low Latency Batches with Unreliable Workers in Volunteer Computing Environments[END_REF], [START_REF] Anglano | Fault-aware scheduling for Bag-of-Tasks applications on Desktop Grids[END_REF], [START_REF] Byun | MJSA: Markov job scheduler based on availability in desktop grid computing environment[END_REF]. These criteria are used to rank processors but also to exclude them from consideration [START_REF] Kondo | Resource Management for Rapid Application Turnaround on Enterprise Desktop Grids[END_REF], [START_REF] Estrada | The Effectiveness of Threshold-Based Scheduling Policies in BOINC Projects[END_REF]. The work in [START_REF] Byun | MJSA: Markov job scheduler based on availability in desktop grid computing environment[END_REF] is particularly related to our own in that it uses a Markov model of processor availability (but without accounting for temporary preemption). Most of these works also advocate for task replication as a way to cope with volatile resources. Expectedly, injecting task replicas is sensible toward the end of application execution. Given the wealth of scheduling approaches, in [START_REF] Estrada | A distributed evolutionary method to design scheduling policies for volunteer computing[END_REF] the authors propose to automatically instantiate the parameters that together define the behavior of a generic scheduling algorithm. Most works published in this area are of a pragmatic nature and few theoretical results have been sought or obtained (one exception is the work in [START_REF] Fujimoto | Near-Optimal Dynamic Task Scheduling of Independent Coarse-Grained Tasks onto a Computational Grid[END_REF]). Note that while in this work we focus on scheduling a single application, other authors have studied the scheduling problem for multiple simultaneous applications [START_REF] Anglano | Scheduling algorithms for multiple Bag-of-Task applications on Desktop Grids: A knowledge-free approach[END_REF], [START_REF] Celaya | A Fair Decentralized Scheduler for Bag-of-Tasks Applications on Desktop Grids[END_REF].

A key aspect of our work is that we seek to develop scheduling algorithms that explicitly manage the master's bandwidth. Limited master bandwidth is a known issue for desktop grid computing [START_REF] Moretti | Challenges in Executing Data Intensive Biometric Workloads on a Desktop Grid[END_REF], [START_REF] Toyoma | A Resource Management System for Data-Intensive Applications in Desktop Grid Environments[END_REF], [START_REF] He | BLAST Application with Data-Aware Desktop Grid Middleware[END_REF] and must therefore be addressed even though it complicates the scheduling problem. To the best of our knowledge, except our previous study for independent tasks [START_REF] Casanova | Scheduling parallel iterative applications on volatile resources[END_REF], no previous work has made such an attempt.

III. MODELS AND ASSUMPTIONS

We assume that time is discretized into a sequence of time-slots of arbitrarily chosen duration. For simplicity when we say "at time t" we imply "at discrete time-slot t". Our approach is agnostic to the time-slot duration. The duration that makes sense in practice depends on the application and/or platform, ranging from seconds to minutes or possibly hours.

A. Application model

We consider an application that performs a sequence of iterations. Each iteration consists of executing m tasks and ends with a global synchronization. All m tasks are identical (in terms of computational cost) and communicate throughout the iteration execution. Therefore, all tasks must make progress at the same rate. If a task is terminated prematurely (due to a worker failure), all computation performed so far for the current iteration is lost, and the entire iteration has to be restarted. If a task is suspended (due to a worker becoming temporarily reclaimed), then the entire execution of the iteration is also suspended. Due to the global synchronization, there is no overlap between communication and computations. We thus consider that an iteration proceeds in two phases: a communication phase and a computation phase. Finally, before being able to compute, a worker must acquire the application code once (e.g., binary executable, byte code), of constant size V prog in bytes, and the input data for each task and iteration, of constant size V data in bytes.

B. Platform model

The platform comprises p processors. Since each processor executes a worker process, we use the terms processor and worker interchangeably. Worker P q , q = 1, . . . , p, can be in one of three states (UP , RECLAIMED or DOWN ), and transitions between these states occur for each processor at each time-slot independently of the other processors. More precisely:

• Any UP processor can become DOWN or RECLAIMED.

• Any UP or RECLAIMED processor can become DOWN . It then loses the application program and all the data for its current tasks. If it was computing some of these tasks, these computations are lost. • Any UP processor can become RECLAIMED.

The processor does not lose any state. If it was receiving the application program or data for a task, the communication is temporarily suspended.

If it was computing a task, the computation on all processors is temporarily suspended. We denote by S q the vector that gives the state of P q at each time-slot starting with time-slot 0.

P q computes a task in w q time-slots if it remains UP . If w q = w for each processor P q , then the processors are homogeneous. The master has network bandwidth BW and communicates with a worker with bandwidth bw, meaning that we assume same capacity links from the master to each worker. Here we equate bandwidth with data transfer rate, acknowledging that in practice the data transfer rate is a fraction of the physical bandwidth. Let n prog be the number of workers receiving the program at time t, and let n data be the number of workers receiving the input data of a task at time t. The constraint on the master's bandwidth writes n prog + n data ≤ n com = BW /bw . Indeed, consider a worker on processor P q that is communicating at time t. Either P q is receiving the program, or it is receiving data for a task. In both cases, it does this at data transfer rate bw. Overall, the master can execute only a limited number n com of such communications simultaneously. The time for a worker to receive the program is T prog = V prog /bw, and the time to receive the data is T data = V data /bw. For simplicity we assume that T prog and T data consist of integral numbers of time-slots. We also assume that the master is always UP , which can be enforced by using, for instance, two dedicated servers with a primary backup mechanism.

C. Application execution model

Let config(t) denote the set of workers enrolled by the master, or configuration, at time t. The configuration is determined by an application scheduler, and in this work we propose algorithms to be used by this scheduler. To complete an iteration, enrolled workers must progress concurrently throughout the computations. One worker may be assigned several tasks and execute them concurrently if it has enough memory to do so. Formally, we define for each worker P q a bound µ q on the maximum number of tasks that it can execute concurrently. We assume that p q=1 µ q ≥ m, otherwise the configuration cannot execute the application. The m tasks are mapped onto k ≤ m workers. Each enrolled worker P q is assigned x q tasks, where k q=1 x q = m. To be able to compute their tasks, the k enrolled workers must have received the application program and all necessary data. More precisely:

• Each enrolled worker P q must receive the program, unless it has received it at some previous time and has not be DOWN since then. • In addition, each worker P q must receive x q data messages (one per task) from the master. Suppose that since the begin of the current iteration P q has received x q data messages. At least (x q -x q )T data time-slots are needed for this communication, and likely more since the master can be engaged in at most n com concurrent communications. Overall, the computation can start at a time t only if each of the k enrolled workers is in the UP state, has the program, has the data of all its allocated tasks, and has never been in the DOWN state since receiving these messages. Because tasks must proceed in locked steps, the execution goes at the pace of the slowest worker. Hence the computation of an iteration requires max q (x q w q ) time-slots of concurrent computations (not necessarily consecutive, due to workers possibly being reclaimed). Consider the interval of time between time t 1 and time t 2 = t 1 + max q (x q w q ) + t -1 for some t . For the iteration to be successfully completed by time t 2 , between t 1 and t 2 there must be max q (x q w q ) timeslots for which all enrolled workers are simultaneously UP , and there may be t time-slots during which one or more workers are RECLAIMED.

The scheduler may choose a new configuration at each time t. If at least one worker in config(t) becomes DOWN , the scheduler must select another configuration and restart the iteration from scratch. Even if all workers in config(t) are UP , the scheduler may decide to change the configuration because more desirable (i.e., faster, more reliable) workers have become available. Let P q be a newly enrolled worker at that point, i.e., P q ∈ config(t + 1) \ config(t). P q needs to receive the program unless it already has a copy of it and has not been DOWN since receiving it. In all cases, P q needs to receive task data, i.e., x q messages of V data bytes. This holds true even if P q had been enrolled at time t < t but was un-enrolled since then. In other words, any interrupted communication must be resumed from scratch if the worker became DOWN or was removed from the configuration.

An example iteration execution with m = 5 tasks and p = 5 processors is shown in Figure 1. For this example the processors are heterogeneous with In the schedule, two tasks are assigned to P 2 and P 3 , and one task is assigned to P 4 , for a workload of 4 time slots on P 2 , 6 time slots on P 3 , and 4 time slots on P 4 . As a result, the iteration computation in this configuration requires 6 time-slots with all processors in the UP state, and during each of these time-slots, 1/6 of each task is executed. This configuration is selected at time 1. At this time, P 1 and P 5 are not UP , so they cannot be included in the configuration.

w i = i for 1 ≤ i ≤ 5, n com = 2, T prog = 2,
The communication phase of this iteration is executed between time 1 and time 7. At time 1, the 3 selected processors can receive data, but because of the bandwidth constraint, P 4 remains idle during the first 3 time-steps. P 3 is temporarily reclaimed after it has downloaded the application program. Due to the bandwidth constraint, processors are idle while others download task data (e.g., P 3 is idle between time 4 and time 5). At time 7, all 3 processors have downloaded the application program and the data for all tasks assigned to them. They all begin computing. At time 10, P 2 is temporarily reclaimed and the computation is suspended with half of the computation of each task completed. When P 2 becomes available again, at time 12, P 3 has been reclaimed and the computation cannot resume immediately. P 3 becomes available again at time 13, P 2 and P 4 are UP , and the computation continues. If a processor had become DOWN , say, at time 14, all the computation would have been lost and the communication phase would have been restarted from scratch. At time 16, the processors synchronize and a new iteration can start. At that time, P 1 and P 5 may be included in the configuration.

IV. OFF-LINE COMPLEXITY

The scheduling problem is to maximize the expected number of completed application iterations before time N , where N is a specified deadline. In this section, we assess the complexity of the off-line version of this problem, assuming full knowledge of future worker states. In other words, S q [j] is known for 1 ≤ q ≤ p and 1 ≤ j ≤ N . We show that the simplest off-line and deterministic versions of the problem are NP-hard.

Fixed number of workers: Consider the problem with no communication (T prog = T data = 0), and identical workers with w q = w and µ q = µ = 1. m workers must be enrolled to complete an iteration. The problem reduces to finding w time-slots such that there exist m workers that are simultaneously UP during all these w time-slots. We call this version of the problem OFF-LINE-COUPLED (µ = 1).

Flexible number of workers: Consider the problem with no communications (T prog = T data = 0), and identical processors with w q = w and µ q = µ = +∞ (in fact µ = m is sufficient). The problem is less constrained than OFF-LINE-COUPLED (µ = 1). Either one finds m processors that are simultaneously UP during w time-slots, or one finds m 2 workers that are simultaneously UP during 2w time-slots, or one finds m 3 workers that are simultaneously UP during 3w time-slots, and so on. We call this version of the problem OFF-LINE-COUPLED (µ = +∞). Proof: We prove both instances with a similar reduction, as detailed below.

(i) NP-completeness of OFF-LINE-COUPLED (µ = 1) . The decision problem associated to OFF-LINE-COUPLED (µ = 1) writes: given a value w and p state vectors S q , can we find m processors that are simultaneously UP during at least w time-steps? This problem clearly belongs to NP: the m × w sub-matrix is a certificate of polynomial (and even linear) size.

For the completeness, we use a reduction from ENCD, the Exact Node Cardinality Decision problem [START_REF] Dawande | On bipartite and multipartitie clique problems[END_REF]. Let I 1 be an instance of ENCD: given a bipartite graph G = (V ∪ W, E) and two integers a and b such that 1 ≤ a ≤ |V | and 1 ≤ b ≤ |W |, does there exist a bi-clique with exactly a nodes in V and b nodes in W ? Recall that a bi-clique

C = U 1 ∪ U 2 is a complete induced sub-graph: U 1 ⊂ V , U 2 ⊂ W , and for every u 1 ∈ U 1 , u 2 ∈ U 2 , the edge (u 1 , u 2 ) ∈ E.
We construct the following instance

I 2 of OFF-LINE- COUPLED (µ = 1): we let p = |V | and N = |W |. Resource R i (which corresponds to vertex v i ∈ V ) is UP at time-step j (which corresponds to vertex w j ∈ W ) if and only if (v i , w j ) ∈ E.
Finally we let m = a and w = b. The size of the instance I 2 is linear in the size of the instance I 1 . We show that I 1 has a solution if and only if I 2 does. Suppose first that I 1 has a solution C = U 1 ∪U 2 . We select the corresponding U 1 processors and the same U 2 time-steps. Because we have a clique, each processor is UP at each time-step, hence I 2 has a solution. Suppose now that I 2 has a solution. The corresponding sub-matrix translates into a bi-clique with a nodes in V and b nodes in W , hence a solution to I 2 .

(ii) NP-completeness of OFF-LINE-COUPLED (µ = ∞). We use the same instance I 1 of ENCD as in (i). We construct the following instance I 2 of OFF-LINE-COUPLED (µ = +∞):

we let p = |V | and N = 2|W |+ 1. Resource R i (which corresponds to vertex v i ∈ V ) is UP at time-step j ≤ N (which corresponds to vertex w j ∈ W ) if and only if (v i , w j ) ∈ E.

V. ANALYTICAL APPROXIMATIONS

In this section, we compute the expectation of the time needed by a configuration to compute a given workload conditioned on this computation being successful (i.e., with no worker becoming DOWN ), as well as the probability of success. Intuitively, these quantities seem relevant for developing scheduling heuristics that account for the need for workers to be UP simultaneously, and for workers that can become temporarily RECLAIMED. To compute the above expectation and probability, we introduce a Markov model of processor availability. The availability of processor P q is described by a 3-state recurrent aperiodic Markov chain, defined by 9 probabilities: P (q) i,j , with i, j ∈ {u, r, d}, is the probability for P q to move from state i at time t to state j at time t + 1, which does not depend on t.

We are aware that the Markov (memory-less) assumption for processor availability does not hold in practice. Several authors have observed that the duration of availability intervals in production desktop grids is often far from being exponentially distributed for a 2-state scenario in which processors are either UP or DOWN [START_REF] Nurmi | Modeling Machine Availability in Enterprise and Wide-area Distributed Computing Environments[END_REF], [START_REF] Kondo | Characterizing Resource Availability in Enterprise Desktop Grids[END_REF], [START_REF] Wolski | An Analysis of Availability Distributions in Condor[END_REF], [START_REF] Javadi | Mining for Statistical Models of Availability in Large-Scale Distributed Systems: An Empirical Study of SETI@home[END_REF], [START_REF] Wingstrom | Probabilistic Allocation of Tasks on Desktop Grids[END_REF]. Unfortunately, there is no consensus in the literature on a realistic model, even though some of these studies may suggest semi-Markov models with approximately Weibull or Log-Normal holding times. Deriving a realistic 3-state statistical model of processor availability is thus an open research question that is outside the scope of this work. Instead, we opt for a Markov model because it is simple and lends itself to tractable analysis. This model gives a framework in which to design heuristics that trade off worker speed for reliability. Furthermore, it allows us to evaluate such trade-off approaches in "laboratory conditions."

A. Probability of success and expected duration of a computation

Consider a set S of workers all in the UP state at time 0. This set is assigned a workload that requires W time-slots of simultaneous computation. To complete this workload successfully, all the workers in S must be simultaneously UP during another W -1 time-slots. They can possibly become RECLAIMED (thereby temporarily suspending the execution) but must never become DOWN in between. What is the probability of the workload being completed? And, if it is successfully completed, what is the expectation of the number of time-slots until completion?

Definition 1: Knowing that all processors in a set S are UP at time-slot t 1 , let P (S) + be the conditional probability that they will all be UP simultaneously at a later time-slot, without any of them going to the DOWN state in between. Formally, knowing that

∀P q ∈ S, S q [t 1 ] = u, P (S)
+ is the conditional probability that there exists a time t 2 > t 1 such that ∀P q ∈ S, S q [t 2 ] = u and S q [t] = d for t 1 < t < t 2 .

Definition 2: Let E (S) (W) be the conditional expectation of the number of time-slots required by a set of processors S to complete a workload of size W knowing that all processors in S are UP at the current time-slot t 1 and none will become DOWN before completing this workload. Formally, knowing that S q [t 1 ] = u, and that there exist W -1 time-slots and E (S) (W ) numerically up to an arbitrary precision ε in fully polynomial time.

t 2 < t 3 < • • • < t W , with t 1 < t 2 , S q [t i ] = u for i ∈ [2, W ], and S q [t] = d for t ∈ [t 1 , t W ], E (S) (W ) is the expectation of t W -t 1 + 1 conditioned on success.
Proof: Consider a set S of processors, all available at time slot 0. Consider the probability P (S) + (t) that all these processors are simultaneously UP again for the first time at time t. This means that for all 0 < t < t, there exists at least one processor RECLAIMED at time t . Also, none of the processors in S goes DOWN between 0 and t.

Let

P (q) u t →u
be the probability that a processor P q that was UP at time 0 is UP again at time t, without having been DOWN in between, and let P (S)

u t →u = Pq∈S P (q) u t →u
. For each processor P q , the value P (q) u t →u can be computed by considering its transition matrix raised to the power t, knowing that the initial state is UP . We form the product to compute P + that all the processors in S will be simultaneously UP again at some point, before the first failure of any of them, is the probability that a processor that was UP at time 0 is UP at time t without having been DOWN . For a processor P q ∈ S, let

P (S) + = t>0 P (S) + (t) = t>0 P (S) u t →u -0<t <t P (S) + (t ) × P (S) u t-t → u = t>0 P (
M q = P (q) u,u P (q) u,r P (q)
r,u P (q) r,r

. Then, P q u t →u = (M t q )[0, 0]. We obtain P q u t →u = µ(λ q 1 ) t + ν(λ q 2 ) t with µ, ν ≥ 0, µ + ν = 1 and λ q 1 > λ q 2 eigenvalues of M q . Then, P q u t →u ≤ (λ q 1 ) t . We obtain P (S)

u t →u ≤ Pq∈S λ q 1 t
and t≥T P (S)

u t →u ≤ Pq∈S λ q 1 T × 1 1-Pq ∈S λ q 1 . Let Λ = Pq∈S λ q 1 . We obtain that T > ln(ε(1-Λ)) ln(Λ) implies E u (S) -0<t<T P (S) u t →u
≤ ε. Thus, we can compute in polynomial time an approximation of E u (S) at ε in polynomial time.

Similarly, we obtain A(S)-0<t<T t×P (S)

u t →u ≤ ε as soon as Λ T T 1-Λ + Λ (1-Λ) 2 ≤ ε.
Therefore A(S) can be approximated with precision ε in polynomial time.

B. Probability of success and expected duration of a communication

The previous section gave approximations for the probability of success and conditional expected duration for computations. Unfortunately, similar approximations cannot be obtained for communications due to complexity added by the n com constraint. Instead, we resort to a coarser approximation as explained hereafter. Let S be a set of enrolled workers. For worker P q ∈ S, let n q be the number of time-slots of communication needed to receive the application program and all the data of its allocated tasks. Suppose first that |S| ≤ n com . In this case, the expected communication time on worker P q , E q , can be estimated precisely reusing the result in the previous section: E q = E (Pq) (n q ). We then estimate the expected communication time of the current configuration as E (S) comm = max Pq∈S {E (Pq) (n q )} . In the case |S| ≥ n com , obtaining an estimate close to the actual expected communication time seems out of reach. Instead, we use a coarser estimation:

E (S) comm = max max Pq∈S E (Pq) (n q ) , Pq ∈S nq ncom . Let P (Pq)
N D (t) denote the probability that worker P q that was UP at time t does not become DOWN between time t and time t +t. The probability of success is then estimated as P (S)

comm = Pq∈S P (Pq) N D (E (S) comm ) .

The expression for P (S)

comm does not take into account the time needed after the end of all communications for all workers to be UP simultaneously. The probability of success of an iteration is estimated by multiplying the probability of success of the communications and the probability of success of the computations.

VI. ON-LINE HEURISTICS

We propose heuristics for solving the on-line version of the scheduling problem, i.e., assuming no knowledge of future processor states. Conceptually, we distinguish between two classes of heuristics. Passive heuristics conservatively keep current processors active as long as possible. In other words, the current configuration is changed only when one of the enrolled processors becomes DOWN . In this case, all previously executed work is lost. However, a worker that has not become DOWN but has already received task data, can reuse that data if the scheduler reassigns tasks to it. Proactive heuristics allow for a complete reconfiguration even if no worker fails, possibly aborting ongoing computation if a better configuration is found. This makes it possible for an iteration to never complete. A criterion must thus be derived to decide whether and when such an aggressive reconfiguration is worthwhile. Our proactive heuristics are defined by a pair (criterion, passive heuristic). When a new configuration is computed using the heuristic, it is compared to the current configuration according to the criterion. If the new configuration is better than the current one, then it is launched, leading to new communications and task allocations. Otherwise, the execution continues with the current configuration for an additional time slot.

We also include results for a baseline RANDOM heuristic that allocates tasks to UP processors randomly using a uniform distribution.

A. Passive heuristics

Passive heuristics assign tasks to workers, which must be in the UP state, one by one until m tasks are assigned. Each task is assigned to a worker according to a criterion that defines the heuristic. As described hereafter, we consider four different criteria: probability of success, expected completion time, estimated yield, and estimated apparent yield. We derive the following heuristics:

• IP (Incremental: Probability of success) -This heuristic attempts to find configurations with high probability of success. The next task is assigned to the worker such that the probability of success of all currently assigned tasks (including the new one) is maximized. More precisely, consider the set S of workers with at least one task already assigned.

For each worker P q , either in S or not, we compute the probability P (S) (q) of success of the communication and the computation if the additional task is assigned to P q , using the results of Section V: P (S) (q) = P (S∪{Pq}) (W q ) × P (S∪{Pq}) comm with W q the maximal load in S ∪ {P q } with an additional task on P q . We assign the next task to worker P q0 , with q 0 = ArgMax P (S) (q) . • IE (Incremental: Expected completion time) -This heuristic attempts to find fast configurations, without considering reliability. The next task is assigned to the worker that minimizes the expected execution time of the iteration. More precisely, consider the set S of workers with at least one task already assigned. For each worker P q , either in S or not, we compute the expected communication time

E (S∪{Pq}) comm
and the expected computation time E (S∪{Pq}) (W q ) with an additional task on P q . We obtain the expected duration of the iteration E (S) (q) = E (S∪{Pq}) comm + E (S∪{Pq}) (W q ). We assign the next task to worker P q0 , with q 0 = ArgMin E (S) (q) . • IY (Incremental: Expected yield) -This heuristic assigns the next task to the worker that maximizes the yield of the configuration. The yield is the expected value of the inverse of the execution time of the current iteration, which we estimate as follows. For a given configuration with probability of success P and expected completion time E for an iteration that has already been running for t time slots, the yield is estimated as Y = P E+t . Intuitively, we expect the yield to achieve a tradeoff between reliability (probability of success) and execution speed. Consider the set S of workers with at least one task already assigned. For each processor P q , either in S or not, we compute the expected yield with an additional task on P q : let P (S) (q) be the probability computed for heuristic IP, E (S) (q) be the expected completion time computed for heuristic IE, and t be the time spent since the beginning of the current iteration. We assign the next task to worker P q0 , with q 0 = ArgMax P (S) (q) t+E (S) (q) . • IAY (Incremental: Expected apparent yield) -

The yield takes into account the time already spent in the current iteration. It could be worthwhile to consider only future work, i.e., the remaining time until iteration completion. To this end we define the apparent yield as AY = P E . Using the same notations as for heuristic IY, we assign the next task to processor P q0 with q 0 = ArgMax P (S) (q) E (S) (q) . B. Proactive heuristics Our proactive heuristics are designed as follows. Consider an application executing on a platform using a passive heuristic H and criterion C at some time t. The configuration config(t -1) was selected by H at time t ≤ t -1 because of a configuration change due to a proactive decision, due to a worker becoming DOWN , or due to the beginning of a new iteration. Let config 1 = config(t ) = config(t -1). At time t , the configuration was measured by criterion C with value c . Suppose that by time t no worker in this configuration has failed. Between t and t, some work may have been done: some communications may be in process or completed, and computations may have started. Consequently, the measure of this configuration given by C should be updated to account for the progress between t and t. Let c be the updated value of criterion C for the current configuration. At step t, a new configuration is computed from scratch using heuristic H, as if no task were allocated to any worker. Let config 2 be this new configuration and c 2 its measure by C. If c ≥ c 2 , then the current configuration at time t -1 is kept for another time-slot: config(t) = config 1 . Otherwise, the current configuration is interrupted, and the new configuration is config(t) = config 2 .

For certain criterion choices, a heuristic could diverge and continually change the configuration, even with workers that are reliably UP . To avoid this divergence, proactive criteria have to respect the following constraint: a given configuration that has been running for t + 1 time-slots must be better for the proactive criterion than the same configuration running for t time slots. With this constraint, all possible configurations are ordered by their value for the selected criterion at the beginning of the iteration, and a lower-ranked configuration in this order cannot be chosen to replace the current configuration. As the number of possible configurations is finite, no proactive heuristic can diverge. The four criteria used to define passive heuristics in the previous section meet this constraint. However, AY (Apparent Yield) leads to many (unnecessary) configuration changes before converging, while the other criteria should be stable. Hence, for the proactive criterion C, we only retain P (Probability of success), E (Expected completion time) and Y (Expected yield). Any passive heuristic H can be used as the building block for a proactive heuristic. We thus obtain 3 × 4 proactive heuristics named C-H where C ∈ {P, E, Y} and H ∈ {IP, IE, IY, IAY}, plus the RANDOM heuristic.

VII. EXPERIMENTAL EVALUATION A. Methodology

To evaluate our heuristics we use a discreteevent simulator, which is publicly available at http: //graal.ens-lyon.fr/ ∼ fdufosse/changing platforms.tar.gz. The input to the simulator are the values for all the parameters listed in Section III. The simulator implements temporal processor availability as Markov processes as described in Section V. All our experiments are for p = 20 processors. The Markov model for processor availability is defined as follows. For each processor P q , we pick a random value uniformly distributed between 0.90 and 0.99 for each P (q)

x,x value (for x = u, r, d). We then set P (q)

x,y to 0.5 × (1 -P (q)

x,x ), for x = y. An experiment is defined by the Markov model for each processor and by three parameters: m, the number of tasks per iteration; n com , the constraint on the master's communication bandwidth; and a third parameter, w min , defined as follows. For each processor P q , we pick w q uniformly between w min and 10×w min . T data is set to w min , meaning that the fastest processor has a computation-communication ratio of 1. T prog is set to 5 × w min , meaning that downloading the program takes 5 times as much time as downloading the data for a task. We define our experimental space as (m, n com , w min ) with m ∈ {5, 10}, n com ∈ {5, 10, 20}, and w min ∈ {1, 2, . . . , 10}. For each possible instantiation of (m, n com , w min ), we create 10 random experimental scenarios so as to obtain different instantiations of the various random parameters. Then, for each experimental scenario, we run 10 trials where each trial uses a different random number generator seed to produce multiple realizations of the Markov chain transitions. The total number of generated problem instances is thus 2 × 3 × 10 × 10 × 10 = 6, 000. We emphasize that our goal here is not to instantiate a representative model for a desktop grid and application, but rather to create arbitrary and simple synthetic experimental scenarios that will make it possible to highlight inherent strengths and weaknesses of the proposed heuristics.

In all experiments, rather than fixing N , the deadline in number of time-slots, we instead fix the number of iterations to 10. The quality of an application execution is then measured by the time needed to complete 10 iterations, or makespan. This problem is equivalent to the problem of maximizing the number of iterations before a deadline. It is also simpler to instantiate since it does not require choosing a meaningful deadline, which would depend on application and platform characteristics. For some experimental scenarios and some heuristics, the time needed to complete 10 iterations successfully can be extremely large, making it impossible to obtain results in a reasonable amount of time. Consequently, we limit the makespan to 1,000,000 seconds and declare that a heuristic fails if it reaches this limit, and succeeds otherwise.

The makespans vary widely between problem instances depending on processor availability patterns, which requires that we define relative metrics for a sound comparison of our heuristics. As expected, with more tasks (larger m), the number of failures increases for each heuristic. With m = 5, no heuristic fails for more than 5 out of the 3,000 problems instances. With m = 10, heuristics fail for up to 6% of the 3,000 problem instances. In both cases, IE is the most robust: It fails only for 1 of the 3, 000 m = 5 instances, and for 81 of the 3, 000 m = 10 instances. Importantly, whenever heuristic IE fails, all the other heuristics also fail. For this reason we use IE as a reference. For a heuristic H, we compute the relative difference (expressed in percentage) between the makespan that it achieves for a given experimental scenario (averaged over all 10 trials) and the one achieved by IE for the same scenario, assuming that heuristic H succeeds. The relative difference is defined as: makespan Hmakespan IE min(makespan H , makespan IE ) .

The denominator is always the makespan of the best performing heuristic, so as to allow consistent comparisons. We denote this metric by %diff. We also count the fraction of trials where heuristic H obtains a makespan smaller than or equal to IE (denoted by %wins) and the fraction of trials where heuristic H obtains a makespan that does not exceed that of IE by more than 30% (denoted by %wins30). We also report the standard deviation over all scenarios (column %stdv). Table I shows results for m = 5 tasks, with heuristics sorted by decreasing %diff. The number of failures for all heuristics is shown in the first column of the table and is at most 5 (recall that IS, the reference heuristic, only fails for 1 instance). Consequently, although some heuristics fail on some scenarios, these failures do not have a large impact on our results.

These results show the efficiency of all our heuristic, when compared with the RANDOM Heuristic. RANDOM is on average more than 20 times worse than IE while all other heuristics have a %diff less than 200%. As seen in the table, only 4 heuristics lead to a %diff value lower than that obtained by IE, with 3 of these heuristics more than 10 points lower. These 4 heuristics are all proactive. We conclude that the best proactive heuristics are significantly better than the best passive heuristics. Several observations can be made on the results in the table. A first one is that using the yield as a heuristic or a criterion is better than using the probability of success. In other terms, heuristic C-IY is better than heuristic C-IP, and heuristic Y-H is marginally better than P-H (in this case an inspection of the simulation traces shows that Y-H and P-H lead to mostly identical executions). Everything else being equal, considering the yield is better than considering the probability of success because it accounts for the host computing speeds in addition to their reliability. A second observation is that heuristic C-IAY is better than heuristic C-IY, thus confirming that the "apparent yield" has merit and is a direct improvement of the yield metric. Anecdotally, while E-IY and E-IAY obtain similar results for %wins and %wins30, E-IAY is significantly better than E-IY on average. A third observation is that although Y-IE and P-IE lead to good results, all other proactive heuristics with the same criteria (i.e., Y-H and P-H) rank last, with %diff values reaching 100%. Finally, the key observation is that the best heuristics (the top 5 heuristics, including the reference IE) all account for expected execution time either as a criterion for selecting a new configuration (E-IAY, E-IY) or as a host selection mechanism (Y-IE, P-IE, IE). A seemingly sensible expectation is thus that E-IE would be very efficient. But instead, E-IE leads to poor results, with on average makespan almost more than 50% longer than that of IE, the fourth lowest %wins value and the fifth lowest %wins30 value. The reason for these poor aggregate results is that E-IE leads to inefficient schedules for problem instances in which the fastest processor is unreliable.

B. Results for m = 10

In this section we discuss results for m = 10 tasks, but only for the reference IE and those 7 heuristics that achieve a %diff value below 50% (the largest such value being in fact below 25%): Y-IE, P-IE, E-IAY, E-IY, IAY, E-IP and IY. Results are shown in Table II. Only two of these heuristics do not consider expected completion time as a criterion: IAY and IY. These two heuristics rank reasonably high in terms of %diff with m = 5 tasks, but are over 130% with m = 10 tasks. The ranking of the heuristics is almost unchanged when compared to the m = 5 results, even if %diff values are lower. When for m = 5, E-IY leads to a negative %diff value, for m = 10 this value becomes positive. For m = 10 only three heuristics achieve positive %diff values: Y-IE, P-IE and E-IAY. With m = 10, most heuristics fail for more than 5% of the problem instances. Given that IE is the most robust heuristic, it should come to no surprise that those proactive heuristics that use IE lead to the lowest number of failures. One conclusion from these results is that Y-IE is only slightly less robust than IE (failing on 4.7% of instances as opposed to 2.7% for IE) but leads to significantly better performance with a %diff value above 11%, leading to a lower makespan for more than 72% of the instances, and leading to a makespan more than 30% larger in less than 8% of the instances.

Figure 2 shows %diff values versus w min for m = 10 tasks. w min is a synthetic parameter defined to instantiate problem instances. Essentially a larger w min value means longer tasks and longer data transfers, leading to a more "difficult" instance. The results show that Y-IE is the best or close to the best heuristic up to w min ≈ 8. For large values of w min , it is outperformed by other several other heuristics, such as P-IE, but also by the reference heuristic IE. IE is the best option for large values of w min ! An intuitive explanation is that when the instance is difficult, meaning that the probability of success is low due to long computations and communications, a good way to obtain a short makespan is to try to find the fastest workers and "hope for the best." When looking at the whole w min range, P-IE appears like a good alternative to Y-IE. For low w min values, it outperforms IE significantly, and for large w min it is outperformed by it only marginally. Recall from Table II that Y-IE and P-IE experience exactly the same number of failures (141 failures out of the 3,000 instances).

Unlike previous work that has considered looselycoupled master-worker applications, in this work a single processor failure can have a dramatic effect on application execution. The requirement that the application can progress only when all enrolled processors are simultaneously available, dramatically complicates all scheduling decisions. Furthermore, our problem formulation includes a limit on the available bandwidth from the master to the workers, and the possibility for processors to be temporarily reclaimed. We have proved the problem to be NP-complete in an off-line setting, i.e., with full knowledge of future processor states. By assuming a Markov model of processor availability, we have proposed polynomial time approximation schemes to compute the expected completion time of a computation, and its probability of success. We have then proposed 16 heuristics that are either passive (change the set of enrolled processors only when a processor failure occurs) or proactive (change the set of enrolled processors when a better set is available even if no failure occurs). These heuristics are easily defined as combinations of two among four sensible metrics: probability of success, expected completion time, expected yield and expected apparent yield.

All these heuristics widely outperform a baseline heuristic that allocates tasks to processors randomly. In addition, simulation results show that four of our 16 heuristics lead to significantly better results than the remaining 12. Passive heuristic IE is the most robust, which is why we have used it as a reference, but it does not lead to the best makespans. Heuristic Y-IE, which attempts to optimize expected execution time while proactively deciding to change the set of enrolled processors based on yield, leads to the best average results. Heuristic P-IE, which changes configuration based on probability of success, leads to more stable performance across our set of experimental problem instances as it is never significantly outperformed by IE. The conclusion is that a proactive heuristic that selects processors to maximize expected execution time and changes configuration based on yield or probability of success is very promising.

We have made it plain that the Markov assumption for processor availability is not meant to be representative of real-world platforms. Nevertheless, faced with the lack of an acknowledged and validated model, we have opted for a Markov model. The advantage of this model is that it provides us with a tractable framework in which we can not only develop heuristics but also evaluate the merit of heuristical ideas in "laboratory conditions". If a more realistic model arises, then a next step in this research would be to determine to which extent the principles from our heuristics can be adapted to the new model. Given that the results in Section V rely on the Markov assumption heavily (to provide accurate expressions for the expected completion time, probability of success and yield), it is unlikely that similar results would hold in a non-Markovian model. However, it may be possible to develop coarser estimates of our four criteria in the new model so as to design meaningful and effective heuristics. If no such new model arises, then an interesting next step would be to simply build a flawed Markov model based on real-world processor availability traces, and investigate how "wrong" the Markov heuristics behave in a real-world setting, when compared to heuristics that have been proposed in previous works and that do not rely on sophisticated probabilistic criteria for making scheduling decisions.

Figure 1 .

 1 Figure 1. Example iteration execution. White means UP , gray means RECLAIMED, and black means DOWN . P means "receiving the program," D means "receiving the data for a task," C means "computing," and I means "idle."

Theorem 4 . 1 :

 41 Problems OFF-LINE-COUPLED (µ = 1) and OFF-LINE-COUPLED (µ = ∞) are NP-hard.

  All processors are up at each step j such that |W | + 1 ≤ j ≤ N . Finally we let m = a and w = b + |W | + 1. Intuitively, this amounts to add |W | + 1 new vertices in W which are interconnected to every vertex in V . The size of the instance I 2 is linear in the size of the instance I 1 . We show that I 1 has a solution if and only if I 2 does. Suppose first that I 1 has a solution C = U 1 ∪ U 2 . We select the corresponding U 1 processors and the same U 2 time-steps, plus the last |W | + 1 time-steps. We have w = b + |W | + 1, hence I 2 has a solution. Suppose now that I 2 has a solution. The corresponding submatrix translates into a bi-clique with x processors and y time-steps. If x < m then at least one processor executes two tasks per iteration, and we need 2w timessteps to perform an iteration. But 2w > N , what is a contradiction. Hence x = m and y = K. At most |W | + 1 of the UP time-steps are greater than |W |, hence at least b of them are smaller than or equal to |W |: this leads to a solution to I 2 .

Theorem 5 . 1 :

 51 It is possible to approximate the values of P (S) +

+ ) W - 1 .

 1 We now explain how we numerically approximate the values of E u (S) and A(S). Let ε be the desired precision. Consider for some value T the difference between E u (S) and 0<t<T P

Figure 2 .

 2 Figure 2. %diff for m = 10 tasks vs. w min .

  E u (S) is the expected number of time slots with all processors UP , before one of these processors fails. Then, P

				converges. We can conclude that E( t>0 A t ) = t>0 P (S) u t →u . Then, (S) + = E u (S) -E u (S) × P (S) + , from which we derive that P (S) + = Eu(S) 1+Eu(S) if , in set
				S, at least one processor has a nonzero probability of
				going DOWN . Otherwise, P + = 1. (S)
				We now consider the expected time E (S) (W ) to
				execute W time slots of computation, conditioned by
				the fact that no processor in S will fail. The first time
				slot of computation is done at t = 0. Let E	(S) c	be
				the expected time of the next time slot of computation.
				Then,
				E	(S) c	= =	(S) + (t) t>0 t × P (S) t>0 t × P u t →u
					-t ×	0<t <t P + (t ) × P (S) u (S) t-t → u
					=	(S) u t>0 t × P t →u
					-	t>0 P + (t) × (S)	(S) u t >0 (t + t )P
	u	S) t →u	-t>0 P + (t) × t >0 P (S) u (S) t →u
	Let E u (S) = t>0 P u (S) t →u sors are UP at time slot 0. Let A t the random variable . Suppose that all proces-
	that is equal to 1 if all processors are UP at time slot
	t without that ant processor goes DOWN in between.
	Then E(A t ) = P we have E( 0≤t ≤t A t ) = 0≤t ≤t P (S) u t . By linearity of the expectation, →u (S) u
			(S)
			u	t →u

t →u . Suppose that, in set S, at least one processor has a nonzero probability of going DOWN . Then, lim t→∞ 0≤t ≤t P t →u Let A(S) = t>0 t × P

Table I RESULTS

 I WITH m = 5 TASKS.

	Heuristic	#fails	%diff	%wins	%wins30	stdv
	Y-IE	2	-11.82	72.58	92.09	0.42
	P-IE	2	-10.50	70.98	91.19	0.44
	E-IAY	4	-10.40	64.75	85.15	0.77
	E-IY	4	-3.40	59.91	81.64	0.80
	IE	1	0.00	100.00	100.00	0.00
	IAY	2	13.59	51.07	76.42	1.93
	E-IP	4	19.35	47.73	69.69	0.98
	IY	2	24.22	45.26	70.85	1.96
	IP	2	52.03	34.79	58.54	2.11
	E-IE	5	53.93	39.57	64.51	2.57
	Y-IAY	3	99.75	53.89	70.77	5.55
	Y-IY	3	113.01	49.22	66.80	5.73
	P-IAY	3	125.27	50.28	67.33	6.08
	Y-IP	2	145.05	38.56	55.54	5.90
	P-IY	3	145.78	42.54	59.66	6.22
	P-IP	2	176.92	36.92	52.00	6.61
	RANDOM	0	2124.42	0.00	0.20	22.54

Table II RESULTS

 II WITH m = 10 TASKS FOR THE BEST EIGHT HEURISTICS.

	Heuristic #fails	%diff	%wins	%wins30	stdv
	Y-IE	141	-10.33	71.35	88.42	0.54
	P-IE	141	-8.62	69.64	87.23	0.55
	E-IAY	178	-6.10	66.62	81.93	1.58
	E-IY	176	8.04	61.90	77.87	3.07
	E-IP	168	29.68	55.12	71.86	3.01
	IAY	152	136.65	46.98	69.31	14.76
	IY	152	147.77	42.06	64.47	14.76

Acknowledgments.: Y. Robert is with the Institut Universitaire de France. This work was supported in part by the ANR RESCUE project.