Edouard Grave
email: edouard.grave@inria.fr

Guillaume Obozinski
email: guillaume.obozinski@inria.fr

Francis Bach
email: francis.bach@inria.fr

Trace Lasso: a trace norm regularization for correlated designs

Using the ℓ1-norm to regularize the estimation of the parameter vector of a linear model leads to an unstable estimator when covariates are highly correlated. In this paper, we introduce a new penalty function which takes into account the correlation of the design matrix to stabilize the estimation. This norm, called the trace Lasso, uses the trace norm, which is a convex surrogate of the rank, of the selected covariates as the criterion of model complexity. We analyze the properties of our norm, describe an optimization algorithm based on reweighted least-squares, and illustrate the behavior of this norm on synthetic data, showing that it is more adapted to strong correlations than competing methods such as the elastic net.

Introduction

The concept of parsimony is central in many scientific domains. In the context of statistics, signal processing or machine learning, it takes the form of variable or feature selection problems, and is commonly used in two situations: first, to make the model or the prediction more interpretable or cheaper to use, i.e., even if the underlying problem does not admit sparse solutions, one looks for the best sparse approximation. Second, sparsity can also be used given prior knowledge that the model should be sparse. Many methods have been designed to learn sparse models, namely methods based on combinatorial optimization [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF][START_REF] Zhang | Adaptive forward-backward greedy algorithm for sparse learning with linear models[END_REF], Bayesian inference [START_REF] Seeger | Bayesian inference and optimal design for the sparse linear model[END_REF] or convex optimization [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF][START_REF] Chen | Atomic decomposition by basis pursuit[END_REF].

In this paper, we focus on the regularization by sparsity-inducing norms. The simplest example of such norms is the ℓ 1 -norm, leading to the Lasso, when used within a least-squares framework. In recent years, a large body of work has shown that the Lasso was performing optimally in highdimensional low-correlation settings, both in terms of prediction [START_REF] Bickel | Simultaneous analysis of Lasso and Dantzig selector[END_REF], estimation of parameters or estimation of supports [START_REF] Zhao | On model selection consistency of Lasso[END_REF][START_REF] Wainwright | Sharp thresholds for high-dimensional and noisy sparsity recovery using ℓ 1 -constrained quadratic programming (Lasso). Information Theory[END_REF]. However, most data exhibit strong correlations, with various correlation structures, such as clusters (i.e., close to block-diagonal covariance matrices) or sparse graphs, such as for example problems involving sequences (in which case, the covariance matrix is close to a Toeplitz matrix [START_REF] Golub | Matrix computations[END_REF]). In these situations, the Lasso is known to have stability problems: although its predictive performance is not disastrous, the selected predictor may vary a lot (typically, given two correlated variables, the Lasso will only select one of the two, at random).

Several remedies have been proposed to this instability. First, the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF] adds a strongly convex penalty term (the squared ℓ 2 -norm) that will stabilize selection (typically, given two correlated variables, the elastic net will select the two variables). However, it is blind to the exact correlation structure, and while strong convexity is required for some variables, it is not for other variables. Another solution is to consider the group Lasso, which will divide the predictors into groups and penalize the sum of the ℓ 2 -norm of these groups [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. This is known to accomodate strong correlations within groups [START_REF] Bach | Consistency of the group Lasso and multiple kernel learning[END_REF]; however it requires to know the group in advance, which is not always possible. A third line of research has focused on sampling-based techniques [START_REF] Bach | Bolasso: model consistent Lasso estimation through the bootstrap[END_REF][START_REF] Liu | Stability approach to regularization selection (stars) for high dimensional graphical models[END_REF][START_REF] Meinshausen | Stability selection[END_REF]]. An ideal regularizer should thus be adapted to the design (like the group Lasso), but without requiring human intervention (like the elastic net); it should thus add strong convexity only where needed, and not modifying variables where things behave correctly. In this paper, we propose a new norm towards this end.

More precisely we make the following contributions:

• We propose in Section 2 a new norm based on the trace norm (a.k.a. nuclear norm) that interpolates between the ℓ 1 -norm and the ℓ 2 -norm depending on correlations.

• We show that there is a unique minimum when penalizing with this norm in Section 2.2.

• We provide optimization algorithms based on reweighted least-squares in Section 3.

• We study the second-order expansion around independence and relate to existing work on including correlations in Section 4.

• We perform synthetic experiments in Section 5, where we show that the trace Lasso outperforms existing norms in strong-correlation regimes.

Notations. Let M ∈ R n×p . The columns of M are noted using superscript, i.e., M (i) denotes the i-th column, while the rows are noted using subscript, i.e., M i denotes the i-th row. For M ∈ R p×p , diag(M) ∈ R p is the diagonal of the matrix M, while for u ∈ R p , Diag(u) ∈ R p×p is the diagonal matrix whose diagonal elements are the u i . Let S be a subset of {1, ..., p}, then u S is the vector u restricted to the support S, with 0 outside the support S. We denote by S p the set of symmetric matrices of size p. We will use various matrix norms, here are the notations we use:

• M * is the trace norm, i.e., the sum of the singular values of the matrix M,

• M op is the operator norm, i.e., the maximum singular value of the matrix M,

• M F is the Frobenius norm, i.e., the ℓ 2 -norm of the singular values, which is also equal to tr(M ⊤ M),

• M 2,1 is the sum of the ℓ 2 -norm of the columns of

M: M 2,1 = p i=1 M (i) 2 .
2 Definition and properties of the trace Lasso

We consider the problem of predicting y ∈ R, given a vector x ∈ R p , assuming a linear model

y = w ⊤ x + ε,
where ε is (Gaussian) noise with mean 0 and variance σ 2 . Given a training set X = (x 1 , ..., x n ) ⊤ ∈ R n×p and y = (y 1 , ..., y n ) ⊤ ∈ R n , a widely used method to estimate the parameter vector w is the penalized empirical risk minimization

ŵ ∈ argmin w 1 n n i=1 ℓ(y i , w ⊤ x i ) + λf (w), ( 1
)
where ℓ is a loss function used to measure the error we make by predicting w ⊤ x i instead of y i , while f is a regularization term used to penalize complex models. This second term helps avoiding overfitting, especially in the case where we have many more parameters than observation, i.e., n ≪ p.

Related work

We will now present some classical penalty functions for linear models which are widely used in the machine learning and statistics community. The first one, known as Tikhonov regularization [START_REF] Tikhonov | Solution of incorrectly formulated problems and the regularization method[END_REF] or ridge regression [START_REF] Hoerl | Ridge regression: Biased estimation for nonorthogonal problems[END_REF], is the squared ℓ 2 -norm. When used with the square loss, estimating the parameter vector w is done by solving a linear system. One of the main drawbacks of this penalty function is the fact that it does not perform variable selection and thus does not behave well in sparse high-dimensional settings. Hence, it is natural to penalize linear models by the number of variables used by the model. Unfortunately, this criterion, sometimes denoted by • 0 (ℓ 0 -penalty), is not convex and solving the problem in Eq. ( 1) is generally NP-hard [START_REF] Davis | Adaptive greedy approximations[END_REF]. Thus, a convex relaxation for this problem was introduced, replacing the size of the selected subset by the ℓ 1 -norm of w. This estimator is known as the Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the Lasso[END_REF] in the statistics community and basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] in signal processing. It was later shown that under some assumptions, the two problems were in fact equivalent (see for example [START_REF] Candes | Decoding by linear programming[END_REF] and references therein).

When two predictors are highly correlated, the Lasso has a very unstable behavior: it may only select the variable that is the most correlated with the residual. On the other hand, the Tikhonov regularization tends to shrink coefficients of correlated variables together, leading to a very stable behavior. In order to get the best of both worlds, stability and variable selection, Zou and Hastie introduced the elastic net [START_REF] Zou | Regularization and variable selection via the elastic net[END_REF], which is the sum of the ℓ 1 -norm and squared ℓ 2 -norm. Unfortunately, this estimator needs two regularization parameters and is not adaptive to the precise correlation structure of the data. Some authors also proposed to use pairwise correlations between predictors to interpolate more adaptively between the ℓ 1 -norm and squared ℓ 2 -norm, by introducing the pairwise elastic net [START_REF] Lorbert | Exploiting covariate similarity in sparse regression via the pairwise elastic net[END_REF] (see comparisons with our approach in Section 5).

Finally, when one has more knowledge about the data, for example clusters of variables that should be selected together, one can use the group Lasso [START_REF] Yuan | Model selection and estimation in regression with grouped variables[END_REF]. Given a partition (S i ) of the set of variables, it is defined as the sum of the ℓ 2 -norms of the restricted vectors w Si :

w GL = k i=1 w Si 2 .
The effect of this penalty function is to introduce sparsity at the group level: variables in a group are selected altogether. One of the main drawback of this method, which is also sometimes one of its quality, is the fact that one needs to know the partition of the variables, and so one needs to have a good knowledge of the data.

The ridge, the Lasso and the trace Lasso

In this section, we show that Tikhonov regularization and the Lasso penalty can be viewed as norms of the matrix X Diag(w). We then introduce a new norm involving this matrix.

The solution of empirical risk minimization penalized by the ℓ 1 -norm or ℓ 2 -norm is not equivariant by rescaling of the predictors X (i) , so it is common to normalize the predictors. When normalizing the predictors X (i) , and penalizing by Tikhonov regularization or by the Lasso, people are implicitly using a regularization term that depends on the data or design matrix X. In fact, there is an equivalence between normalizing the predictors and not normalizing them, using the two following reweighted ℓ 2 and ℓ 1 -norms instead of the Tikhonov regularization and the Lasso:

w 2 2 = p i=1 X (i) 2 2 w 2 i and w 1 = p i=1 X (i) 2 |w i |. (2)
These two norms can be expressed using the matrix X Diag(w):

w 2 = X Diag(w) F and w 1 = X Diag(w) 2,1 ,
and a natural question arises: are there other relevant choices of functions or matrix norms? A classical measure of the complexity of a model is the number of predictors used by this model, which is equal to the size of the support of w. This penalty being non-convex, people use its convex relaxation, which is the ℓ 1 -norm, leading to the Lasso.

Here, we propose a different measure of complexity which can be shown to be more adapted in model selection settings [START_REF] Hastie | The elements of statistical learning[END_REF]: the dimension of the subspace spanned by the selected predictors. This is equal to the rank of the selected predictors, or also to the rank of the matrix X Diag(w). As for the size of the support, this function is non-convex, and we propose to replace it by a convex surrogate, the trace norm, leading to the following penalty that we call "trace Lasso":

Ω(w) = X Diag(w) * .
The trace Lasso has some interesting properties: if all the predictors are orthogonal, then, it is equal to the ℓ 1 -norm. Indeed, we have the decomposition:

X Diag(w) = p i=1 X (i) 2 w i X (i) X (i) 2 e ⊤ i ,
where e i are the vectors of the canonical basis. Since the predictors are orthogonal and the e i are orthogonal too, this gives the singular value decomposition of X Diag(w) and we get

X Diag(w) * = p i=1 X (i) 2 |w i | = X Diag(w) 2,1 .
On the other hand, if all the predictors are equal to X (1) , then

X Diag(w) = X (1) w ⊤ ,
and we get X Diag(w) * = X (1) 2 w 2 = X Diag(w) F , which is equivalent to the Tikhonov regularization. Thus when two predictors are strongly correlated, our norm will behave like the Tikhonov regularization, while for almost uncorrelated predictors, it will behave like the Lasso. Always having a unique minimum is an important property for a statistical estimator, as it is a first step towards stability. The trace Lasso, by adding strong convexity exactly in the direction of highly correlated covariates, always has a unique minimum, and is much more stable than the Lasso.

Proposition 1. If the loss function ℓ is strongly convex with respect to its second argument, then the solution of the empirical risk minimization penalized by the trace Lasso, i.e., Eq. ( 1), is unique.

The technical proof of this proposition is given in appendix B, and consists of showing that in the flat directions of the loss function, the trace Lasso is strongly convex.

A new family of penalty functions

In this section, we introduce a new family of penalties, inspired by the trace Lasso, allowing us to write the ℓ 1 -norm, the ℓ 2 -norm and the newly introduced trace Lasso as special cases. In fact, we note that Diag(w) * = w 1 and p -1/2 1 ⊤ Diag(w) * = w ⊤ * = w 2 . In other words, we can express the ℓ 1 and ℓ 2 -norms of w using the trace norm of a given matrix times the matrix Diag(w). A natural question to ask is: what happens when using a matrix P other than the identity or the line vector p -1/2 1 ⊤ , and what are good choices of such matrices? Therefore, we introduce the following family of penalty functions: Definition 1. Let P ∈ R k×p , all of its columns having unit norm. We introduce the norm Ω P as

Ω P (w) = P Diag(w) * .
Proof. The positive homogeneity and triangle inequality are direct consequences of the linearity of w → P Diag(w) and the fact that • * is a norm. Since all the columns of P are not equal to zero, we have P Diag(w) = 0 ⇔ w = 0, and so, Ω P separates points and is a norm. As stated before, the ℓ 1 and ℓ 2 -norms are special cases of the family of norms we just introduced. Another important penalty that can be expressed as a special case is the group Lasso, with nonoverlapping groups. Given a partition (S j ) of the set {1, ..., p}, the group Lasso is defined by

w GL = Sj w Sj 2 .
We define the matrix P GL by

P GL ij = 1/ |S k | if i and j are in the same group S k , 0 otherwise.
Then,

P GL Diag(w) = Sj 1 Sj |S j | w ⊤ Sj . (3)
Using the fact that (S j ) is a partition of {1, ..., p}, the vectors 1 Sj are orthogonal and so are the vectors w Sj . Hence, after normalizing the vectors, Eq. ( 3) gives a singular value decomposition of P GL Diag(w) and so the group Lasso penalty can be expressed as a special case of our family of norms:

P GL Diag(w) * = Sj w Sj 2 = w GL .
In the following proposition, we show that our norm only depends on the value of P ⊤ P. This is an important property for the trace Lasso, where P = X, since it underlies the fact that this penalty only depends on the correlation matrix X ⊤ X of the covariates.

Proposition 2. Let P ∈ R k×p , all of its columns having unit norm. We have

Ω P (w) = (P ⊤ P) 1/2 Diag(w) * .
We plot the unit ball of our norm for the following value of P ⊤ P (see figure (

    1 1 0 1 1 0 0 0 1  
We can lower bound and upper bound our norms by the ℓ 2 -norm and ℓ 1 -norm respectively. This shows that, as for the elastic net, our norms interpolate between the ℓ 1 -norm and the ℓ 2norm. But the main difference between the elastic net and our norms is the fact that our norms are adaptive, and require a single regularization parameter to tune. In particular for the trace Lasso, when two covariates are strongly correlated, it will be close to the ℓ 2 -norm, while when two covariates are almost uncorrelated, it will behave like the ℓ 1 -norm. This is a behavior close to the one of the pairwise elastic net [START_REF] Lorbert | Exploiting covariate similarity in sparse regression via the pairwise elastic net[END_REF].

Proposition 3. Let P ∈ R k×p , all of its columns having unit norm. We have

w 2 ≤ Ω P (w) ≤ w 1 .

Dual norm

The dual norm is an important quantity for both optimization and theoretical analysis of the estimator. Unfortunately, we are not able in general to obtain a closed form expression of the dual norm for the family of norms we just introduced. However we can obtain a bound, which is exact for some special cases: Proposition 4. The dual norm, defined by Ω * P (u) = max

Ω P (v)≤1
u ⊤ v, can be bounded by:

Ω * P (u) ≤ P Diag(u) op .
Proof. Using the fact that diag(P ⊤ P) = 1, we have

u ⊤ v = tr Diag(u)P ⊤ P Diag(v) ≤ P Diag(u) op P Diag(v) * ,
where the inequality comes from the fact that the operator norm • op is the dual norm of the trace norm. The definition of the dual norm then gives the result.

As a corollary, we can bound the dual norm by a constant times the ℓ ∞ -norm:

Ω * P (u) ≤ P Diag(u) op ≤ P op Diag(u) op = P op u ∞ .
Using proposition (3), we also have the inequality Ω * P (u) ≥ u ∞ .

Optimization algorithm

In this section, we introduce an algorithm to estimate the parameter vector w when the loss function is equal to the square loss: ℓ(y, w ⊤ x) = 1 2 (yw ⊤ x) 2 and the penalty is the trace Lasso. It is straightforward to extend this algorithm to the family of norms indexed by P. The problem we consider is min

w 1 2 y -Xw 2 2 + λ X Diag(w) * .
We could optimize this cost function by subgradient descent, but this is quite inefficient: computing the subgradient of the trace Lasso is expensive and the rate of convergence of subgradient descent is quite slow. Instead, we consider an iteratively reweighted least-squares method. First, we need to introduce a well-known variational formulation for the trace norm [START_REF] Argyriou | Multi-task feature learning[END_REF]:

Proposition 5. Let M ∈ R n×p .
The trace norm of M is equal to:

M * = 1 2 inf S 0 tr M ⊤ S -1 M + tr (S) ,
and the infimum is attained for S = MM ⊤ 1/2 .

Using this proposition, we can reformulate the previous optimization problem as min

w inf S 0 1 2 y -Xw 2 2 + λ 2 w ⊤ Diag diag(X ⊤ S -1 X) w + λ 2 tr(S).
This problem is jointly convex in (w, S) [START_REF] Boyd | Convex optimization[END_REF]. In order to optimize this objective function by alternating the minimization over w and S, we need to add a term λµi 2 tr(S -1 ). Otherwise, the infimum over S could be attained at a non invertible S, leading to a non convergent algorithm. The infimum over S is then attained for S = X Diag(w)

2 X ⊤ + µ i I 1/2 .
Optimizing over w is a least-squares problem penalized by a reweighted ℓ 2 -norm equal to w ⊤ Dw, where D = Diag diag(X ⊤ S -1 X) . It is equivalent to solving the linear system

(X ⊤ X + λD)w = X ⊤ y.
This can be done efficiently by using a conjugate gradient method. Since the cost of multiplying (X ⊤ X + λD) by a vector is O(np), solving the system has a complexity of O(knp), where k ≤ p is the number of iterations needed to converge. Using warm restarts, k can be much smaller than p, since the linear system we are solving does not change a lot from an iteration to another. Below we summarize the algorithm:

Iterative algorithm for estimating w Input: the design matrix X, the initial guess w 0 , number of iteration N , sequence µ i . For i = 1...N :

• Compute the eigenvalue decomposition U Diag(s k )U ⊤ of X Diag(w i-1 ) 2 X ⊤ . • Set D = Diag(diag(X ⊤ S -1 X)), where S -1 = U Diag(1/ √ s k + µ i )U ⊤ .
• Set w i by solving the system (X ⊤ X + λD)w = X ⊤ y.

For the sequence µ i , we use a decreasing sequence converging to ten times the machine precision.

Choice of λ

We now give a method to choose the regularization path. In fact, we know that the vector 0 is solution if and only if λ ≥ Ω * (X ⊤ y) [START_REF] Bach | Convex optimization with sparsityinducing norms[END_REF]. Thus, we need to start the path at λ = Ω * (X ⊤ y), corresponding to the empty solution 0, and then decrease λ. Using the inequalities on the dual norm we obtained in the previous section, we get

X ⊤ y ∞ ≤ Ω * (X ⊤ y) ≤ X op X ⊤ y ∞ .
Therefore, starting the path at λ = X op X ⊤ y ∞ is a good choice.

Approximation around the Lasso

In this section, we compute the second order approximation of our norm around the special case corresponding to the Lasso. We recall that when P = I ∈ R p×p , our norm is equal to the ℓ 1 -norm. We add a small perturbation ∆ ∈ S p to the identity matrix, and using Prop. 6 of the appendix A, we obtain the following second order approximation:

(I + ∆) Diag(w) * = w 1 + diag(∆) ⊤ |w|+ |wi|>0 |wj |>0 (∆ ji |w i | -∆ ij |w j |) 2 4(|w i | + |w j |) + |wi|=0 |wj |>0 (∆ ij |w j |) 2 2|w j | + o( ∆ 2 ).
We can rewrite this approximation as

(I + ∆) Diag(w) * = w 1 + diag(∆) ⊤ |w| + i,j ∆ 2 ij (|w i | -|w j |) 2 4(|w i | + |w j |) + o( ∆ 2 ),
using a slight abuse of notation, considering that the last term is equal to 0 when w i = w j = 0.

The second order term is quite interesting: it shows that when two covariates are correlated, the effect of the trace Lasso is to shrink the corresponding coefficients toward each other. Another interesting remark is the fact that this term is very similar to pairwise elastic net penalties, which are of the form |w| ⊤ P|w|, where P ij is a decreasing function of ∆ ij .

Experiments

In this section, we perform synthetic experiments to illustrate the behavior of the trace Lasso and other classical penalties when there are highly correlated covariates in the design matrix.

For all experiments, we have p = 1024 covariates and n = 256 observations. The support S of w is equal to {1, ..., k}, where k is the size of the support. For i in the support of w, we have

w i = 2 (b i -1/2),
where each b i is independently drawn from a uniform distribution on [0, 1].

The observations x i are drawn from a multivariate Gaussian with mean 0 and covariance matrix Σ. For the first experiment, Σ is set to the identity, for the second experiment, Σ is block diagonal with blocks equal to 0.2I + 0.811 ⊤ corresponding to clusters of eight variables, finally for the third experiment, we set Σ ij = 0.95 |i-j| , corresponding to a Toeplitz design. For each method, we choose the best λ for the estimation error, which is reported. Overall all methods behave similarly in the noiseless and the noisy settings, hence we only report results for the noisy setting. In all three graphs of Figure 2, we observe behaviors that are typical of Lasso, ridge and elastic net: the Lasso performs very well on sparse models, but its performance is rather poor for denser models, almost as poor as the ridge regression. The elastic net offers the best of both worlds since its two parameters allow it to interpolate adaptively between the Lasso and the ridge. In experiment 1, since the variables are uncorrelated, there is no reason to couple their selection. This suggests that the Lasso should be the most appropriate convex regularization. The trace Lasso approaches the Lasso as n goes to infinity, but the weak coupling induced by empirical correlations is sufficient to slightly decrease its performance compared to that of the Lasso. By contrast, in experiments 2 and 3, the trace Lasso outperforms other methods (including the pairwise elastic net) since variables that should be selected together are indeed correlated. As for the penalized elastic net, since it takes into account the correlations between variables it is not surprising that in experiment 2 and 3 it performs better than methods that do not. We do not have a compelling explanation for its superior performance in experiment 1.

Conclusion

We introduce a new penalty function, the trace Lasso, which takes advantage of the correlation between covariates to add strong convexity exactly in the directions where needed, unlike the elastic net for example, which blindly adds a squared ℓ 2 -norm term in every directions. We show on synthetic data that this adaptive behavior leads to better estimation performance. In the future, we want to show that if a dedicated norm using prior knowledge such as the group Lasso can be used, the trace Lasso will behave similarly and its performance will not degrade too much, providing theoretical guarantees to such adaptivity. Finally, we will seek applications of this estimator in inverse problems such as deblurring, where the design matrix exhibits strong correlation structure.

We note A and B the first two terms of the right hand side of this equation. We have

tr(A) = j,k tr(w j w ⊤ j ∆w k w ⊤ k ) 1 2iπ C λdλ (λ -s j )(λ -s k ) = j tr(w ⊤ j ∆w j ) 1 2iπ C λdλ (λ -s j ) 2 = j tr(w ⊤ j ∆w j ) = j tr(u ⊤ j ∆v j ),
and

tr(B) = j,k,l tr(w j w ⊤ j ∆w k w ⊤ k ∆w l w ⊤ l ) 1 2iπ C λdλ (λ -s j )(λ -s k )(λ -s l ) = j,k tr(w j ∆w k w k ∆w j ) 1 2iπ C λdλ (λ -s j ) 2 (λ -s k )
.

If s j = s k , the integral is nul. Otherwise, we have

λ (λ -s j ) 2 (λ -s k ) = a λ -s j + b λ -s k + c (λ -s j ) 2 , where a = -s k (s k -s j ) 2 , b = s k (s k -s j ) 2 , c = s j s j -s k .
Therefore, if s j and s k are both inside or outside the interior of C, the integral is equal to zero. So

tr(B) = sj >0 s k ≤0 -s k (w ⊤ j ∆w k ) 2 (s j -s k ) 2 + sj ≤0 s k >0 s k (w ⊤ j ∆w k ) 2 (s j -s k ) 2 = sj >0 s k >0 s k (w ⊤ j ∆w -k ) 2 (s j + s k ) 2 + sj >0 s k >0 s k (w ⊤ -j ∆w k ) 2 (s j + s k ) 2 + sj =0 s k >0 (w ⊤ j ∆w k ) 2 s k = sj >0 s k >0 (w ⊤ -j ∆w k ) 2 s j + s k + sj =0 s k >0 (w ⊤ j ∆w k ) 2 s k .
For s j > 0 and s k > 0, we have

w ⊤ -j ∆w k = 1 2 u ⊤ j ∆v k -u ⊤ k ∆v j ,
and for s j = 0 and s k > 0, we have

w ⊤ j ∆w k = 1 2 ±u ⊤ k ∆v j + u ⊤ j ∆v k . So tr(B) = sj >0 s k >0 (u ⊤ j ∆v k -u ⊤ k ∆v j ) 2 4(s j + s k ) + sj =0 s k >0 (u ⊤ k ∆v j ) 2 + (u ⊤ j ∆v k ) 2 2s k .
Now, let C 0 be the circle of center 0 and radius s r /2. We can study the perturbation of the singular values of M equal to zero by computing the trace norm of Π C0 ( M + ∆) -Π C0 ( M). We have

Π C0 ( M + ∆) -Π C0 ( M) = 1 2iπ C0 λ(λI -M) -1 ∆(λI -M) -1 dλ + 1 2iπ C0 λ(λI -M) -1 ∆(λI -M) -1 ∆(λI -M) -1 dλ + 1 2iπ C0 λ(λI -M) -1 ∆(λI -M) -1 ∆(λI -M -∆) -1 dλ.
Then, if we note the first integral C and the second one D, we get

C = j,k w j w ⊤ j ∆w k w ⊤ k 1 2iπ C0 λdλ (λ -s j )(λ -s k ) .
If both s j and s k are outside int(C 0 ), then the integral is equal to zero. If one of them is inside, say s j , then s j = 0 and the integral is equal to

C0 dλ λ -s k Then this integral is non nul if and only if s k is also inside int(C 0 ). Thus C = j,k w j w ⊤ j ∆w k w ⊤ k ✶ sj ∈int(C0) ✶ s k ∈int(C0) = sj =0 s k =0 w j w ⊤ j ∆w k w ⊤ k = W 0 W ⊤ 0 ∆W 0 W ⊤ 0
, where W 0 are the eigenvectors associated to the eigenvalue 0. We have

D = j,k,l w j w ⊤ j ∆w k w ⊤ k ∆w l w ⊤ l 1 2iπ C0 λdλ (λ -s j )(λ -s k )(λ -s l )
.

The integral is not equal to zero if and only if exactly one eigenvalue, say s i , is outside int(C 0 ). The integral is then equal to -1/s i . Thus

D = -W 0 W ⊤ 0 ∆W 0 W ⊤ 0 ∆WS -1 W ⊤ -WS -1 W ⊤ ∆W 0 W ⊤ 0 ∆W 0 W ⊤ 0 -W 0 W ⊤ 0 ∆WS -1 W ⊤ ∆W 0 W ⊤ 0 ,
where S = Diag(-s, s). Finally, putting everything together, we get

Proposition 6. Let M = U Diag(s)V ⊤ ∈ R n×p , the singular value decomposition of M, with U ∈ R n×r , V ∈ R p×r . Let ∆ ∈ R n×p . We have M + ∆ * = M * + Q * + tr(VU ⊤ ∆)+ sj >0 s k >0 (u ⊤ j ∆v k -u ⊤ k ∆v j ) 2 4(s j + s k ) + sj =0 s k >0 (u ⊤ k ∆v 0j ) 2 + (u ⊤ 0j ∆v k ) 2 2s k + o( ∆ 2 ),
where

Q = U ⊤ 0 ∆V 0 -U ⊤ 0 ∆V 0 V ⊤ 0 ∆ ⊤ U Diag(s) -1 -Diag(s) -1 V ⊤ ∆ ⊤ U 0 U ⊤ 0 ∆V 0 -U ⊤ 0 ∆V Diag(s) -1 U ⊤ ∆V 0 .

B Proof of proposition 1

In this section, we prove that if the loss function is strongly convex with respect to its second argument, then the solution of the penalized empirical risk minimization is unique.

Let ŵ ∈ argmin w n i=1 ℓ(y i , w ⊤ x i ) + λ X Diag(w) * . If ŵ is in the nullspace of X, then ŵ = 0 and the minimum is unique. From now on, we suppose that the minima are not in the nullspace of X.

Let u, v ∈ argmin w n i=1 ℓ(y i , w ⊤ x i ) + λ X Diag(w) * and δ = vu. By convexity of the objective function, all the w = u+tδ, for t ∈]0, 1[ are also optimal solutions, and so, we can choose an optimal solution w such that w i = 0 for all i in the support of δ. Because the loss function is strongly convex outside the nullspace of X, δ is in the nullspace of X.

Let X Diag(w) = U Diag(s)V ⊤ be the SVD of X Diag(w). We have the following development around w:

X Diag(w + tδ) * = X Diag(w) * + tr(Diag(tδ)X ⊤ UV ⊤ )+ si>0 sj >0 tr(Diag(tδ)X ⊤ (u i v ⊤ j -u j v ⊤ i )) 2 4(s i + s j ) + si>0 sj =0 tr(Diag(tδ)X ⊤ u i v ⊤ j ) 2 2s i + o(t 2 ).
We note S the support of w. Using the fact that the support of δ is included in S, we have X Diag(tδ) = X Diag(w) Diag(tγ), where γ i = δi wi for i ∈ S and 0 otherwise. Then:

X Diag(w + tδ) * = X Diag(w) * + tγ ⊤ diag(V Diag(s)V ⊤ )+ si>0 sj >0 t 2 tr (s i -s j ) Diag(γ)v i v ⊤ For small t, w + tδ is also a minimum, and therefore, we have:

∀ s i > 0, s j > 0, (s i -s j ) tr Diag(γ

)v i v ⊤ j = 0, (4)
∀ s i > 0, s j = 0, tr Diag(γ)v i v ⊤ j = 0.

This could be summarized as

∀ s i = s j , v ⊤ i (Diag(γ)v j ) = 0. ( 6
)
This means that the eigenspaces of Diag(w)X ⊤ X Diag(w) are stable by the matrix Diag(γ). Therefore, Diag(w)X ⊤ X Diag(w) and Diag(γ) are simultaneously diagonalizable and so, they commute. Therefore: ∀ i, j ∈ S, σ ij γ i = σ ij γ j [START_REF] Zhao | On model selection consistency of Lasso[END_REF] where σ ij = [X ⊤ X] ij . We define a partition (S k ) of S, such that i and j are in the same set S k if there exists a path i = a 1 , ..., a m = j such that σ an,an+1 = 0 for all n ∈ {1, ..., m -1}. Then, using equation ( 7), γ is constant on each S k . δ being in the nullspace of X, we have:

0 = δ ⊤ X ⊤ Xδ (8) = S k S l δ ⊤ S k X ⊤ Xδ S l (9) = S k δ ⊤ S k X ⊤ Xδ S k (10) = S k Xδ S k 2 2 . (11)
So for all S i , Xδ Si = 0. Since a predictor X i is orthogonal to all the predictors belonging to other groups defined by the partition (S k ), we can decompose the norm Ω:

X Diag(w) * = S k X Diag(w S k ) * . (12)
We recall that γ is constant on each S k and so δ S k is colinear to w Si , by definition of γ. If δ Si is not equal to zero, this means that w Si , which is not equal to zero, is in the nullspace of X.

Replacing w Si by 0 will not change the value of the data fitting term but it will strictly decreases the value of the norm Ω. This is a contradiction with the optimality of w. Thus all the δ Si are equal to zero and the minimum is unique.

C Proof of proposition 3

For the first inequality, we have The first equality is the fact that the dual norm of the trace norm is the operator norm and the second inequality uses the fact that all matrices of operator norm smaller than one have columns of ℓ 2 norm smaller than one.

Figure 1 :

 1 Figure 1: Unit balls for various value of P ⊤ P. See the text for the value of P ⊤ P. (Best seen in color).

Figure 2 :

 2 Figure 2: Experiment for uncorrelated variables (Best seen in color. en stands for elastic net, pen stands for pairwise elastic net and trace stands for trace Lasso.)

j 2 4

 2 (s i + s j ) + si>0 sj =0 t 2 tr s i Diag(γ)v i v ⊤ j 2 2s i + o(t 2 ).

w 2 =|M

 2 P Diag(w) F ≤ P Diag(w) * .For the second inequality, we haveP Diag(w) * = max M op ≤1 tr M ⊤ P Diag(w) (i)⊤ P (i) | |w i | ≤ w 1 .

Acknowledgements

This paper was partially supported by the European Research Council (SIERRA Project ERC-239993).

A Perturbation of the trace norm

We follow the technique used in [START_REF] Bach | Consistency of trace norm minimization[END_REF] to obtain an approximation of the trace norm.

A.1 Jordan-Wielandt matrices

Let M ∈ R n×p of rank r. We note s 1 ≥ s 2 ≥ ... ≥ s r > 0, the strictly positive singular values of M and u i , v i the associated left and right singular vectors. We introduce the Jordan-Wielandt matrix

The singular values of M and the eigenvalues of M are related: M has eigenvalues s i and s -i = -s i associated to eigenvectors

The remaining eigenvalues of M are equal to 0 and are associated to eigenvectors of the form

A.2 Cauchy residue formula

Let C be a closed curve that does not go through the eigenvalues of M. We define

We have

A.3 Perturbation analysis

Let ∆ ∈ R n×p be a perturbation matrix such that ∆ op < s r /4, and let C be a closed curve around the r largest eigenvalues of M and M + ∆. We can study the perturbation of the strictly positive singular values of M by computing the trace of Π

λ(λI -M) -1 ∆(λI -M) -1 dλ

λ(λI -M) -1 ∆(λI -M) -1 ∆(λI -M) -1 dλ

λ(λI -M) -1 ∆(λI -M) -1 ∆(λI -M -∆) -1 dλ.


Keywords: Processus de Decision Markovien, Politique à seuils critiques, Systèmes de gestion de stock

Politiques à seuils critiques dans des systèmes de gestion de stock avec pertes et différentes classes de clients Résumé : On s'intéresse ici à un système de gestion de stock manipulant un seul type d'objet et pour lequel on autorise des pertes de ventes. Ce système manipule différentes classes de clients chacune ayant un coût de perte de vente différent. On suppose que les demandes suivent un processus de Poisson et nous nous intéressons au cas à un seul serveur de réassortiment dont la loi est hypoexponentielle. Nous décrivons le Processus de Décision Markovien associé à un tel problème de contrle optimal et montrons quelques propriétés structurelles de son opérateur de programmation dynamique. Nous en déduisons que la politique optimale est du type politique à seuils critiques. Nous discutons ensuite d'extensions lorsque les distributions de temps de réassortiment sont hyperexponentielles. Nous donnons alors quelques simulations numériques qui tendent à montrer l'existence de propriétés structurelles similaires.

Introduction

We consider an inventory system with single item type and a finite stock capacity. Upon a demand, we can decide to either deliver an item or reject the demand. We assume different demand classes with different rejection penalty costs. The objective is to find a sequence of optimal decisions that minimizes the costs. This problem can be easily modeled as a Markov Decision Process (MDP). We refer the reader to Puterman [START_REF] Puterman | Markov Decision Processes Discrete Stochastic Dynamic Programming[END_REF] for general theory on MDPs and to Porteus [START_REF] Porteus | Foundations of Stochastic Inventory Theory[END_REF] for inventory systems and MDPs.

Finding a sequence of optimal decisions (an optimal policy) is in general a very difficult problem. When the state space and the set of possible decisions (actions) are finite, one can solve the problem numerically. However, this can be done effectively only for very small instances. The aim is therefore to prove some structural properties of an optimal policy, which then simplifies the computation by reducing significantly the search space of all admissible policies.

We are interested here in the class of critical level policies. A critical level policy can be seen as a set of thresholds (a threshold by class) which determine the optimal action to perform (either rejection or acceptance of a demand). Such policies and their applications have been largely studied in the literature We will focus here only on the results closely related to our problem; for a survey with applications to communication networks see Altman [START_REF] Altman | Applications of Markov Decision Processes in communication networks[END_REF] or Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF], and for applications to web services see Mazzucco et al. [START_REF] Mazzucco | Web service hosting and revenue maximization[END_REF].

In the case of an exponential replenishment server, Ha [START_REF] Ha | Inventory rationing in a make-to-stock production system with several demand classes and lost sales[END_REF] has shown the optimality of a critical level policy. This result was further extended by de Véricourt, Karaesmen and Dallery [START_REF] Véricourt | Optimal stock allocation for a capacitated supply system[END_REF], by authorizing the backordering of unsatisfied demands, where critical level policy are also shown to be optimal. Both previous results assume exponential service times. Dekker et al. [START_REF] Dekker | On the (S-1,S) lost sales inventory model with priority demand classes[END_REF], and Kranenburg and van Houtum [START_REF] Kranenburg | Cost optimization in the (S-1, S) lost sales inventory model with multiple demand classes[END_REF] consider general service times and they exhibit exact and heuristic methods to compute the real value of the optimal levels, under the assumption of a critical level policy. Nevertheless, there is no warranty that the critical level policies are optimal under the general service distributions.

In the case of Erlang distributed service times, the optimality of a critical level policy was proven by Ha [START_REF] Ha | Stock rationing in an M/E k /1 make-to-stock queue[END_REF]. However, his model relies on the assumption that the service operator has a constant rate in any state of the model.

We consider here the case of hypoexponential service (see Svoronos [START_REF] Svoronos | Evaluation of one-for-one replenishment policies for multiechelon inventory systems[END_REF] for examples of complex services in inventory problems). We prove the optimality of critical level policies. Note that this represents the simplest case of a phase-type service where the service rate depends on the service phase. Also, we consider a two-dimensional state space that is more suitable to possible extensions to other distributions. This new state-space representation allows us to prove that the switching curve is monotone in k (service phase). We expect that our results represent a first step towards general phase-type service. Our numerical experiments suggest that the result may still be preserved in the hyperexponential case.

This paper is organized as follows. In Section 2 we give the formal description of the MDP problem considered in the paper. Section 3 is devoted to the hypoexponential service case. We prove first in Section 3.3 that the dynamic programming operator preserves convexity, which then has two different applications: i) that the optimal policy is of threshold type and ii) that the thresholds are ordered for different classes. The preservation of submodularity, proved in Section 3.5, implies the monotonicity of the switching curve for a fixed demand class. In Section 4 we discuss the possible extensions of the model. Section 5 contains conclusions and some future directions.

Problem formulation 2.1 Model description

We consider a single-item inventory model. We denote by S the maximum stock level. This single-item is demanded by J different types of customers (or demand classes). Demands arrive according to a Poisson process of parameter λ and the probability that an arriving demand is of class j is equal to p j . Let λ j be the arrival rate of class j. Then λ j = p j λ and J j=1 λ j = λ. When an item is not delivered, the item is lost and a penalty cost is incurred. We denote by c j the penalty cost of demand class j and we assume without loss of generality that the costs are ordered such that c 1 ≤ c 2 ≤ . . . ≤ c J . When a demand arrives, one has to decide whether to satisfy or to refuse it. Since the stock level is finite, rejecting a demand causes a penalty cost but it also saves an item that can be used for future demands (of possibly higher class).

We assume a single replenishment server and the replenishment times can follow different distributions. In this paper we will consider hypoexponential and hyperexponential replenishments.

We detail now how this inventory problem can be expressed as an admission control problem in a queuing system with different classes of customers. We represent the items in replenishment for the inventory system by a queue of capacity S. An arriving demand of class j is an arriving customer of class j. Delivering an item is the equivalent to admitting the customer into the queue. The penalty cost due to a lost demand is equivalent to the cost of rejection. From now on, we will use this queuing model formalism.

Optimal control problem

The problem of finding the best decision (rejection or acceptance), to perform when customers arrive is an optimal control problem that can be modeled as a (time-homogeneous) MDP. We detail now the elements of this MDP.

The system behaves in continuous time, but we can restrict our attention to the transition epochs. This can be done without loss of generality since we only consider costs that depend on the state of the system. At these transition instants an event occurs and induces a change in the system. When the event is an arrival, then we have to choose the best action to trigger. Although the system behaves continuously, we consider here a discrete time MDP. Indeed, the use of the uniformization standard method makes the problem discrete (see Lippman [START_REF] Lippman | Applying a new device in the optimization of exponential queueing systems[END_REF] or Puterman [START_REF] Puterman | Markov Decision Processes Discrete Stochastic Dynamic Programming[END_REF]).

The state space of our MDP is X = {1, . . . , S} × {1, . . . , N } ∪ {(0, 1)}. For any (x, k) ∈ X , the parameter x denotes the number of customers in the system, while k describes the phase of the customer being served (hypoexponential and hyperexponential have different phases). Actually, in state (0, 1) no customer is served and the indication of the phase is useless but we keep it for readability. Note that the state space does not depend on the classes of customers.

We denote by A = {0, 1} the action set where action 1 is an acceptance and action 0 is a rejection. A policy π is a sequence of decision rules that maps the information history (past states and actions) to the action set A. Since in our case the action set and the state space are finite, and since the costs are thus trivially bounded, we restrict our attention to Markov deterministic policies. Indeed, standard results (see Puterman [START_REF] Puterman | Markov Decision Processes Discrete Stochastic Dynamic Programming[END_REF]) insure that an optimal policy belongs to the set of Markov deterministic policies. A Markov deterministic policy is of the form (a(•), a(•), . . .) where a(•) is a single deterministic decision rule that maps the current state to a decision (hence, in our case a(•) is a function from X to A). Under a policy π the evolution of the system generates a random sequence of states and actions y n and a n .

We denote by C(y, a) with y ∈ X and a ∈ A the instantaneous cost in state y when action a is triggered. Here C(y, 0) = c j according to the class of the customer we reject, and C(y, 1) = 0 in case of acceptance.

Since we are interested in the average cost criteria, we define the objective, known as the minimal long-run average costs, by v * = min

π lim n→∞ 1 n E π y n-1 ℓ=0 C(y ℓ , a ℓ ) ,
where the initial state y 0 = y. Our aim is to characterize the optimal policy π * that reaches v * . This means to determine the optimal decision rule, denoted by a * , with π * = (a * , a * , . . .).

Total cost criteria.

We adopt the framework of Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF] and we will work with the total cost criteria to characterize the optimal average policy. We define the minimal (expected) n-stage total cost by

V n (y) = min π(n) E π(n) y n-1 ℓ=0 C(y ℓ , a ℓ ) , ∀y ∈ X , (1)
where π(n) is a policy over n steps. We assume that the initial state y 0 = y and that, by convention, V 0 (y) = 0, ∀y ∈ X . We call V n the n-stage total cost value function.

From Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF] and Puterman [START_REF] Puterman | Markov Decision Processes Discrete Stochastic Dynamic Programming[END_REF]Chapter 8], the minimal n-stage total cost value function V n does not converge when n tends to infinity. Instead, the difference V n+1 (y) -V n (y) converges to the minimal long run average cost. On the other hand, still from Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF], the n-stage policy which minimizes V n tends to the optimal average policy π * when n tends to infinity.

Finally, since the instantaneous costs are positive, from standard results on MDPs (see Puterman [START_REF] Puterman | Markov Decision Processes Discrete Stochastic Dynamic Programming[END_REF]), V n is well defined for all n and satisfies the dynamic programming equation (Bellman equation):

V n+1 = T V n , (2)
where T is the dynamic programming operator. More precisely, T operates on the total cost value functions and is defined by:

(T f )(y) = min a ( T f )(y, a) = min a   C(y, a) + y ′ ∈X P y ′ |(y, a) f (y ′ )   , (3)
with P (y ′ |(y, a)) the probability to move to y ′ from state y ∈ X with action a ∈ A, and f a function that belongs to the total cost value function set. The optimal decision rule can then be deduced from Eq. ( 3):

a(y) = min arg min a∈A ( T f )(y, a) . (4)

Critical level policies

In this work, we do not want to numerically compute the optimal policy but rather to characterize some "good" properties that the optimal policy should verify. Indeed, very often the optimal policy admits a simple form. The aim of our work is to prove that the optimal policy for our model is a critical level policy.

Definition 1 (Critical level policy). A policy is called a critical level policy if for any fixed k and any customer class j it exists a level in x: t k,j , such that in state (x, k):

-for all 0 ≤ x < t k,j it is optimal to accept any customer of class j, -for all x ≥ t k,j it is optimal to reject any customer of class j.

Structured policy framework.

In order to show that the optimal policy exhibits some special properties we use a standard framework presented in the monograph of Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF]. Most of the results presented now follow from this monograph in which detailed proofs can be found. This framework is that of property preservation by the Dynamic Programming operator. It consists of three steps: First, identify two related sets. The one of the value functions that exhibit a special form (say V) and the one one of the decision rules D that exhibit a special form. The relation between these two sets being such that when the value function belongs to V, then the optimal decision belongs to D. Next, show that the properties of v ∈ V are conserved by the operator T . Finally, check that these properties are kept when passing to the limit.

We present now some special forms for value functions. Let f be a function from X → R + .

Definition 2 (Monotonicity). We say that a function f is increasing in x (denoted by Incr(x)) if, for all y = (x, k), we have f (x + 1, k) ≥ f (x, k). Similarly, we say that a function f is decreasing in k (denoted by Decr(k)) if, for all y = (x, k), we have f (x, k + 1) ≤ f (x, k).

Definition 3 (Convexity). We say that a function f is convex in x (denoted by Convex(x)) if, for all y = (x, k), we have:

2f (x + 1, k) ≤ f (x, k) + f (x + 2, k) . (5)
Definition 4 (Submodularity). We say that a function f is submodular in x and k (denoted by Sub(x, k) if, for all y = (x, k), we have:

f (x + 1, k + 1) + f (x, k) ≤ f (x + 1, k) + f (x, k + 1) . (6)
We give now the relationship between the properties of the value functions and the optimal decisions. Proposition 5 (Th 8.1 [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF]). Let a(y), with y = (x, k), be the optimal decision rule defined by Eq. (4

). i) If f ∈ Convex(x), then a(y) is decreasing in x. ii) If f ∈ Sub(x, k), then a(y) is increasing in k.
Assume a critical level policy and consider a decision for a fixed demand class j. Definition 6 (Switching curve). For every k, we define a level t(k) such that when we are in state (x, k) decision 1 is taken if and only if x ≤ t(k) and 0 otherwise. The mapping k → t(k) is called a switching curve.

Let us notice that a switching curve separates the state space X in two regions X 0 , in which decision 0 is taken, and X 1 , in which decision 1 is taken.

Definition 7 (Monotone switching curve). We say that a decision rule is of the monotone switching curve type if the mapping k → t(k) is monotone.

Hypoexponential service

In order to show that the optimal level policy is of critical level type we will prove increasingness (Section 3.2) and convexity (Section 3.3) properties, while submodularity (proved in Section 3.5) implies monotonicity of the switching curve. Furthermore, although the problem was given under its most general formulation in Section 2.2, for the sake of readability and in order to simplify the proofs, in the rest of the paper we will decompose the dynamic programming operator according to the event-based approach presented in Glasserman and Yao [START_REF] Glasserman | Monotone Structure in Discrete-Event Systems[END_REF], and Koole [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF]. This means a decomposition of the dynamic programming operator into operators related to a dedicated event and an observation of the system after decision epoch.

Model description

Let us detail the event-based decomposition. The dynamic programming equation (2) becomes:

V n = T unif J i=1 p i T CA(i) (V n-1 ), T D (V n-1 ) , (7)
where V 0 (x, k) ≡ 0, and T unif , T CA(i) and T D are different event operators. If the service (replenishment) time is hypoexponential, then the random value R that describes the service time is defined by R = N k=1 E k , where E k is an exponential random value of parameter µ k , 1 ≤ k ≤ N . Let α = max k µ k . We will use the uniformization rate equal to λ + α. The uniformization operator, T unif , is defined by:

T unif (f (x, k), g(x, k)) = λ λ + α f (x, k) + α λ + α g(x, k).
Let µ ′ k = µ k /α, the departure operator, T D f , is defined by:

T D f (x, k) = µ ′ k f (x, k + 1) if (k < N ) and (x > 0), f ((x -1) + , 1) if (k = N ) or (x = 0 and k = 0) + (1 -µ ′ k )f (x, k).
The controlled arrival operator of a customer of class i, T CA(i) , is defined by:

T CA(i) f (x, k) = min{f (x + 1, k), f (x, k) + c i } if x < S, f (x, k) + c i if x = S.

Increasing property

We first notice that the state space admits a linear order (denoted ) which is a lexicographic order defined to be increasing with respect to the first dimension and decreasing with respect to the second dimension:

Definition 8 (Linear order ). For every (x, i) and (y, j) in X ,

(x, i) (y, j) ⇔ (x < y) or (x = y, i ≥ j) .
Let us define the increasing property of functions with respect to .

Definition 9 (Increasing property). For any function f :

X → R we say that f ∈ Incr( ) if (x, i) (y, j) ⇒ f (x, i) ≤ f (y, j).
We now show that operators T CA(i) and T D preserve this increasing property.

Lemma 10. Let f : X → R such that f ∈ Incr( ), then T CA(i) f ∈ Incr( ).
Proof. Assume that f ∈ Incr( ). In order to prove that T CA(i) f ∈ Incr( ), it is sufficient to prove that:

T CA(i) f (x, k) ≤ T CA(i) f (x, k -1), ∀x > 0, ∀k > 1, (8)
T CA(i) f (x, 1) ≤ T CA(i) f (x + 1, N ), ∀x, ∀k. (9)
Let us consider first [START_REF] Koole | Monotonicity in Markov reward and decision chains: Theory and applications[END_REF]. For 0 < x < S and k > 1 we get:

T CA(i) f (x, k) = min{f (x, k) + c i , f (x + 1, k)} ≤ min{f (x, k -1) + c i , f (x + 1, k -1)} = T CA(i) f (x, k -1),
For x = S and k > 1 we have:

T CA(i) f (S, k) = f (S, k) + c i ≤ f (S, k -1) + c i = T CA(i) f (S, k -1)
.

We now prove [START_REF] Kranenburg | Cost optimization in the (S-1, S) lost sales inventory model with multiple demand classes[END_REF]. For x < S -1 we get:

T CA(i) f (x, 1) = min{f (x, 1) + c i , f (x + 1, 1)} ≤ min{f (x + 1, k) + c i , f (x + 2, k)} = T CA(i) f (x + 1, k) .
For x = S -1,

T CA(i) f (S -1, 1) = min{f (S -1, 1) + c i , f (S, 1)} ≤ f (S, k) + c i = T CA(i) f (S, k) . Lemma 11. Let f : X → R such that f ∈ Incr( ), then T D f ∈ Incr( ).
Proof. Assume that f is Incr( ). In order to prove that T D f ∈ Incr( ), it is sufficient to prove that:

T D f (x, k) ≤ T D f (x, k -1), ∀x > 0, ∀k > 1, (10)
T D f (x, 1) ≤ T D f (x + 1, N ), ∀x, ∀k. (11)
Let us consider first [START_REF] Lippman | Applying a new device in the optimization of exponential queueing systems[END_REF]. For x > 0 and k = N ,

T D f (x, N ) = µ ′ N f (x -1, 1) + (1 -µ ′ N )f (x, N ) ≤ µ ′ N -1 f (x, N ) + (1 -µ ′ N -1 )f (x, N -1) = T D f (x, N -1)
. Indeed, using Definition 9, one has that both f (x-1, 1) and f (x, N ) are smaller or equal than f (x, N ) and f (x, N -1) and the same holds for any of their convex combinations. Similarly, for x > 0 and 1 < k < N we have:

T D f (x, k) = µ ′ k f (x, k + 1) + (1 -µ ′ k )f (x, k) ≤ µ ′ k-1 f (x, k) + (1 -µ ′ k-1 )f (x, k -1) = T D f (x, k -1)
. Let us now consider [START_REF] Mazzucco | Web service hosting and revenue maximization[END_REF]. For x > 0:

T D f (x, 1) = µ ′ 1 f (x, 2) + (1 -µ ′ 1 )f (x, 1) ≤ µ ′ N f (x, 1) + (1 -µ ′ N )f (x + 1, N ) = T D f (x + 1, N ) For x = 0: T D f (0, 1) = f (0, 1) = T D f (1, N ).
Remark 12. From Definition 8, it follows that Incr( ) implies Incr(x) and Decr(k).

Convexity

In this part we prove that the convexity with respect to the first dimension is preserved by operators T CA(i) and T D . However, due to the state space form, this requires to introduce an additional property (called augmented convexity) at state (0, 1).

Definition 13 (Augmented convexity). We say that a function f : X → R is augmented convex in x (denoted by AConvex(x)) if for all k ∈ {1, .., N } we have f (0, 1)

+ f (2, k) ≥ 2f (1, k). Lemma 14. Let f be a function such that f ∈ Incr( ) ∩ Convex(x) ∩ AConvex(x), then T CA(i) f ∈ Incr( ) ∩ Convex(x) ∩ AConvex(x). Proof. Assume that f is Incr( ) ∩ Convex(x) ∩ AConvex(x).
(i) Preservation of Convex(x) shall be proved for x > 0. For x = S -2,

T CA(i) f (x,k)+T CA(i) f (S, k) = min f (x+ 1, k)+ f (S, k) + c i f (x, k)+ c i +f (S, k)+ c i ≥ 2f (x+1, k)+2c i ≥2 min{f (S, k), f (x+1, k)+c i } = 2T CA(i) f (x+1, k). For 0 < x < S -2, the sum T CA(i) f (x, k) + T CA(i) f (x + 2, k) is equal to min        f (x + 1, k) + f (x + 3, k) f (x + 1, k) + f (x + 2, k) + c i f (x, k) + c i + f (x + 3, k) + f (x + 2, k) -f (x + 2, k) f (x, k) + c i + f (x + 2, k) + c i ≥ min        2f (x + 2, k) f (x + 1, k) + f (x + 2, k) + c i f (x + 1, k) + f (x + 2, k) + c i 2f (x + 1, k) + 2 c i ≥ 2 min{f (x + 2, k), f (x + 1, k) + c i } = 2T CA(i) f (x + 1, k).

RR n°7732

Note that for the third case we used:

f (x, k) + c i + f (x + 3, k) + f (x + 2, k) -f (x + 2, k) ≥ 2f (x + 1, k) + f (x + 3, k) -f (x + 2, k) + c i .
(ii) For x = 0, we show that augmented convexity is preserved:

T CA(i) f (0, 1)+T CA(i) f (2, k) = min                f (1, 1)+f (3, k) ≥ f (1, k)+f (3, k) ≥ 2f (2, k) f (1, 1) + f (2, k) + c i f (0, 1) + c i + f (3, k) + f (2, k) -f (2, k) ≥ 2f (1, k) + f (3, k) -f (2, k) + c i ≥ f (1, k) + f (2, k) + c i f (0, 1) + c i + f (2, k) + c i ≥ 2f (1, k) + 2 c i ≥ 2 min{f (2, k), f (1, k) + c i } = 2T CA(i) f (1, k),
Preservation of Incr( ) follows from Lemma 10.

Lemma 15. Let f be a function such that f ∈ Incr( ) ∩ Convex(x) ∩ AConvex(x), then T D f ∈ Incr( ) ∩ Convex(x) ∩ AConvex(x). Proof. Assume that f is in Incr( )∩ Convex(x) ∩ AConvex(x).
(i) Preservation of Convex(x) shall be proved for x > 0. For x > 0 and k < N :

T D f (x, k) + T D f (x + 2, k) = µ ′ k f (x, k + 1) + (1 -µ ′ k )f (x, k) + µ ′ k f (x + 2, k + 1) + (1 -µ ′ k )f (x + 2, k) ≥ 2µ ′ k f (x + 1, k + 1) + 2(1 -µ ′ k )f (x + 1, k) = 2T D f (x + 1, k)
For x > 0 and k = N :

T D f (x, N ) + T D f (x + 2, N ) = µ ′ N f (x -1, 1) + (1 -µ ′ N )f (x, N ) + µ ′ N f (x + 1, 1) + (1 -µ ′ N )f (x + 2, N ) ≥ 2µ ′ N f (x, 1) + 2(1 -µ ′ N )f (x + 1, N ) = 2T D f (x + 1, N )
(ii) For x = 0 augmented convexity is preserved. For x = 0 and k < N :

T D f (0, 1)+T D f (2, k) = µ ′ 1 f (0, 1)+(1-µ ′ 1 )f (0, 1)+µ ′ k f (2, k+1) + (1-µ ′ k )f (2, k) = µ ′ k f (0, 1)+(1-µ ′ k )f (0, 1)+µ ′ k f (2, k+1)+(1-µ ′ k )f (2, k) ≥ 2µ ′ k f (1, k + 1) + 2(1 -µ ′ k )f (1, k) = 2T D f (1, k)
For x = 0 and k = N :

T D f (0, 1)+T D f (2, N ) = µ ′ 1 f (0, 1)+(1-µ ′ 1 )f (0, 1)+µ ′ N f (1, 1)+(1-µ ′ N )f (2, N ) = µ ′ N f (0, 1)+µ ′ N f (1, 1)+(1-µ ′ N )f (0, 1)+(1-µ ′ N )f (2, N ) ≥ 2µ ′ N f (0, 1) + 2(1 -µ ′ N )f (1, N ) = 2T D f (1, N )
Preservation of Incr( ) follows from Lemma 11.

We aggregate now all the event-based operators to get the preservation of structural properties passing through the whole dynamic programming operator.

Proposition 16. Let V n be a n-steps total cost value function that satisfy Equation 7. Then, for all n ≥ 0, V n is in Incr( ) ∩ AConvex(x) ∩ Convex(x).

Proof. First, let us note that if for all s.t. i ∈ {1, ..., n} we have f i ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1) then we have

n i=1 p i f i ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1) . ( 12
)
This follows from the fact that we have a direct sum of functions each of them satisfying the properties. Furthermore, since the operator T unif is a convex sum of two functions, we get that if f, g ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1) then

T unif (f, g) ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1) (13)
Thus, combination of Lemmas 14, 15 with observations ( 12) and ( 13) yields

V n ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1) ⇒ V n+1 ∈ Incr( ) ∩ AConvex(1) ∩ Convex(1)
which, along with the fact that V 0 ≡ 0 is in Incr( ) ∩ AConvex(1) ∩ Convex(1), proves the result.

Critical level policy

The next step is to show that the optimal policy, as defined in Section 2.2.1, is a critical level policy (see Definition 1).

Theorem 17. The optimal policy is a critical level policy.

Proof. First, any value function V n satisfying Eq. ( 7) belongs to Convex(x) this follows from Proposition 16. From Proposition 5 it follows that, for any n, the action minimizing V n when V n ∈ Convex(x) defines a critical level policy. Finally, from the discussion in 2.2.1 it follows that the policy minimizing V n converges (with a simple convergence) to the optimal policy when n tends to infinity. Since convexity is kept with simple convergence, then the optimal policy is of critical level policy type. Numerical computation of the optimal policy using value iteration method results in a critical level policy as presented in Figure 1. It can be seen that for every customer class j and for every phase k there exists a unique threshold in x: t k,j . These thresholds are represented on the figure as the transition between acceptance for different classes (blue circle -all classes are accepted, green triangle -classes 2 and 3 are accepted, pink square -only class 3 is accepted, red asterisk -rejection of any class). This example and Figure 1 also serve as an illustration for Theorems 19 and 24. Now we show that when the rejection costs are ordered for consecutive customer classes, the critical levels in x are ordered. Theorem 19. For any critical level policy (see Definition 1), if the rejection costs are nondecreasing (c 1 ≤ • • • ≤ c J ), then the levels t k,j are nondecreasing with respect to customer class j, i.e. t k,j ≤ t k,j+1 .

Proof. First recall that the operator T CA(j) is defined as: T CA(j) f (x, k) = min {f (x+1, k), f (x, k) + c j } for x < S (for x = S the only allowed action is rejection).

Observe now that by definition t k,j is a number such that: in state (x, k) for all 0 ≤ x < t k,j it is optimal to accept any customer of class j and to reject it for all x ≥ t k,j . This, by definition of T CA(j) , means that:

V n (x + 1, k) -V n (x, k) < c j ∀x < t k,j (14)
V

n (x + 1, k) -V n (x, k) ≥ c j ∀x ≥ t k,j
Owing to the assumption c j ≤ c j+1 , it yields

V n (x + 1, k) -V n (x, k) < c j ≤ c j+1
which, by definition of t k,j+1 holds for all x < t k,j+1 . But, by Eq. ( 14), this also holds for all x < t k,j . This means that the set of x such that x < t k,j is a subset of the set of x such that x < t k,j+1 , which, by definition of t k,j implies that t k,j ≤ t k,j+1 .

Submodularity

In this part we prove the preservation by operators T CA(i) and T D of submodularity with respect to the first and second dimension (Sub(x, k)). Similarly to convexity, we need an additional property (boundary-submodularity) and subsequently its preservation along with submodularity shall be proved too.

Definition 20 (Boundary-submodularity). We say that a function

f is in BSub(x, k) if ∀0 < x < S we have f (x, 1) + f (x, N ) ≤ f (x -1, 1) + f (x + 1, N ). Lemma 21. Let f be a function such that f ∈ Sub(x, k) ∩ BSub(x, k). Then T CA(i) f ∈ Sub(x, k) ∩ BSub(x, k).
Proof. The proof is postponed in Appendix A.2.

Lemma 22. Let f be a function such that f ∈ Sub(x, k)∩BSub(x, k). Then T D f ∈ Sub(x, k)∩ BSub(x, k).
Proof. The proof is postponed in Appendix A. [START_REF] Dekker | On the (S-1,S) lost sales inventory model with priority demand classes[END_REF]. Also here we prove that the submodularity properties are kept by passing through the whole dynamic programming operator.

Proposition 23. Let V n be a n-step total cost value function which satisfies Eq. [START_REF] Kleywegt | The dynamic and stochastic knapsack problem[END_REF]. Then for all (n ≥ 0) V n is in Sub(x, k) ∩ BSub(x, k).

Proof. Using the same arguments than in the proof of Proposition 16, we obtain that if for all i ∈ {1, . . . , n} we have

f i ∈ BSub(1,2) ∩ Sub(1,2) then n i=1 p i f i ∈ BSub(1,2) ∩ Sub(1,2) . (15)
We also obtain that if f, g ∈ BSub(1,2) ∩ Sub(1,2), then

T unif (f, g) ∈ BSub(1,2) ∩ Sub(1,2) (16)
Combination of Lemmas 21, 22 with observations (15) and ( 16) yields

V n ∈ BSub(1,2) ∩ Sub(1,2) ⇒ V n+1 ∈ BSub(1,2) ∩ Sub(1,2)
which, along with the fact that V 0 ≡ 0 is in BSub(1,2) ∩ Sub(1,2) proves the result.

Monotone switching curve

Henceforth we can now refine the optimal policy using Definition 6.

Theorem 24. The optimal policy defines a monotone switching curve.

Proof. For all n ≥ 0, a value function V n that checks Eq. ( 7) belongs to Sub(x, k) (this follows from Proposition 23). Thus, by Theorem 17, we know that action is decreasing in x and that there exists a critical level policy. Since the action set consists of two elements and since the switching curve is the boundary between inverse images of these two elements with respect to the Markov deterministic policy, then Theorem 17 yields that for a fixed customer class the boundary (switching curve) is a function k → t k,j . Furthermore, from Prop. 5 it follows that the action minimizing V n for V n ∈ Sub(x, k) is increasing in k. This, combined with the previous observation, implies that the switching curve is increasing. Observe now, that (from the discussion in 2.2.1) the policy minimizing V n converges to the optimal policy when n tends to infinity. Since submodularity is kept by simple convergence, thus the optimal policy is a monotone switching curve.

At last, it has to be noted that submodularity is not well defined for the state (0, 1) (indeed when x = 0 only one phase is possible). Therefore the same reasoning based on Proposition 5 cannot be applied to the zero state. However, Theorem 17 ensures that an unique level t 1,j for k = 1 exists. Submodularity, by application of Proposition 5, implies the fact that the switching curve is increasing for x > 0. Finally, the only event not covered by the previous discussion (i.e. the existence of the switching point at (0, 1)) does not impair the monotonicity of the switching curve, since it is the smallest possible state.

Model extensions

We discuss the possible extensions of the model in two different directions: hyperexponential service times and including different instantaneous costs.

Hyperexponential service

We now assume that the service times follow a hyperexponential distribution (an exponential distribution of rate µ k with a probability γ k ). Let V n be the n-step value function, the dynamic programming equation is similar:

V n (x, k) = T unif J i=1 p i T CA(i) (V n-1 ), T D (V n-1
) , with V 0 (x, k) ≡ 0. However, since the dynamics of the system changes compared to the hypoexponential case, the definitions of controlled arrivals T CA(i) and departure operators T D differ.

Let us recall that α is the maximum over all µ k and that µ ′ k = µ k α . We assume, without loss of generality, that µ 1 ≤ µ 2 ≤ . . . ≤ µ N . The departure operator T D is now defined by:

T D f (x, k) = µ ′ k N i=1 a i f (x -1, i) if x > 1 f (0, 1) if x ≤ 1 + (1 -µ ′ k )f (x, k) ,
and the controlled arrivals T CA(i) by: (rejection costs). The optimal policy is given in Figure 2. For this example, and many others that we have tested, the optimal policy is still of critical level type. This suggests that even in the hyperexponential service case, the optimal policies are critical level policies. However, since the arrival operator is now more complex at the zero state, the proof technique used for hypoexponential case cannot be used directly, and were not able to prove any results yet about critical level optimality.

T CA(i) f (x, k) =      f (x, k) + c i if x = S, min{f (x + 1, k), f (x, k) + c i } if 0 < x < S, min{f (0, 1) + c i , N j=1 a j f (1, j)} if x = 0.

Adding holding costs

We now discuss the possibility of adding holding costs to Eq. [START_REF] Kleywegt | The dynamic and stochastic knapsack problem[END_REF]. As mentionned in De Véricourt et al. [START_REF] Véricourt | Optimal stock allocation for a capacitated supply system[END_REF], this breaks the queueing model and the inventory system similarities. Thus the two models should be studied separately. The admission control in a queueing system model with holding costs has a wide range of applications (see Kleywegt et al. [START_REF] Kleywegt | The dynamic and stochastic knapsack problem[END_REF]). In the context of inventory systems, it models the make-to-stock problem (see Ha [START_REF] Ha | Inventory rationing in a make-to-stock production system with several demand classes and lost sales[END_REF]).

We obtain the similar results as before when we add linear and increasing holding costs in admission control model.

Corollary 26 (Critical level policy with holding costs). Let T cost be the operator defined by

T cost f (x, k) = c h x λ + α + f (x, k) ,
with c h the per capita holding cost. The value function is given by V

n (x, k) = T cost (T p (V n-1 ))
where T p is the operator given by Eq. ( 7). The optimal policy is of critical level policy type. The levels are increasing with respect to the customer class and the optimal policy defines a monotone switching curve.

Proof. We only give a sketch of proof by noticing that the increasing, convexity, and submodularity properties are kept by the operator. Then the results follow using similar arguments as in the proofs of Theorems 17, 19, and 24.

Conclusions

The main contribution of this paper is the proof that in the case of a hypoexponential service case we still have an optimal policy of critical level type. This extends the previous results from the literature to the simplest case of a phase-type service where the service rate depends on the service phase. We expect that our result is a first step towards more general results for phase-type service. Our numerical experiments suggest that the result is still preserved for hyperexponential case. However, the particularity of the zero state (i.e. when the queue is empty) becomes more challenging because of the possibility of choice in the arrival operator when the queue is empty. Technically, on the boundary around zero state convexity and submodularity must be replaced with some more elaborate properties. Defining the boundary properties for general phase type service and finding the conditions under which they are preserved by the dynamical programming operator represents our first direction for future extensions of this work. The second direction concerns the extensions of the result for different types of immediate costs. Finally, it would be interesting to consider other types of arrival processes or adding the possibility of decision for service.

Acknowledgements.

This work was supported by INRIA ARC OCOQS grant.

A Proofs A.1 Preliminary Lemma

The following observation is used extensively within the proofs:

Lemma 27. For any real numbers A ≤ B ≤ C and 0 ≤ θ 1 , θ 2 ≤ γ we have θ 1 A + (γ -θ 1 )B ≤ θ 2 B + (γ -θ 2 )C.

A.2 Proof of Lemma 21

Assume that f is Sub(x, k) ∩ BSub(x, k). Propagation of submodularity (i) will be proved first, followed by the propagation of boundary-submodularity (ii).

(i) Sub is preserved by T CA(i) . For x < S -1

T CA(i) f (x, k) + T CA(i) f (x + 1, k + 1) = min                        f (x,

A.3 Proof of Lemma 22

Assume that f is Incr( ) ∩ Sub(x, k) ∩ BSub(x, k). Propagation of Incr( ) follows from Lemma 11. Propagation of submodularity (i) will be proved first, followed by the propagation of boundary submodularity (ii). Both properties need to be considered in two cases which makes the proof consist of four parts. All the parts use the same reasoning and make use of Lemma 27 feeding in different inequalities.

(i) Sub is preserved by T D . For x > 0 and k < N -1 in order to prove propagation of Sub it has to be shown that:

T D f (x, k) + T D f (x + 1, k + 1) = µ ′ k f (x, k + 1) + (1 -µ ′ k )f (x, k) + µ ′ k+1 f (x + 1, k + 2) + (1 -µ ′ k+1 )f (x + 1, k + 1) ≤ µ ′ k+1 f (x, k + 2) + (1 -µ ′ k+1 )f (x, k + 1) + µ ′ k f (x + 1, k + 1) + (1 -µ ′ k )f (x + 1, k) = T D f (x, k + 1) + T D f (x + 1, k) (17)
Submodularity, given as assumption, at f (x, k) and f (x, k + 1) respectively can be rewritten as:

f (x, k) -f (x + 1, k) ≤ f (x, k + 1) -f (x + 1, k + 1) (18)
f (x, k + 1) -f (x + 1, k + 1) ≤ f (x, k + 2) -f (x + 1, k + 2) (19)

Example 18 .

 18 Consider a model with the following parameters as an example: J = 3 (number of customer classes), N = 5 (number of phases), S = 10 (stock size), λ = 3, µ = [2, 6, 9, 4, 7] (service rates in different phases), p = [0.3, 0.4, 0.3] (probabilities of J customer classes), c = [30, 40, 50] (rejection costs).

Figure 1 :

 1 Figure 1: Acceptance points for different customer classes. Blue circle -all classes are accepted, green triangle -classes 2 and 3 are accepted, pink square -only class 3 is accepted, red asterisk -rejection of any class.

Example 25 .

 25 Consider a model with the following parameters: J = 3 (number of customer classes), N = 5 (number of phases), S = 10 (stock size), λ = 3, µ = [2, 4, 6, 7, 9] (service rates for different phases), p = [0.3, 0.4, 0.3] (probabilities of J customer classes), c = [0.5, 1, 3]

Figure 2 :

 2 Figure 2: Acceptance points for different customer classes. Blue circle -all classes are accepted, green triangle -classes 2 and 3 are accepted, pink square -only class 3 is accepted, red asterisk -rejection of any class.

1 T 1 Tf 1 Tf

 111 k) + f (x + 1, k + 1) + 2c i ≤ f (x, k + 1) + f (x + 1, k) + 2c i f (x, k) + f (x + 2, k + 1) + c i = f (x, k) + f (x + 2, k + 1) + f (x + 1, k + 1) -f (x + 1, k + 1) + c i ≤ f (x + 1, k) + f (x, k + 1) + f (x + 2, k + 1) -f (x + 1, k + 1) + c i ≤ f (x + 2, k) + f (x + 1, k + 1) + f (x, k + 1) -f (x + 1, k + 1) + c i ≤ f (x, k + 1) + f (x + 2, k) + c i f (x + 1, k) + f (x + 1, k + 1) + c i f (x + 1, k) + f (x + 2, k + 1) ≤ f (x + 1, k + 1) + f (x + 2, k) = min{f (x, k + 1) + c i , f (x + 1, k + 1)} + min{f (x + 1, k) + c i , f (x + 2, k)} = T CA(i) f (x, k + 1) + T CA(i) f (x + 1, k),and for x = S -CA(i) f (x, k) + T CA(i) f (S, k + 1) = min f (x, k) + f (S, k + 1) + 2c i ≤ f (x, k + 1) + f (S, k) + 2c i f (S, k) + f (S, k + 1) + c i = f (S, k) + c i + min{f (S, k + 1) + c i , f (x, k + 1)} = T CA(i) f (S, k) + T CA(i) f (x, k + 1). RR n°7732 (ii) BSub is preserved by T CA(i) . For x < S -CA(i) f (x, N ) + T CA(i) f (x(x, N ) + f (x, 1) + 2c i ≤ f (x -1, 1) + f (x + 1, N ) + 2c i f (x, N ) + f (x + 1, 1) + c i = f (x, N ) + f (x + 1, 1) + f (x, 1) -f (x, 1) + c i ≤ f (x, 1) + f (x + 2, N ) + f (x -1, 1) -f (x, 1) + c i = f (x + 2, N ) + f (x -1, 1) + c i f (x + 1, N ) + f (x, 1) + c i f (x + 1, N ) + f (x + 1, 1) ≤ f (x, 1) + f (x + 2, N ) = min{f (x -1, 1) + c i , f (x, 1)} + min{f (x + 1, N ) + c i , f (x + 2, N )} = T CA(i) f (x -1, 1) + T CA(i) f (x + 1, N ), and for x = S -CA(i) f (x, N )+T CA(i) f (x, 1) = min (x, N ) + f (x, 1) + 2c i ≤ f (x -1, 1) + f (S, N ) + 2c i f (x, N ) + f (S, 1) + c i f (S, N ) + f (x, 1) + c i f (S, N ) + f (S, 1)≤ min f (x -1, 1) + f (S, N ) + 2c i f (x, 1) + f (S, N ) + c i = f (S, N ) + c i + min{f (x -1, 1) + c i , f (x, 1)} = T CA(i) f (S, N ) + T CA(i) f (x -1, 1).

RR n°7732

Now, since the right hand side of inequality (18) is equal to the left hand side of inequality (19), Lemma 27 can be used with θ 1 = (1 -µ ′ k ) and θ 2 = (1 -µ ′ k+1 ) and γ = 1. Application of this lemma yields:

which is equivalent to the inequality (17) to be proved (after rearrangement of terms so that every value of f is positive). For x > 0 and k = N -1 in order to prove propagation of Sub it has to be shown that:

Submodularity, given as assumption, at f (x, k) and boundary-submodularity, also given as assumption, at f (x, N ) can be rewritten as:

Now, since the right hand side of inequality ( 21) is equal to the left hand side of inequality (22), Lemma 27 can be used with θ 1 = (1 -µ ′ k ) and θ 2 = (1 -µ ′ N ) and γ = 1. Application of this lemma yields:

which is equivalent to the inequality (20) to be proved.

(ii) BSub is preserved by T D . For x > 1 in order to prove propagation of BSub it has to be shown that:

Boundary-submodularity, given as assumption, at f (x, N ) and submodularity, also given as assumption, at f (x -1, 1) respectively can be rewritten as:

Now, since the right hand side of inequality ( 24) is equal to the left hand side of inequality (25), Lemma 27 can be used with θ 1 = (1 -µ ′ N ) and θ 2 = (1 -µ ′ 1 ) and γ = 1. Application of this lemma yields:

which is equivalent to the inequality (23) to be proved. For x = 1 in order to prove propagation of BSub it has to be shown that:

Boundary-submodularity, given as assumption, at f (1, N ) and increasing property, also given as assumption, at f (1, 1) (f (1, 1) ≥ f (1, 2)) respectively can be rewritten as:

Now, since the right hand side of inequality ( 27) is equal to the left hand side of inequality (28), Lemma 27 can be used with θ 1 = (1 -µ ′ N ) and θ 2 = (1 -µ ′ 1 ) and γ = 1. Application of this lemma yields:

which is equivalent to the inequality (26) to be proved.