Julie Delon

Antoine Houdard

Gaussian Priors for Image denoising

come

Gaussian Priors for Image denoising

Introduction

This chapter focuses on patch priors for image denoising. In the last decade, patch-based models (also known as Non-local models) have created a new paradigm in image processing, leading to very significant improvements both for classical image restoration problems (denoising, inpainting, interpolation) or for image synthesis and editing. These models represent images by a set of local neighborhoods or patches, and make them collaborate regardless of their spatial position in the image, relying on the observation that most natural images present a remarkable redundancy at a semi-local scale. A patch y i (v) is a piece (most of the time a square) of an image v centered at the pixel i. As pointed out by Mumford and Desolneux [START_REF] Mumford | Pattern theory: the stochastic analysis of real-world signals[END_REF], patches are "the analogs of the phonemes of speech".

Patch-based models have been the subject of numerous works, especially in the context of image denoising. Assuming that the noise is additive, image denoising amounts to estimate an image u from its noisy version v ∈ R m (m is the image size) such that

v = u + ε, (1)
with ε a noise with known statistics (not necessarily Gaussian). In digital cameras, the two major sources of noise during the acquisition process are the thermal agitation, which produces an almost white and Gaussian noise, and the discrete nature of light, which is behind the photon shot noise, modeled as a Poisson variable (for a complete description of the sources of noise in a digital camera, see [START_REF] Aguerrebere | Study of the digital camera acquisition process and statistical modeling of the sensor raw data[END_REF]). Stabilizing the noise variance by a generalized Anscombe transform [START_REF] Makitalo | A closed-form approximation of the exact unbiased inverse of the Anscombe variance-stabilizing transformation[END_REF] results in a noise model well approximated by a white Gaussian noise ε ∼ N (0, σ 2 I m ). The vast majority of works on image denoising focus on this simplified model and it is also our assumption in this chapter.

In this framework, patch-based methods usually attempts at rewriting (1) into a degradation model that can be expressed for each patch separately. All patches {y i , i = 1 . . . , m} of size p = s × s are first extracted Figure 1: Image patches can be seen as vectors in a high-dimensional space. Most of the patch-based methods uses the patch-space of an image which is the set of all the sliding patches of size p = s × s extracted from the image. from the image v and seen as noisy vectors in a high dimensional space, as illustrated by Figure 1 (in the whole chapter, when writing patches as vectors, we assume that the patches are read column-wise). Then the noisy patches are restored sequentially, before reconstructing the whole image. The degradation model on the patches becomes

y i = x i + ε i , i ∈ {1, . . . , m} (2)
where x i is the patch centered at pixel i in u, y i the same patch in v, and ε i the additive noise. In practice, it is almost always assumed that the {ε i , i = 1 . . . , m} are independent samples from the Gaussian distribution N (0, σ 2 I p ), although this hypothesis is obviously wrong since patches are overlapping. We will briefly discuss this issue in Section 4, along with the aggregation of the restored patches to reconstruct the whole image.

The first denoising methods relying on patches appear in 2004 [START_REF] Ordentlich | A discrete universal denoiser and its application to binary images[END_REF][START_REF] Weissman | Universal discrete denoising: Known channel[END_REF][START_REF] Awate | Image denoising with unsupervised information-theoretic adaptive filtering[END_REF][START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF]. Among these methods, one of the most popular remains the Non-Local Means [START_REF] Buades | A review of image denoising algorithms, with a new one[END_REF], which sees similar patches as independent realizations of the same distribution and averages these repeated structures to reduce noise variance. If numerous approaches have built on the same core ideas since 2004, the recent and most convincing approaches in patch-based denoising rely on a Bayesian reformulation of the denoising problem, using local or global statistical priors for the distribution of each patch [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: from gaussian mixture models to structured sparsity[END_REF][START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF][START_REF] Aguerrebere | A bayesian hyperprior approach for joint image denoising and interpolation, with an application to hdr imaging[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF]. Under the white Gaussian noise model (2), the conditional distribution of a noisy patch y knowing its original version x (we omit the index i in the following) can be written

p(y|x) ∝ e -x-y 2 2σ 2 .
(

) 3
The Bayesian model assumes that the original patch x is a realization of a random vector X with a probability distribution p(x) called the prior distribution. Therefore, the noisy patch y is a realization of the random vector

Y = X + N, (4)
with N ∼ N (0, σ 2 I p ). Under these hypotheses, and assuming that N and X are independent, we can compute the posterior distribution

p(x|y) ∝ p(y|x)p(x) ∝ e -x-y 2 2σ 2 p(x). (5)
Ideally, in order to reconstruct the (unknown) original patch x from the degraded version y, we would like to compute the conditional expectation E[X|Y ] (i.e. the mean of the posterior distribution), which minimizes the quadratic risk under the previous model. This estimator is also called the minimum mean square error (MMSE) estimator. In practice, computing this conditional expectation is often complex, and it is classical to compute instead the affine function (called linear MMSE) of Y minimizing the quadratic risk, i.e. the affine estimator DY + α (with D a p × p real matrix and α a vector in R p ) minimizing the risk

E[ DY + α -X 2 ].
This affine estimator, is called the Wiener estimator and will be denoted E Wiener [X|Y ] in the following. It can be easily shown by deriving the previous risk that (assuming that the following quantities exist),

E Wiener [X|Y ] = E[X] + Σ X,Y Σ -1 Y (Y -E[Y ]), (6)
where

Σ X,Y := E [(X -E[X])(Y -E[Y ]) t ] and Σ Y := E [(Y -E[Y ])(Y -E[Y ]) t ].
This affine estimator only relies on second-order moments of the signal and noise. Under model ( 4) and assuming that N and X are independent, the Wiener estimator becomes

E Wiener [X|Y ] = E[X] + Σ X (Σ X + σ 2 I p ) -1 (Y -E[Y ]), (7)
with Σ X the covariance matrix of the random vector X.

Another classical solution to reconstruct x is to compute the maximum (MAP) of the a posteriori distribution p(y|x), which yields: xy 2 2σ 2log p(x).

From this point of view, restoring each patch is equivalent to solve a variationnal problem, with a quadratic fidelity term and a smoothness term derived from the prior.

The most convenient prior for computing the previous estimators is the Gaussian distribution. Indeed, on the one hand, Gaussian priors are well suited to encode patch structures with some kind of contrast invariance, as we will see in Section 2. On the other hand, under a Gaussian prior, the conditional expectation, Wiener estimator and MAP coincide, as we will see in Section 3. For these reasons, these priors are favored in most recent works on patch-based image denoising [START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Aguerrebere | A bayesian hyperprior approach for joint image denoising and interpolation, with an application to hdr imaging[END_REF]. A slightly more involved prior used in the literature is the Gaussian Mixture Model (GMM) [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF][START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: from gaussian mixture models to structured sparsity[END_REF][START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF]. In this case, computing the conditional expectation remains simply tractable. All these works differ among other things in the way they infer the parameters of the Gaussian or GMM distributions. These distributions live in R p and estimation in such high-dimensional spaces is complex. We will see in Section 5 the different possibilities to infer these parameters and how some of these works tackle the curse of dimensionality. Figure 2 illustrates the main steps common to all these patch based denoising methods, and each of these steps is described in the following sections. 2 What is encoded in Gaussian and GMM priors ?

Before going into the details of estimation under Gaussian priors, we provide in this section a few insights on the actual structures they encode. Assume a Gaussian model N (µ, Σ) for p = s × s patches (µ ∈ R p and Σ ∈ M p (R)). The diagonal coefficients of the covariance matrix Σ represent the variance of each pixel in the patch, while the non-diagonal coefficients represent the covariances between pixels. A positive covariance coefficient means that the two pixels tend to be either both greater or smaller than their means, while a negative coefficient implies that they tend to be on opposite sides of their means. Clearly, if Σ is purely diagonal, patches drawn from the model N (µ, Σ) will only be noisy versions of the mean patch µ. In this case, the only structure information is contained in µ. More interesting models contain geometric information directly in the covariance matrix Σ. To illustrate this point, we propose to create models encoding different patch structures. For instance, in order to model a vertical edge, we define a Gaussian distribution with constant mean µ = (0.5, . . . , 0.5) and a covariance matrix with coefficient 1 in the second and third quarter of Σ, and coefficient 0 in the first and fourth quarters of Σ (see Figure 3). In this simplistic example, the matrix Σ has rank two, with (non trivial) eigenvectors (1, . . . , 1, 0, . . . , 0) and (0, . . . , 0, 1, . . . , 1), so all the patches drawn from this distribution can be written 0.5 + (α, . . . , α, β , . . . , β ) with α ∼ N (0, 1) and β ∼ N (0, 1). These patches all contain a vertical edge in their middle, with grey levels α and β on both sides of the edge. In this example, we see that the model encodes a structure and authorizes different contrasts on both sides of the structure. With the same mechanic, we can create a covariance matrix encoding any desired shape, see for instance Figure 4. Again, the samples from the corresponding distribution exhibit all possible grey levels in the different regions defined by the covariance matrix, even if all these grey levels are not all equally likely. Now, although these models authorize contrast changes or contrast inversions, they are not well suited to encode geometric invariances on patches. For instance, if we try to learn a model encoding different vertical edges with invariance to translation, we end up with an average model encoding a vertical gradient image (see Figure 5).

How to derive estimators under Gaussian and GMM priors

Now that we have seen more precisely what could be contained in Gaussian priors, we will now see more precisely how they can be used to derive estimators under the Bayesian model described in the introduction.

In the whole section, we assume that we work with the model ( 4)

Y = X + N,
with N ∼ N (0, σ 2 I p ) independent from X. We wish to estimate X knowing Y . We first recall some classical results on the conditioning of Gaussian vectors, and on the links between the conditional expectation, Wiener estimator and MAP for Gaussian and GMM priors. These different estimators will serve in the rest of the chapter as denoising strategies for image patches.

Estimation with Gaussian priors

We first assume that X follows a Gaussian distribution N (µ X , Σ X ) and that the noise N is independent from X. The classical properties of Gaussian vectors make it possible to show that in this case the estimator E[X|Y ] is an affine function of Y (thus equivalent in this case to the Wiener estimator). Indeed, recall that if (T,V ) is a Gaussian vector, then the conditional expectation

E[T |V ] is the affine function of V E[T |V ] = E[T ] + Σ T,V Σ -1 V (V -E[V ]), (8)
where

Σ V is the covariance matrix of V and Σ T,V = E[(T -E[T ])(V -E[V ]) t ] (if Σ V is not full rank, the
result is still true by taking the Moore-Penrose pseudo-inverse of Σ V ). Now, if X and N are independent Gaussian random vectors, the concatenated vector (X,Y ) = (X, X + N) is also Gaussian. We directly deduce the following result. Proposition 1. Assume that X and Y follow the model (4), with X ∼ N (µ X , Σ X ) and N ∼ N (0, σ 2 I p ) independent, then the conditional expectation and Wiener estimator of X knowing Y coincide and can be written Proof. On the one hand, since (X,Y ) is a Gaussian vector, the conditional expectation E[X|Y ] can be written

E[X|Y ] = E Wiener [X|Y ] = µ X + Σ X (Σ X + σ 2 I p ) -1 (Y -µ X ).
E[X|Y ] = E[X] + Σ X,Y Σ -1 Y (Y -E[Y ]) = E[X] + E[(X -E[X])(X + N -E[X + N]) t ](Σ X + σ 2 I p ) -1 (Y -E[Y ]). = E[X] + Σ X (Σ X + σ 2 I p ) -1 (Y -E[Y ]) = µ X + Σ X (Σ X + σ 2 I p ) -1 (Y -µ X ).
Under the same hypothesis, if we try to maximize the a posteriori probability on the patch X, we obtain arg max

X log P[X|Y ] = arg max X (log P[Y |X] + log P[X]) = arg min X (X -Y ) t (X -Y )/σ 2 + (X -E[X]) t Σ -1 X (X -E[X]) .
We check easily that the solution of this minimization problem is also given by

ψ(Y ) = µ X + Σ X (Σ X + σ 2 I p ) -1 (Y -µ X ).
Said otherwise, for a Gaussian prior, the MMSE, linear MMSE and MAP all coincide and all these estimators only require linear operations. This property makes Gaussian priors particularly convenient in practice and explains their success in the restoration literature.

We can illustrate the interest of this estimator on the Gaussian model N (µ X , Σ X ) presented on Figure 3 and representing a vertical edge. If X is an (unknown) realization of this model and Y = X + N with N ∼ N (0, σ 2 I p ) independent from X, then E[X|Y ] will also be a patch (α, . . . , α, β , . . . , β ) with

α = 0.5 + 1 p/2+σ 2 ∑ p/2 k=1 (Y k -0.5) and β = 0.5 + 1 p/2+σ 2 ∑ p k=p/2+1 (Y k -0.5
)(assuming p is even for the sake of simplicity). Said otherwise, the denoised patch E[X|Y ] represents the same vertical edge as X and its values α and β on both sides of the edge are (if σ 2 << p/2) the averages of Y on these two half patches. Figure 6 represents three denoising experiments with the previous estimator. On the first line, a vertical edge is denoised with the Gaussian model of Figure 3. On the second line, a "duck" patch is denoised with the Gaussian model of Figure 4. In both cases, using the conditional expectation works extremely well because the Gaussian model used in the estimator fits perfectly the image to be denoised. On the third line, the noisy edge is denoised with the Gaussian model of Figure 5. In this case, the denoised patch is constant on each column (since the model is learned from a set of translated vertical edges). Although the model imposes a strong correlation between columns of the first half of the patch on the one hand, and between columns of the second half of the patch on the other hand, this is not enough to restore the patch perfectly. 3 is used to denoise (PSNR = 37.17). Second line, the duck Gaussian model of Figure 4 is used to denoise (PSNR = 34.29). Third line, the gradient model of Figure 5 is used to denoise (PSNR = 29.68). In this last case, the image to be denoised is not well represented by the model and the result is less convincing.

Estimation with Gaussian Mixture Models

The case of Gaussian Mixture Models is a bit more involved but remains globally simple. Assume that X follows a Gaussian Mixture Model

X ∼ K ∑ k=1 π k N (µ k , Σ k ), (9)
with ∑ K k=1 π k = 1. There exists a latent random variable Z on {1, . . . , K} such that P[Z = k] = π k and such that X|Z = k ∼ N(µ k , Σ k ). In the following, we note ψ k (y) the Wiener estimator for the k th Gaussian, i.e.

ψ k (y) = µ k + Σ k (Σ k + σ 2 I p ) -1 (y -µ k ).
Under this model, we have the following proposition.

Proposition 2. Assume that X and Y follow the model (4), with X ∼ ∑ K k=1 π k N (µ k , Σ k ) and N ∼ N (0, σ 2 I p ) independent, then the conditional expectation of X knowing Y can be written

E[X|Y ] = K ∑ k=1 ψ k (Y )P[Z = k|Y ]. (10)
Proof. To compute the conditional expectation, we can start by noting that if Z = k, (X,Y ) is a Gaussian vector and the results of the previous section apply. We can now compute the conditional expectation

E[X | Y, Z] = ψ Z (Y ) = K ∑ k=1 ψ k (Y )1 Z=k .
It follows that

E[X|Y ] = E[E[X | Y, Z] | Y ] because σ (Y ) ⊂ σ (Y, Z) = E[ψ Z (Y ) | Y ] = K ∑ k=1 E[ψ k (Y )1 Z=k | Y ] = K ∑ k=1 ψ k (Y )E[1 Z=k | Y ] because ψ k (Y ) is σ (Y )-measurable.
We deduce that

E[X|Y ] = K ∑ k=1 ψ k (Y )E[1 Z=k | Y ] = K ∑ k=1 ψ k (Y )P[Z = k | Y ].
The conditional expectation E[X|Y ] can be seen as a linear combination of affine functions of Y , with weight P[Z = k|Y ] representing the probability that the patch belongs to the class k. However, the weights

P[Z = k | Y ] are not linear functions of Y .
The expression of the Wiener estimator E Wiener [X|Y ] can be deduced directly from Equation [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF], by replacing E[X] by ∑ K k=1 π k µ k and Σ X by the complete covariance of the GMM. Finally, computing the MAP arg max X log P[X|Y ] under a GMM prior on X is much less convenient and does not lead to a closed-form solution. Indeed, it boils down to compute the maximum of the posterior distribution, which is another Gaussian Mixture distribution.

In other words, the linear MMSE, MMSE and MAP do not coincide for Gaussian Mixture priors. In practice, the conditional expectation is favored since it is much simpler to compute than the MAP.

Other estimation strategies

Estimation under Gaussian or GMM models has several links with other estimation strategies found in the literature. For a noisy patch y, and a Gaussian model N (µ, Σ), we have seen that the conditional expectation strategy consists in computing the denoised patch

x(y) = µ + Σ(Σ + σ 2 I p ) -1 (y -µ).
Now, if we consider the eigendecomposition Σ = Q∆Q t with ∆ = diag(λ 1 , . . . , λ p ), this can be rewritten

x(y) = µ + Qdiag λ 1 λ 1 + σ 2 , . . . , λ p λ p + σ 2 Q t (y -µ). (11)
More generally, denoting Q 1 , . . . , Q p the columns of Q representing the eigenvectors, we can write

x(y) = µ + p ∑ k=1 η k ( Q k |y -µ ) Q k , (12)
with η k (z) = λ k λ k +σ 2 z. Although the previous Wiener estimator is used in numerous recent patch-based denoising methods [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF], other choices are obviously possible for η k , such as hard or soft thresholding [START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF], or all estimators classically used in diagonal estimation.

Writing x = Q t (xµ), we can see that the conditional expectation x(y) is also solution of the optimization problem

argmin x Q x -(y -µ) 2 + σ 2 xt ∆ -1 x = argmin x Q x -(y -µ) 2 + σ 2 p ∑ k=1 x2 j λ k .
This permits to see the link between the previous approach and the dictionary-based approaches, the dictionary here being given by Q and the second term corresponding to a regularization of the solution x. Figure 7 represents the denoising of a noisy patch with the same Gaussian model and two different denoising strategies: the conditional expectation (Wiener) and hard thresholding at 2.7σ (as recommended in [START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF]). 4 From patches to images: aggregation procedures

In the previous sections, we have seen how to derive bayesian estimators to perform denoising on each patch separately. In this framework, each observed patch y i from a noisy image v is denoised into x i , which is an estimate of the unknown patch x i . Each pixel of the image v is contained in p patches, which provide p denoised versions for this pixel. Most aggregation procedures consists in defining a reprojection function ψ : R m×p → R m which reconstructs an image from the set of its denoised patches. Observe that since denoised patches usually do not coincide on their overlap, this operation is not invertible. Moreover, since the noise on overlapping patches is not independent, the p denoised versions of the pixel carry this dependence under the form of low-frequency noise. In the literature, we find three main strategies for this reprojection step:

• Central pixel reprojection. The idea is to keep only the central pixel of each denoised patch.

• Uniform reprojection. All the estimators coming from the different patches containing the pixel are averaged with uniform weights. This strategy is the most commonly used in practice, and this is the one we use in this chapter for the sake of simplicity.

• Weighted reprojection. All the estimators coming from the different patches containing the pixel are averaged with weights representing the precision of the corresponding estimator. For some details see [START_REF] Salmon | From patches to pixels in non-local methods: Weighted-average reprojection[END_REF][START_REF] Pierazzo | Multi-scale dct denoising[END_REF][START_REF] Dabov | Image denoising by sparse 3-d transform-domain collaborative filtering[END_REF].

A more involved strategy is explored in [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. The authors propose to reconstruct the denoised image u as the solution of

argmin u λ 2 u -v 2 2 -∑ j log p(x j ),
where the {x j } are the patches extracted from the unknown image u and p is a GMM prior on the image patches. This formulation includes both the denoising and aggregation step into a single variational problem.

Inference of Gaussian and GMM priors

Gaussian models and GMMs appear to be well suited for patch based denoising. However, the quality of the restoration strongly depends on the relevance of the model. Unfortunately, in real denoising problems the perfect model is never known and the most challenging step is to find a good prior for each patch. In the literature, we find essentially two strategies to learn these models. The first one consists in learning the model on some external set of patches that represent the diversity of natural images [START_REF] Zoran | From learning models of natural image patches to whole image restoration[END_REF]. The second one consists in learning the model directly on the noisy patches [START_REF] Afonso M Teodoro | Single-frame image denoising and inpainting using gaussian mixtures[END_REF][START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF]. In this section, we discuss different approaches adopting the second strategy. Before going further, we recall some basics about statistical inference. Given a set of patches {y 1 , . . . , y n } ∈ R p extracted from an image, we consider them as independent realizations of a random variable Y with density φ depending on some parameters θ . The parameters θ of the model are inferred by maximizing the likelihood of the data w.r.t. θ , where the likelihood is defined as

(y; θ ) = n ∏ i=1 φ (y i ; θ ). ( 13
)
Maximizing the likelihood is equivalent to minimize the negative log-likelihood

L (y; θ ) = -log ( (y; θ )) = - n ∑ i=1 log (φ (y i ; θ )) , (14)
which is usually more convenient for computation.

In the context of denoising, we put a prior model on the random vector X representing the clean patches. When X follows a Gaussian model of parameters (µ X , Σ X ), resp. a Gaussian mixture model of parameters {π k , µ k , Σ k } k=1...K , then Y = X + N also follows a Gaussian model of parameters {µ X , Σ X + σ 2 I} k , resp. a GMM of parameters (π k , µ k , Σ k + σ 2 I). Since Σ X (resp. Σ k ) is positive semi-definite and σ > 0, Σ X + σ 2 I (resp. Σ k + σ 2 I) is always positive definite. Thus, the random vector Y always has a probability density function φ and the likelihood is always defined.

Gaussian models

In the case of a Gaussian prior X ∼ N (µ X , Σ X ) on the clean patches, the set of parameters on the noisy patches is given by θ = {µ Y , Σ Y } where Σ Y = Σ X + σ 2 I and µ X = µ Y . The negative log-likelihood for a set of noisy data {y 1 , . . . , y n } becomes

L (y; θ ) = 1 2 n ∑ i=1 (y -µ Y ) T Σ Y -1 (y -µ Y ). ( 15
)
The computation of the maximum likelihood estimators (MLE) of the parameters, i.e. argmin θ L (x; θ ), for µ Y and Σ Y yields the sample mean

µ Y (n) = 1 n n ∑ i=1 y i , (16)
and the sample covariance matrix

Σ Y (n) = 1 n n ∑ i=1 (y i -µ Y ) T (y i -µ Y ). (17)
Theses estimators depend on the number n of samples and from the strong law of large numbers

µ Y (n) a.s. -→ n→∞ µ Y and Σ Y (n) a.s. -→ n→∞ Σ Y . (18)
This gives us an estimator

Σ X := Σ Y -σ 2 I for Σ X satisfying Σ X (n) a.s. -→ n→∞ Σ X . (19)
In summary, for a given set of noisy patches {y 1 , . . . , y n } we can easily compute the MLE of the parameters (µ X , Σ X ) for the Gaussian model on the underlying clean patches. Now, since we showed in Section 2 that Gaussian models are representing really precise structures, the most challenging part is to choose the set of noisy patches from which the model can be derived.

How to group patches to infer Gaussian priors?

In this section, we discuss how patches can be grouped in order to learn the previous Gaussian models directly from a noisy image.

Global Gaussian prior

The first really basic idea is to model the set of all image patches with a unique Gaussian prior. In this case, we are modeling the whole "patch-space" by a unique Gaussian model of mean µ X and covariance Σ X . This model poorly represents the complexity of the patch-space but still encodes some proper image information. This modeling is adopted in [START_REF] Deledalle | Image denoising with patch based pca: local versus global[END_REF] to perform a basic denoising by performing the eigendecomposition Σ X = Q∆Q t and denoising the patches with an estimator of the form [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF]. Figure 8 illustrates the fact that the eigenvectors of the covariance matrix learned on the whole patch space encode some proper information about the image. In this case, since the Gaussian model is very broad, we do not expect the Wiener estimator to yield good results. But since the eigenbasis seems to encode some proper information about the image patches, the hard thresholding strategy manages surprisingly good denoising. The second line of Figure 9 shows the denoising result for this global grouping with the two denoising strategies and shows that in this case, the hard-thresholding strategy is better than the Wiener one.

Spatially local Gaussian priors

To derive more precise prior models, it is necessary to group "similar" patches and to restrict the inference to each of these groups. A first possibility is to group patches based on their spatial proximity in the image. This makes sense in homogeneous regions, but the risk is high to group patches representing really different structures. The third line of Figure 9 shows that the result of this strategy is not really better, PSNR-wise, than the result of the global strategy. However, the Wiener strategy for this local approach seems nicer than in the global approach, while the result of the hard-thresholding strategy does not really change.

Local Gaussian priors in the space of patches

In order to learn more precise models, patches can be clustered directly in the patch space and a Gaussian model can be inferred for each cluster. All patches from the cluster can then be denoised using this model. This clustering implies to use an appropriate similarity measure between patches. The fourth line of Figure 9 shows such a denoising experiment with a K-means clustering relying on the Euclidean distance, with K = 256 clusters (Figure 10 shows the corresponding clustering). This usually yields a better denoising than the global and the local grouping strategies.

This way of grouping patches in the patch space together with a Wiener filtering is also one of the main ideas behind the two steps of the NL-Bayes algorithm [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF]. In this algorithm, each patch y i is associated with the group of all its ε-close patches for the Euclidean norm. A Gaussian model is inferred from this group and the whole group is denoised using this model. The final estimator for each patch is the average of all its denoised versions. The NL-Bayes algorithm uses this strategy twice: in the first step, distances are computed directly between noisy patches in R p ; in the second step, distances between patches are computed between the versions which have denoised during the first step. Grouping ε-close patches presents the advantage of putting together patches representing the same structures. However, a straightforward one-step implementation (fifth row of Figure 9) of this idea shows that it does not work as well as expected in practice. Two major issues arise in this context:

• The high dimensionality of the patch space makes the estimation of the covariance matrix difficult;

• The use of the Euclidean distance for grouping does not allow similar patches with different contrast to be in the same group, which is a loss because we saw in Section 2 that a Gaussian model can encode information up to contrast changes.

The first issue, discussed in Section 5.4, is crucial and related to the curse of dimensionality. Unfortunately, it is hardly taken into account in the image denoising literature.

To tackle the second issue, other norms were investigated in the literature [START_REF] Deledalle | How to compare noisy patches? patch similarity beyond gaussian noise[END_REF]. Another idea is to use the Gaussian models previously learned for recalculating new clusters. Indeed, each covariance matrix of the different Gaussian models provides a semi-norm that can be used to recompute the ε-nearest patches of each group.

Inference for Gaussian Mixture Models

The inference in the case of a mixture model is slightly more challenging since a direct maximization of the likelihood is not possible. The negative log-likelihood of the noisy data {y 1 , . . . , y n } is given by

L (y; θ ) = n ∑ i=1 log K ∑ k=1 π k φ (y i ; θ k ) (20)
and the minimization of this function w.r.t θ is a complex problem. However, if we know to which group each sample x i belongs, the log-likelihood becomes

L (y, z; θ ) = n ∑ i=1 K ∑ k=1 z ik log (π k φ (y i ; θ k )) (21)
with z ik = 1 if y i belongs to the group k and 0 otherwise. L (y, z; θ ) is the log-likelihood of the data completed with the latent random variable Z that determines the group from which the observations come from, that is

Y i |(Z i = k) ∼ N (µ k , Σ k ) and p(Z i = k) = π k .
The EM algorithm consists in iterating two steps ; the expectation (E) step that calculates the expected value of ( 21) with respect to the conditional distribution of Z given Y for the current value of the parameters θ . And the maximization (M) step that consists in the update of the parameters by minimizing the expectation of the complete log-likelihood from the E-step:

E (L (y, z; θ )) = n ∑ i=1 K ∑ k=1 E(z ik |x i , θ ) log (π k φ (y i ; θ k )) (22)
which leads to tractable expressions for the MLE of the parameters. It can be shown (see for example [START_REF] Christopher | Pattern recognition and machine learning[END_REF]) that this algorithm converges to a local minimum of the log-likelihood [START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF].

In the precise case of a Gaussian mixture model, the two steps of the algorithm become

• E-step, computation of t ik := E(z ik |y i , θ ) t ik = π k φ (y i ; θ k ) ∑ K l=1 π l φ (y i ; θ l ) (23)
• M-step, update of the parameters

π k = 1 n n ∑ i=1 t ik , (24)
µ k = ∑ n i=1 t ik y i ∑ n i=1 t ik , (25)
Σ k = ∑ n i=1 t ik (y i -µ k )(y i -µ k ) T ∑ n i=1 t ik . (26)
Observe that if we impose the t ik to be 1 when the patch i belongs to the group k and 0 otherwise, the M-step consists in inferring the parameters of the Gaussian models for the groups, while the E-step uses the knowledge of the inferred model to update the groups themselves. This model provides a better clustering of the patches than a K-means clustering with the Euclidean norm (which only produces isotropic clusters) and consequently should yield better denoising results. This idea is used in [START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF][START_REF] Yang | Compressive sensing by learning a gaussian mixture model from measurements[END_REF][START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF] and the GMM model on patches is also used in [START_REF] Yu | Solving inverse problems with piecewise linear estimators: from gaussian mixture models to structured sparsity[END_REF]. A straightforward implementation of the denoising with a GMM model on the patches gives the result in the first line of Figure 11. However, this inference of a GMM also strongly suffers from the curse of the dimensionality and algorithms such S-PLE [START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF] or HDMI [START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF][START_REF] Houdard | Clustering en haute dimensions pour le débruitage d'image[END_REF] propose to use Gaussian Mixture models with intrinsic lower dimensions in order to reduce the number of parameters to estimate, as detailed in the following section.

Inference in high dimension

The dimensions of the patch spaces are usually high, from p = 9 (for 3 × 3 patches) to p = 100 for 10 × 10 patches, or even higher. Estimating the parameters of Gaussian models (or GMM) in such high dimensional spaces is complex. When p is large, patches seen as points in R p are essentially isolated, the euclidean distance and the notion of nearest neighbor become much less reliable than in low dimensional spaces [START_REF] Giraud | Introduction to high-dimensional statistics[END_REF]. These phenomena, known as the curse of dimensionality, cause difficulties to decide which patches should be grouped together in a common Gaussian model. Besides, parametric models such as Gaussian Mixture Models in high-dimension are usually over-parametrized: the covariance matrix of a Gaussian model in dimension p = 100 contains 5050 different coefficients. They necessitate huge quantities of data to be estimated correctly. Indeed, the convergence of the sample covariance matrices to the true covariance matrix depends on the ratio between the number n of samples and the dimension p. More precisely, if n and p both tend toward infinity while n p tends toward a constant c > 0, the eigenvalues of the sample covariance matrix Σ(n) do not necessarily converge towards the eigenvalues of the model covariance matrix (Marcenko-Pastur Theorem [START_REF] Vladimir | Distribution of eigenvalues for some sets of random matrices[END_REF] describes the limit law of the empirical distribution of these eigenvalues).

A consequence of the curse of dimensionality is that clustering methods such as K-means of GMM are often disappointing in high dimension, or do not converge at all if p is too large. Solutions to circumvent these problems usually rely on dimension reduction, or regularization of the model parameters. For instance, if the sample covariance matrix Σ is singular of ill-conditioned, or is not definite positive, it is usual to add a small εI p to it. This is the strategy followed by [START_REF] Lebrun | A Nonlocal Bayesian Image Denoising Algorithm[END_REF][START_REF] Yu | Solving inverse problems with piecewise linear estimators: from gaussian mixture models to structured sparsity[END_REF]. In the case of Gaussian Mixture Models, another approach consists in assuming that the intrinsic dimension of the Gaussian is lower than p. This is the idea adopted in [START_REF] Wang | SURE Guided Gaussian Mixture Image Denoising[END_REF], where the groups intrinsic dimensions are heuristically fixed to 1 (flat regions), p 2 or p-1. A more involved method consists in inferring for each group its own intrinsic dimension [START_REF] Houdard | High-Dimensional Mixture Models For Unsupervised Image Denoising[END_REF] (see Figure 11). The corresponding parsimonious model assumes that each Gaussian of the mixture lives in its own specific subspace.

Conclusion

In this chapter, we have focused on patch priors for image denoising. As we have seen, assuming Gaussian and GMM priors on image patches is now quite common in the restoration literature. We have tried to provide a unified point of view for all of these methods, in order to underline their similarities and differences. Table 1 summarizes the main features of the methods mentioned in this chapter. We have also described some of their limitations, such as the inference difficulties in high dimension or the absence of invariance properties to geometric transformations. We did not discuss the computational cost of these approaches, but this point is clearly a critical issue for industrial applications.

  x(y) = arg max x∈R p p(x|y) = arg max x∈R p p[y|x] p(x) = arg min x∈R plog p(y|x)log p(x) = arg min x∈R p

Figure 2 :

 2 Figure 2: The whole process of patch-based image denoising with Gaussian prior models. First, patches are extracted from the noisy image. Next, these noisy patches are grouped and modeled with local Gaussians or Gaussian Mixture Models, whose parameters are inferred by maximum likelihood (Section 5). Each patch is then denoised with an estimator derived from the model (Section 3). Finally, the clean patches are aggregated to recover the denoised image (Section 4).

Figure 3 :

 3 Figure 3: Left: a covariance matrix Σ with 1 (white) on the second and third quarters, and 0 (black) on the first and fourth quarters. Right: patches drawn from the Gaussian distribution N (µ, Σ) with µ a constant patch equal to 0.5.

Figure 4 :

 4 Figure 4: Left: a covariance matrix Σ composed of 1 (white) and 0 (black). Right: patches drawn from the Gaussian distribution N (µ, Σ) with µ a constant patch equal to 0.5.

Figure 5 :

 5 Figure 5: Left: a covariance matrix Σ learned as the sample covariance matric of a set of vertical edges at different spatial positions, and with also different choices of grey levels on both sides of the edge. Right: patches drawn from the corresponding Gaussian distribution N (µ, Σ) with µ a constant patch equal to 0.5.

Figure 6 :

 6 Figure 6: For each line, from left to right, clean patch, noisy patch (σ = 10%), denoised patch with the Wiener estimator. First line, the edge Gaussian model of Figure 3 is used to denoise (PSNR = 37.17).Second line, the duck Gaussian model of Figure4is used to denoise (PSNR = 34.29). Third line, the gradient model of Figure5is used to denoise (PSNR = 29.68). In this last case, the image to be denoised is not well represented by the model and the result is less convincing.

Figure 7 :

 7 Figure 7: Clean patch, noisy patch (10% noise), denoised patch with gradient model (from Fig. 5) and Wiener estimator (PSNR = 29.68dB), and denoised patch with gradient model and hard thresholding (PSNR = 31.12dB, th = 2.7σ )

Figure 8 :

 8 Figure 8: Visualization of the first 16 eigenvectors of the sample covariance matrix of the whole patch space for two different images. Left: original images. Middle: the 16 first eigenvectors. Right: patches generated with the low rank covariance matrix created from these eigenvectors.

Figure 10 :

 10 Figure 10: Left: the local grouping used in the local strategy. Middle and Right: the grouping used in the K-means strategy for the two images Simpson and Alley.

Figure 11 :

 11 Figure 11: First line: Denoising with a full GMM model (50 groups) on all the patches. The clustering (left) is quite noisy and the denoising result (right) is not very good (PSNR: 28.50dB). Second line: Denoising with a GMM model (50 groups) with intrinsic dimension regularization as in [11]. The clustering (left) is smoother and the denoising yields quite good results (PSNR: 31.23dB)


Tianyu Wang
email: tianyu.wang@emn.fr

• Odile Bellenguez-Morineau

Odile Bellenguez-Morineau

Complexity of parallel scheduling unit-time jobs with in-tree precedence constraints while minimizing the mean flow time

Keywords: parallel scheduling, in-tree, precedence constraints, complexity theory

This paper deals with a particular scheduling problem. We consider unit-time jobs and in-tree precedence constraints while minimizing the mean flow time. This problem is observed as P |p j = 1, in-tree| C j with the use of the 3-filed notation. To the best of our knowledge, its complexity is still open. Through a reduction from 3-Partition, we show that this problem is N P-complete.

to be different for specific type of precedence graphs. In this paper, we focus on the problem with in-trees (or in-forest in some literature) precedence graph. This particular graph suggests that each job has no more than one successor. We use C j to denote the completion time of a job j. The problem involves finding an optimal schedule, respecting the precedence constraints and minimizing the total completion time C j . This criterion is termed as mean flow time(noted as MFT in the following text). When the number of machines m is not fixed, this problem is noted as P |p j = 1, in-tree| C j in the 3-field notation of [START_REF] Ronald L Graham | Optimization and approximation in deterministic sequencing and scheduling: a surveyhaha[END_REF].

This problem was among the minimal open problems according to Sigrid [START_REF] Brucker | Complexity results for scheduling problems[END_REF]. To the best of our knowledge, it is still open [START_REF] Prot | A survey on how the structure of precedence constraints may change the complexity class of scheduling problems[END_REF]. Herein, we aim at providing a proof of its N P-completeness.

Organization of this paper is as hereunder: in the subsequent section, we present the state of the art. In section 3, we prove that the problem P |p j = 1, in-tree| C j is N P-complete. In subsection 3.1, we present how we reduce the 3-Partition problem to the decision version of the scheduling problem. We prove the N P-completeness in detail in subsection 3.2 and 3.3. Finally, section 4 provides our conclusion.

State of the art

A study involving in-tree precedence constraints in a parallel machine scheduling problem has been performed: Hu (1961) focus on the makespan and suggests that the problem is polynomially solvable using the HLF(highest level first) strategy. Nevertheless, [START_REF] Garey | Scheduling opposing forests[END_REF] prove that the problem is N P-complete when the precedence graph is an opposing forest, which corresponds to a set of in-trees and out-trees.

HLF can also be used to minimize the MFT if the precedence graph is an out-tree. Furthermore, a polynomial algorithm is put forward by [START_REF] Brucker | A polynomial algorithm for P |p j = 1, r j , outtree| C j[END_REF] for the solution of this problem in case of allowance of preemption.

Nevertheless, [START_REF] Huo | Minimizing mean flow time for UET tasks[END_REF] show HLF cannot optimally solve the in-tree version. In addition, [START_REF] Baptiste | Ten notes on equal-processingtime scheduling[END_REF] proved that it can be solved in O(n m ) time, i.e. if the number of machines m is a fixed parameter, the problem is polynomially solvable.

Finally, [START_REF] Garey | Scheduling opposing forests[END_REF] prove that the problem is N P-complete for a profile scheduling. Accordingly, the number of available machines varies along the time rather than parallel machines. However, the complexity of the given parallel scheduling problem P |p j = 1, in-tree| C j is still open.

Proof of N P-completeness

In a bid to prove N P-completeness of P |p j = 1, in-tree| C j , we employed a reduction from a known N P-complete problem: 3-Partition. An instance π 1 of the 3-Partition problem (noted Π 1 ) is defined as hereunder [START_REF] Michael | Computers and intractability[END_REF]: Definition 1. Let α be a set of 3q elements and B an integer, we have ∀a ∈ α, ∃w a ∈ Z + , B 4 < w a < B 2 and a∈α w a = qB. The question is: is there a partition of α, which is α 1 , α 2 , . . . , α q , such that ∀i ≤ q, a∈αi w a = B ?

As Π 1 is N P-complete in the strong sense, the theorem provided hereunder allows us to employ a pseudopolynomial transformation (Garey and Johnson, 2002, p. 101):

Theorem 1. If Π is N P-complete in the strong sense, Π ∈N P and there exists a pseudo-polynomial transformation from Π to Π , then Π is N P-complete in the strong sense.

Firstly, we define Π 2 as P |p j = 1, in -tree| C j ≤ M F T * , which represents the decision version of our scheduling problem. The decision question deals with whether there exists a feasible schedule such that C j ≤ M F T * or not, where M F T * is a given integer. In the subsequent part, we demonstrate how to reduce any instance π 1 of Π 1 to an instance π 2 of Π 2 .

Transformation

In order to begin our transformation, we define some constants:

-B = q 5 B -∀a ∈ α, w a = q 5 w a -K = B 5 m = 3q + 1 + B, which will be used as the number of machines

-M F T * = q i=1 i(3(q + 1 -i) + 1) + K 3 i=q+1 i + B q i=1 i + a∈α Kwa i=1 i + KB q i=1 i -∀k ≤ q + 1, m k = B + 3k -3
Notice that an optimal partition for π 1 , such that a∈αi w a = B, is totally equivalent to a partition such that a∈αi w a = B , since the problem does not change when multiplying w a by an integer q 5 . Thereafter, we create the following jobs:

u 1 a u 2 a u 3 a u wa a v 1 a v 2 a v 3 a v 4 a v 5 a v Kwa a
Fig. 1: in-tree of u-jobs and v-jobs -For each a ∈ α, we build two sets of jobs, which include u-jobs and v-jobs: u 1 a , u 2 a , . . . , u wa a , and

v 1 a , v 2 a , . . . , v Kwa a . The u-jobs precede directly v 1 a , i.e. ∀i ≤ w a , u i a ≺ v 1 a . Subsequent to that, we set v i a ≺ v i+1
a , ∀i < Kw a , and they form a chain. The corresponding tree of the u-jobs and v-jobs can be observed in the Figure 1. They form 3q in-trees -In addition, we consider a set of d-jobs, which are defined as hereunder:

d 1 1 , . . . , d m-m1 1 d 1 2 , . . . , d m-m2 2 . . . d 1 q , d 2 q , d 3 q , d 4 q d q+1 d q+2 . . . d K 3
Those d-jobs also design a tree. In that tree, we fix for any k ≤ q + 1,

d 1 k-1 ≺ d 1 k , d 2 k-1 ≺ d 1 k and d 3 k-1 ≺ d 1 k . ∀i ∈ {4, 5, . . . , m -m k-1 }, then, we set: d i k-1 ≺ d i-3 k . Thereafter, for ∀k ∈ {q + 2, q + 3, . . . , K 3 }, d k-1 ≺ d k , accordingly forming a chain.
The corresponding tree is presented in the Figure 2. Thus, we transformed π 1 to an instance π 2 of the scheduling problem Π 2 . The construction is made in polynomial time of B and q (pseudo-polynomial). In accordance to the Theorem 1, we can establish the N Pcompleteness of Π 2 through the provision of theorem provided as hereunder:

d 1 1 d 2 1 d 3 1 d 4 1 d 5 1 d 6 1 d 7 1 d 8 1 d 9 1 d 10 1 d m-m1 1 d 1 2 d 2 2 d 3 2 d 4 2 d 5 2 d 6 2 d 7 2 d m-m2 2 d 1 3 d 2 2 d 3 2 d 4 2 d m-m3 3 d 4 q d 1 q d q+1 d q+2 d q+3 d q+4 d q+5 d K 3 d 3 q d 2 q ... ... ... ...

Fig. 2: in-tree of d-jobs

Theorem 2. There exists a partition for π 1 ∈ Π 1 , if and only if there exists a feasible schedule for π 2 such that We firstly show that if there exists a partition α 1 , α 2 , . . . , α q for π 1 , such that a∈αi w a = B , then, we can build a feasible schedule σ * for π 2 and C j ≤ M F T * . w.l.o.g: moreover, we assume that, for each subset α k in the partition for k = 1, . . . , q, the three elements in it are a k 1 , a k 2 and a k 3 . We make use of notation M F T * d (resp. M F T * v and M F T * u ) for representing the mean flow time of d-jobs(resp. v-jobs and u-jobs) in σ * . In the subsequent part, we show how jobs are scheduled in σ * together with calculating the MFT.

C j ≤ M F T * . 3.2 Proof of Theorem 2: ⇒ 0 1 2 3 4 q -1 q q + 1 t 5 6 dummy jobs v-jobs u-jobs
Firstly, d-jobs are executed as soon as they are liberated. Formally, ∀i ≤ q, ∀j ≤ m -m i , C d j i = i; ∀i ≤ K 3 , C di = i. This enforces a profile, wherein, the k th timeslot for k ≤ q + 1 has m k available machines for both u-jobs and v-jobs.

Accordingly, we consider MFT for the d-jobs:

M F T * d = q i=1 m-m i j=1 C d j i + K 3 i=q+1 C di (1) as ∀i ≤ K 3 , ∀j ≤ m -m i , C d j i = i, the first part is equivalent to: q i=1 m-mi j=1 C d j i = q i=1 m-mi j=1 i = q i=1 i(m -m i )
and as m -m i = 3(q + 1 -i) + 1, we have:

q i=1 m-mi j=1 C d j i = q i=1 i(3(q + 1 -i) + 1)
Then, for i > q, we also have C di = i, so the second part of (1) is:

K 3 i=q+1 C di = K 3 i=q+1 i
Accordingly, we can obtain the MFT of all d-jobs:

M F T * d = q i=1 i(3(q + 1 -i) + 1) + K 3 i=q+1 i (2)
In the second step, all the u-jobs corresponding to the subset α k are executed in the same time-slot, which is the k th time-slot. They need exactly B machines by definition of 3-Partition. Formally, ∀k ≤ q, ∀a ∈ α k , ∀i ∈ {1, 2, 3}, and ∀j ≤ w a k i , C u j a k i = k. Accordingly, there are exactly B u-jobs executed at each time slot, i.e. |{C u = t}|= B, ∀t ≤ q. Thus, we are able to deduce that all the u-jobs are accomplished before t = q, and we have:

M F T * u = q t=1 t|{u|C u = t}|= B q i=1 i (3)
Thereafter, every v-job is executed without any delay in accordance with the precedence constraints. Formally, ∀a ∈ α, C v 1 a = C u 1 a + 1 and ∀i ∈ {2, 3, . . . , Kw a },

C v i a = C v i-1 a + 1. So, for each v-chain corresponding to the element a ∈ α, we have C v i a = C v i-1 a + 1, ∀i ∈ {2, 3, . . . , Kw a }, that is C v i a = C v 1 a + i -1. Let u a be any u-job preceding v 1 a , as C v 1 a = C ua + 1, we have C v i a = C ua + i
. Now, we calculate the MFT for this given chain:

Kwa i=1 C v i a = Kwa i=1 C ua + Kwa i=1 i = Kw a C ua + Kwa i=1 i
Then, we deduce the MFT for all these v-chains, i.e. all the v-jobs:

M F T * v = ∀a∈α Kwa i=1 C v i a = K ∀a∈α w a C ua + ∀a∈α Kwa i=1 i
We are aware of the fact that, for each subset α k , we have C u a i j = i, where j ∈ {1, 2, 3}, and i ≤ q. So, we are able to deduce that the first part is equivalent to:

K ∀a∈α w a j≤wa C u j a = K q i=1 ∀a∈αi w a j≤wa C u j a = K q i=1 i ∀a∈αi w a = K q i=1 iB So, M F T * v = K q i=1 iB + ∀a∈α Kwa i=1 i (4)
So, in k th time slot of σ * , k < q+1, there are exactly m-m k d-jobs, B v-jobs and 3k -3 v-jobs. By definition of m k , there are exactly m jobs in progress. After t = q, there are only one d-job and at most 3q v-jobs executed in parallel. Thus, there are no more than m jobs. So, this schedule respects the availability of machines.

Then, we deduce that the total MFT of σ * , which is a feasible schedule for π 2 , is:

C j = M F T * v + M F T * u + M F T * d = M F T *

Proof of Theorem 2: ⇐

In this section, we prove that if π 2 has a feasible schedule, then π 1 has a partition. A schedule without unforced idle time is termed as non-delay schedule by [START_REF] Weiss | Scheduling: Theory, algorithms, and systems[END_REF]. The authors also prove that there always exists an optimal non-delay schedule. Let σ be a feasible non-delay schedule. We prove in the following part that σ takes the same shape as σ * and followed by building a partition from σ.

For that purpose, we at first, define the MFT of u-jobs, v-jobs and d-jobs in σ as M F T u , M F T v and M F T d . Subsequent to that, define the following gaps:

∆ u = M F T u -M F T * u ∆ v = M F T v -M F T * v ∆ d = M F T d -M F T * d
In order to respect the M F T * , we understand that:

∆ d + ∆ u + ∆ v ≤ 0 (5)
In addition to this, we would like to throw emphasis on the following observations:

Observation 1. When a chain of n c jobs is rightshifted by one time slot, then its MFT will be increased by n c .

Observation 2. The sizes of different sets of jobs exhibit the following relation:

K 3 Kw a B q 5
First of all, let us consider the d-jobs. We show that d-jobs must have the same profile as in σ * , which is the following lemma: Lemma 1. In σ, the d-jobs are scheduled as soon as they are liberated (as in σ * ), i.e. ∆ d = 0.

Proof. By contradiction: assume that there is at least one delayed d-job. σ is a non-delay schedule; accordingly, this d-jobs can be delayed only by insertion of u-jobs or v-jobs.

Notice that the d-jobs can take no more than m-m 1 machines in parallel because of the precedence constraints. At any moment, there are at least m 1 = B available machines for both the v-jobs and u-jobs. Moreover, in a non-delay schedule, the makespan(C max ) is not larger than the number of jobs. As the total number of v-jobs and u-jobs is KB + qB, the last v-job or u-job finishes before KB + qB. Consequently, the delay must happen before KB + qB.

As the length of the d-chain is K 3 , at least K 3 -KB + qB jobs are right-shifted. Following the Observation 1, we have

∆ d ≥ K 3 -KB + qB Remind that M F T * u K 2 , M F T * v K 2 , so we have ∆ d M F T * u + M F T * v
In accordance with the definition, we have

∆ u > -M F T * u and ∆ v > -M F T * v , then ∆ d + ∆ v + ∆ u > ∆ d -M F T * u -M F T * v 0
This is a contradiction to (5). So, our assumption is false, there is no delay during the execution of djobs.

Since ∆ d = 0, (5) becomes

∆ u + ∆ v ≤ 0 (6)
Here we conclude that the profile of available machines in σ is the same as in σ * to schedule u-jobs and v-jobs. Now, let us deal with u-jobs.

Definition 2. We define the u-jobs in the same component a u-set and define U t as the number of u-jobs in u-sets completed at t. Lemma 2. At time t * ≤ q, the completed u-sets contain at most t * B u-jobs, i.e.

t * t=1 U t ≤ t * B.
Proof. The number of jobs finished at t * ≤ q is at most t * m which is:

t * m = t * B + t * (3q + 1) ≤ t * B + 3q 3
As we are aware of the fact, the number of u-jobs in a u-set is w a , which is integer times q 5 , so, the number of u-jobs in completed u-sets is at most t * B. Now, let us move on to the v-jobs.

Definition 3. We define V t as the number of v-jobs in the chains starting at the time t for t < q in σ, and V q as the number of v-jobs in the chains starting at or after t = q. Lemma 3. In σ, at any time t ≤ q, the total length of the v-chains beginning at t in σ is exactly KB, as in σ * , i.e. V t = KB.

Proof. At first, we prove that

q i=1 iV i ≥ KB q i=1 i.
According to the precedence constraints, we have a relation between V t and U t :

t * t=1 V t ≤ K t * t=1 U t ∀t * < q
So, from Lemma 2, we can deduce:

t * t=1 V t ≤ t * BK ∀t * < q which is: (7) V 1 ≤ KB V 1 + V 2 ≤ 2KB V 1 + V 2 + V 3 ≤ 3KB . . . . . . V 1 + V 2 + V 3 + • • • + V q-1 ≤ (q -1)KB
If we sum them up, we get:

(q -1)V 1 + (q -2)V 2 + • • • + V q-1 ≤ KB q 2 -q 2 (8)
By definition, we understand that q t=1 V t corresponds to the norm of whole set of v-jobs and it is equal to KqB:

q t=1 V t = KqB (9)
Then, by multiplying ( 9) by q, we get:

qV 1 + qV 2 + . . . qV q = KBq 2 (10)
With ( 10)-( 8), we get

V 1 + 2V 2 + • • • + (q -1)V q-1 + qV q ≥ KB q 2 + q 2 which is exactly q i=1 iV i ≥ KB q i=1 i
It requires observation that the two terms are equal if and only if every decomposition of group ( 7) is written with an equality, which indicates that

V 1 = V 2 = • • • = V q = KB. Now, we demonstrate that q i=1 iV i = KB q i=1 i. Sup- pose by absurd that q i=1 iV i > KB q i=1
i. By definition, V i is the total length of some v-chains, which implies that V i is integer times K. So, we have:

q i=1 iV i ≥ KB q i=1 i + K (11)
In the best case of σ, all the v-jobs in one chain are executed without any delay, so, we have:

M F T v ≥ a∈α Kwa i=1 i + q i=1 iV i (12)
Reminding the expression of M F T * v in (4), we can deduce that:

∆ v ≥ q i=1 iV i -KB q i=1 i (13)
By ( 13) and ( 11), we have∆ v ≥ K. Then, by (6), we have ∆ u ≤ -K. This contradicts the fact that ∆ u ≥ -M F T * u -K. Finally, we can conclude that ∆ v = 0 and V t = KB, ∀t ≤ q. This suggests that, at any t, the length of vchains started is exactly KB.

The Lemma above allows us to build a partition from V i by selecting the elements corresponding to its 3 v-chains of KB jobs. Thus, we accomplished the proof for the Theorem 2.

Conclusion

In this paper, we put forward a proof of N P-completeness for P |p j = 1, in-tree| C j . Based on that, we refine the knowledge about the complexity of parallel scheduling problems subjected to precedence constraints in specific graphs.

Nevertheless, some open problems still exist requiring investigation. For example, when preemption is allowed, the complexity of the problem P |p j = 1, in-tree, pmtn| C j is unknown. We conjecture that, with the same reduction as in subsection 3.1, the Theorem 2 still holds even if preemption is allowed, and the problem shall be N Pcomplete as well.

Another possible research avenue involves exploration of how to solve the problem when m is a fixed number. The complexity of the algorithm proposed by [START_REF] Baptiste | Ten notes on equal-processingtime scheduling[END_REF] is O(n m ), but it can be certainly improved.

Fig. 3 :

 3 Fig. 3: Gantt graph of σ *

Acknowledgments

This work was supported by the China Scholarship Council [grant numbers 201404490037].


Shiqi Wang
email: s_wang@insa-toulouse.fr

Lucien Baldas
email: baldas@insa-toulouse.fr

Stéphane Colin
email: colin@insa-toulouse.fr

Stéphane Orieux
email: orieux@insa-toulouse.fr

Nicolas Laurien
email: laurien@insa-toulouse.fr

Azeddine Kourta
email: azeddine.kourta@univ-orleans.fr

Nicolas Mazellier
email: nicolas.mazellier@univ-orleans.fr

Stéphane Loyer
email: stephane.loyer@univ-orleans.fr

A Marie-Curie-ITN within H2020 migrate2017:154498 ACTIVE FLOW CONTROL OF RAMP FLOW BY FLUIDIC OSCILLATORS

Keywords: Separation control, hot wire, PIV, momentum coefficient

ABSTRACT

The study of actuators for active flow control has been in rapid expansion in the last several decades, pursuing different goals such as reducing drag on bluff bodies [START_REF] Aubrun | Separation flow control on a generic ground vehicle using steady microjet arrays[END_REF] , increasing lift of airfoils [START_REF] Melton | Active flow control at low Reynolds numbers on a NACA 0015 airfoil[END_REF][START_REF] Lin | Experimental Study of Flow Separation over NACA63 3 018 Wing with Synthetic Jet Control at Low Reynolds Numbers[END_REF] or enhancing mixing in combustion chambers [START_REF] Fanping | Air Flow Control by Fluidic Diverter for Low NOx Jet Engine Combustion, 1st Flow Control Conference[END_REF][START_REF] Daniel | Active Combustion Control Using a Fluidic Oscillator for Asymmetric Fuel Flow Modulation[END_REF] . Compared to traditional passive control methods or steady blowing method, the active flow control based on periodic fluidic excitations is much more efficient, with a gain of two orders of magnitude in terms of added momentum coefficient, as demonstrated by numerous researches (e.g. Greenblatt and Wygnanski 6 ). These periodic fluidic disturbances can be provided by various kinds of actuators such as ZNMF (Zero Net Mass Flow) actuators, plasma actuators or MEMS (Micro-Electro-Mechanical-Systems) [START_REF] Cattafesta | Actuators for active flow control[END_REF] . Among them, fluidic oscillators can emit oscillating jets in a large operating frequency and velocity range when supplied with a pressurized fluid without requiring any moving part, since their oscillations are totally self-induced and self-sustained and only depend on the internal flow dynamics, which is a great advantage in terms of reliability and robustness [START_REF] Ciro | An Experimental and Numerical Investigation on Fluidic Oscillators For Flow Control[END_REF][START_REF] Gregory | Characterization of the Microfluidic Oscillator[END_REF][START_REF] Bobusch | Experimental study of the internal flow structures inside a fluidic oscillator[END_REF] .

A typical fluidic oscillator is basically composed of an inlet nozzle N, two feedback loops F1 and F2 and two outlets O1 and O2, as shown in Figure 1a. Its behavior is based on the Coanda effect: the jet issuing from nozzle N attaches one of the two walls W1 or W2. The attachment either to wall W1 or wall W2 depends on the initial conditions or is the result of specific actions on the jet. If there was no feedback loop and if the outlet sections were large, the attachment to wall W1 or wall W2 would be stable and the flow would exit through the corresponding outlet, O1 or O2, respectively. With feedback loops, when the jet is attached to wall W1, part of the flow fills in the feedback loop F1 and a pressure increase in the left side of the device is observed, due to the hydraulic restriction at outlet O1. This pressure increase forces the jet to switch toward the right side. Following the jet switching, the same phenomenon develops in the right side of the oscillator and results in a self-sustained oscillating behavior, with a pulsed flow alternatively exiting outlets O1 and O2. In the current design, the feedback loops are plastic tubes as shown in Figure 1b, connected perpendicularly to the base plate. The channels of the oscillator's central part are milled in the base plate in a depth of 370 µm while the outlet slot is milled in the cover plate with an area of 0.5 mm 2 . The throat section of inlet nozzle has a width of about 200 µm. The outlet jet has the same direction as the inlet air.

After examining by hot wire the frequency and velocity response of an isolated oscillator to inlet pressure, an array of 12 identical fluidic oscillators is integrated in a ramp with a 25° slant angle by assembling the cover plate, the base plate and the ramp together as shown in Figure 1c. The centers of outlet  corresponding author The 2D mean velocity fields are measured in 4 cases, which are the baseline flow without control and the flow controlled with the fluidic oscillators supplied at 3 different pressure levels. In all cases, the main flow velocity is 20 m/s, corresponding to a Reynolds number Re =3.8×10 4 based on the ramp height (h=30 mm) as a characteristic length. In the controlled cases, inlet pressures are P i = 0.2, 0.25 and 0.3 MPa respectively, and their corresponding control parameters are presented in Table 1.

f is the measured frequency of the synchronized fluidic oscillator array, while b U and rms b U are the mean and the root mean square values of the blowing velocity in each outlet slot. Assuming that the jet sinusoidally oscillates, the momentum coefficient defined as a ratio of the momentum added to that in the free-stream can be calculated as 6 :

2 2 22 () 11 22 rms b b b b m U N U U A C U L U L        (1)
where b m and b U are the instantaneous blowing mass flow rate and velocity of the outlet jets, N is the total number of outlet slots which is N = 24, A is the outlet slot area, ω is the width of the ramp, L is a characteristic length, here chosen as the length between the ramp slant edge and the separation reattachment point in the baseline flow without control and is the free stream flow velocity. The mean velocity fields in the streamwise direction obtained in both baseline case and controlled cases are shown in Figure 2. In the baseline case ( 0 C   ), the red line 0

P i (
x U  highlights the upper limit of the A Marie-Curie-ITN within H2020

  of the array are at a distance of 200 µm upstream from the ramp slant edge. The efficiency of this fluidic oscillator array on the ramp separation flow is measured in the S2 wind tunnel, in PRISME laboratory, Orléans. The 2D mean velocity fields are captured by using a 2D PIV system. PIV images are focused on the centerline plane of the ramp.

Figure 1 :

 1 Figure 1: a) sketch of a typical fluidic oscillator; b) photo of a real fluidic oscillator during test; c) sketch of the designed ramp

.

  where the flow is reversed. This recirculation region terminates at approximately / 4flow separation has been totally eliminated and no reversed flow can be observed. When  C continues increasing to 0.27%, the situation becomes even better as the minimumx U velocity becomes higher. In case of The obtained results show that this fluidic oscillator array is very promising for active separation flow control considering the low momentum coefficient needed to totally eliminate the separation.

Figure 2 :

 2 Figure 2: average U x field contours under various conditions

Table 1 :

 1 control parameters of three controlled flow cases

		MPa)	0.2	0.25	0.3
	f (Hz)	716	660	660
	b U (m/s)	44	57	70
	U	rms b	(m/s)	6.4	7.7	9.2
	C 		0.16%	0.27%	0.41%

Acknowledgements

The authors acknowledge the scholarship offered by China Scholarship Council (CSC).