Sophie Mercier
email: sophie.mercier@univ-pau.fr

I T Castro

Stochastic comparisons of imperfect maintenance models for a gamma deteriorating system

Keywords: Reliability, non homogeneous gamma process, (increasing) convex order, virtual age model

This paper compares two imperfect repair models for a degrading system, with deterioration level modeled by a non homogeneous gamma process. Both models consider instantaneous and periodic repairs. The first model assumes that a repair reduces the degradation of the system accumulated from the last maintenance action. The second model considers a virtual age model and assumes that a repair reduces the age accumulated by the system since the last maintenance action. Stochastic comparison results between the two resulting processes are obtained. Furthermore, a specific case is analyzed, where the two repair models provide identical expected deterioration levels at maintenance times. Finally, two optimal maintenance strategies are explored, considering the two models of repair.

Introduction

Safety and dependability are crucial issues in many industries, which have lead to the development of a huge literature devoted to the so-called reliability theory. In the oldest literature, the lifetimes of industrial systems or components were usually directly modeled through random variables, see, e.g., [START_REF] Barlow | Mathematical theory of reliability[END_REF] for a pioneer work on the subject. In case of repairable systems, successive lifetimes of a system then appear as the points of a counting process leading to so-called recurrent events. Based on the development of on-line monitoring which allows the effective measurement of a system deterioration (length of a crack, thickness of a cable, intensity of vibrations, temperature, ...), numerous papers nowadays model the degradation in itself, which is often considered to be (mostly) monotonous with respect to the time. This is done through the use of stochastic processes such as Wiener processes with trend [START_REF] Liu | A condition-based maintenance policy for degrading systems with age-and state-dependent operating cost[END_REF]; [START_REF] Zhang | Degradation-based maintenance decision using stochastic filtering for systems under imperfect maintenance[END_REF]; [START_REF] Hu | Optimum step-stress accelerated degradation test for wiener degradation process under constraints[END_REF]), inverse gaussian models [START_REF] Chen | Condition-based maintenance using the inverse gaussian degradation model[END_REF]), transformed Beta degradation processes [START_REF] Giorgio | A new state-dependent degradation process and related model misidentification problems[END_REF]) or gamma processes [START_REF] Huynh | On the use of mean residual life as a condition index for condition-based maintenance decision-making[END_REF]), among others (see also [START_REF] Mercier | A preventive maintenance policy for a continuously monitored system with correlated wear indicators[END_REF] in case of a bivariate deterioration indicator). This paper focuses on gamma processes, which seem the most popular, see [START_REF] Van Noortwijk | A survey of the application of Gamma processes in maintenance[END_REF] with a large number of references therein.

To mitigate the effect of the system degradation and to extend the system lifetime, a large volume of maintenance models have been proposed in the literature. Most of these models are limited to perfect repairs [START_REF] Caballé | A condition-based maintenance of a dependent degradation-threshold-shock model in a system with multiple degradation processes[END_REF]; [START_REF] Hong | Optimal condition-based maintenance decisions for systems with dependent stochastic degradation of components[END_REF]). However, imperfect maintenance actions describe more realistic situations than perfect repairs. Some advances have been made to include imperfect repairs in a degrading system [START_REF] Alaswad | A review on condition-based maintenance optimization models for stochastically deteriorating system[END_REF]; [START_REF] Giorgio | A new state-dependent degradation process and related model misidentification problems[END_REF]). However, as [START_REF] Zhang | Degradation-based maintenance decision using stochastic filtering for systems under imperfect maintenance[END_REF] claimed, the issue of treating imperfect maintenance in the context of degrading systems remains widely open nowadays.

Stochastic orders and related inequalities play an important role in reliability theory and maintenance policies, as they allow to, e.g., obtain bounds for system reliability or availability, or to compare different maintenance strategies [START_REF] Barlow | Comparison of replacement policies and renewal theory implications[END_REF], [START_REF] Ohnishi | Stochastic Models in Reliability and Maintenance, Osaki Edition[END_REF]). There is a huge reliability literature on the use of stochastic orders which compare locations of the lifetime, residual lifetime or inactivity time of the systems [START_REF] Khaledi | Ordering conditional lifetimes of coherent systems[END_REF]). However, there exist other types of stochastic orders which measure variability and spread. Though their use has become classical in insurance literature [START_REF] Denuit | Some new classes of stochastic order relations among arithmetic random variables, with applications in actuarial sciences[END_REF] and [START_REF] Denuit | Lorenz and excess wealth orders, with applications in reinsurance theory[END_REF] amongst others), they are not so common in the reliability literature apart from a few exceptions [START_REF] Kochar | Comparisons of parallel systems according to the convex transform order[END_REF], [START_REF] Fang | On the right spread ordering of series systems with two heterogeneous weibull components[END_REF]).

Following the spirit showed in [START_REF] Mercier | On the modelling of imperfect repairs for a continuously monitored gamma wear process through age reduction[END_REF] and [START_REF] Castro | Performance measures for a deteriorating system subject to imperfect maintenance and delayed repairs[END_REF] (see also [START_REF] Giorgio | A new state-dependent degradation process and related model misidentification problems[END_REF]), two models of imperfect repair are analyzed in this paper for a gamma deteriorating system. The first model, called Arithmetic Reduction of of Deterioration of order 1 (ARD1), assumes that the repair removes the ρ 1 % of the degradation accumulated by the system from the last maintenance action. The second model is based on the notion of virtual age as introduced by [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] in the context of recurrent events (where only lifetime data are available). The idea is that an imperfect repair rejuvenates the system, namely puts it back to a similar state as it was before the repair (details further). Following [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF] in the context of recurrent events, an Arithmetic Reduction of Age of order 1 (ARA1) is here considered, which assumes that the repair removes the ρ 2 % of the (virtual) age accumulated by the system since the last maintenance action. An ARD1 repair hence lowers the deterioration level, without rejuvenating the system. On the contrary, by an ARA1 repair, the system is put back to the exact situation where it was some time before, which entails the lowering of both its deterioration level and (virtual) age. The two models may hence correspond to different maintenance actions in an application context. As an example, [START_REF] Sadeghi | Improvement of railway ballast maintenance approach, incorporating ballast geometry and fouling conditions[END_REF] show that the tamping or cleaning of the ballast of a railway track may have different consequences: The tamping mainly improves the track geometry conditions (short term impact) but has mostly no impact on the ballast mechanical conditions (long term impact), whereas the cleaning also improves the latter. Then, one could think that the tamping corresponds to some reduction of the track deterioration level (such as an ARD1 repair) whereas the cleaning also acts on its potential of future degradation, which could be represented by some reduction of age model (such as an ARA1 repair).

The choice between the two models may however not always be so clear in an applied context. For a better understanding of their differences, this paper focuses on their comparison, from a probabilistic point of view. Assuming that the degradation of the system is modeled by a non homogeneous gamma process, stochastic comparisons of both location and spread of the two resulting processes are given. Moreover, a specific case is analyzed, where the two models provide identical expected deterioration levels at repair times ("equivalent" case).

Going back to the general setup, two maintenance strategies are next proposed. Both strategies consider periodic imperfect repairs (period T ) based on either one of the two models (ARD1 or ARA1). In the first strategy ((n, T ) policy), the system is replaced at the time of the n-th repair. The second strategy ((M, T ) policy) considers a control limit rule, with replacement when the degradation level exceeds a preventive threshold M . For both maintenance strategies, the objective function is the expected profit rate, which takes into account some reward produced by the system, with unitary reward (or cost) per unit time depending on the degradation level of the system. (The lower the deterioration level, the higher the unitary reward per unit time). This reward function is based on classical utility functions used in insurance literature [START_REF] Rolski | Stochastic Processes for Insurance and Finance[END_REF]). The use of this reward function represents an advance in the reliability literature where the cost objective function is usually developed considering that the system state is binary (up or down), with some fixed unavailability cost per unit time, independent on the deterioration level. Theoretical results are obtained for the comparison of the objective functions of the (n, T ) policy under the two types of imperfect repairs (ARD1 or ARA1).

The paper is organized as follows: Section 2 provides some technical reminders. The two imperfect repair models are described in Section 3. The corresponding moments are compared in Section 4 whereas Section 5 is devoted to stochastic comparison results. The "equivalent" case is studied in Section 6. Section 7 deals with the reward function and the two maintenance strategies. Concluding remarks are provided in Section 8, together with possible extensions.

Technical reminders

The definition of two stochastic orders is first recalled, which allows to compare the location of random variables.

Definition 1. Let X and Y be two non negative random variables with probability density functions (p.d.f.) f X and f Y with respect to the Lebesgue measure, cumulative distribution functions (c.d.f.) F X and F Y and survival functions FX and FY , respectively. Then:

1. X is said to be smaller than Y in the usual stochastic order (X ≺ sto Y ) if FX ≤ FY (or F X ≥ F Y , equivalently). 2.
X is said to be smaller than Y in the likelihood ratio order (X ≺ lr Y ) if f Y f X is non-decreasing on the union of the supports of X and Y .

We recall that the likelihood ratio order implies the usual stochastic order. The definition of other stochastic orders is next provided, which allows to compare the variability of two random variables.

Definition 2. Let X and Y be two non negative random variables where the support of X is assumed to be included in the support of Y and the support of Y to be an interval (for the log-concavity). Then:

1. X is said to be smaller than Y in the log-concave order (X ≺ lc Y ) if the ratio f X f Y is log-concave over the support of Y . 2. X is said to be smaller than Y in the convex (concave) order (X ≺ cx(cv) Y ) if E (φ (X)) ≤ E (φ (Y )
) for all convex functions φ (provided the expectations exist).

3. X is said to be smaller than Y in the increasing convex (concave) order (X

≺ icx(icv) Y ) if E (φ (X)) ≤ E (φ (Y )
) for all increasing convex (concave) functions φ (provided the expectations exist).

Following [START_REF] Shaked | Stochastic Orders[END_REF], X ≺ icx (≺ icv )Y roughly means that E(X) ≤ E(Y ) (location condition) plus the fact that X is less (more) "variable" than Y , in a stochastic sense. Also,

X ≺ cx (≺ cv )Y is equivalent to X ≺ icx (≺ icv )Y plus E(X) = E(Y ).
Setting X and Y to be two non negative random variables, we recall (Shaked and Shanthikumar, 2007, page 182)

that X ≺ icx Y if and only if ∫ +∞ x FX (u) du ≤ ∫ +∞ x FY (u) du for all x ≥ 0 (1)
and that X ≺ icv Y if and only if

∫ x 0 F X (u) du ≥ ∫ x 0 F Y (u) du for all x ≥ 0. (2)
Finally, the usual stochastic order (and hence the likelihood ratio order as well) implies both increasing convex and concave orders (Müller and Stoyan, 2002, p.61).

We next come to reminder on gamma distribution. Let a, b > 0. We recall that the gamma distribution Γ (a, b) with parameters (a, b) admits

f a,b (x) = b a Γ (a) x a-1 e -bx 1 R + (x)
as p.d.f. (with respect to Lebesgue measure) and that the corresponding mean and variance are a b and a b 2 , respectively. We shall also make use of the following well-known facts repeatedly, without further notification:

If X is gamma distributed Γ (a, b), then c X is gamma distributed Γ (a, b/c) for all c > 0. If X 1 , • • • , X n are independent gamma distributed random variables with respective distributions Γ (a 1 , b) , • • • , Γ (a n , b), then ∑ n i=1 X i is gamma distributed Γ ( ∑ n i=1 a i , b
). Finally, the following technical result may be found in (Müller and Stoyan, 2002, p. 62).

Lemma 1. Let X and Y be gamma distributed random variables with parameters (a 1 , b 1 ) and (a 2 , b 2 ), respectively, where a i , b i > 0 for i = 1, 2. Then:

1. If a 1 ≤ a 2 and b 1 ≥ b 2 , then X ≺ lr Y ; 2. If a 1 ≥ a 2 and a 1 /b 1 ≤ a 2 /b 2 , then X ≺ icx Y ; 3. If a 1 ≤ a 2 , b 1 ≤ b 2 and a 1 /b 1 ≤ a 2 /b 2 , then X ≺ icv Y .

The two models of imperfect repairs

In the sequel of this work, we set (X t ) t≥0 to be the intrinsic (out of repair) degradation process of the system. We assume that (X t ) t≥0 follows a non homogeneous gamma process with parameters A(•) and b, where A(•) : R + -→ R + is continuous and non-decreasing with A (0) = 0, and b > 0. We recall that (X t ) t≥0 is a process with independent increments such that X 0 = 0 almost surely (a.s.) and such that each increment X t+s -X t is gamma distributed Γ(A(t + s) -A(t), b) for all s, t > 0.

The system is periodically and instantaneously maintained each T units of time. For modeling purpose, we set X (i) , i ∈ N * to be i.i.d. copies of X = (X t ) t≥0 , where X (i) describes the evolution of the deterioration level between the i-th and (i+1)-th maintenance actions. For each imperfect repair model, the maintenance efficiency is measured by an Euclidian parameter ρ ∈ (0, 1).

First model: Arithmetic Reduction of Deterioration of order 1 (ARD1)

In this model, the maintenance action instantaneously removes the ρ% of the degradation accumulated by the system from the last maintenance action (or from the origin). Let (Y t ) t≥0 be the process that describes the degradation level of the maintained system under this model of repair.

The ARD1 model is developed as follows: At the beginning, the system deteriorates according to X (1) and it is first maintained at time T . This provides:

Y t = X (1) t t < T, Y T = (1 -ρ) X (1)
T . Between T and 2T , the system deteriorates according to X (2) . The age of the system is unchanged at time T and we simply have

Y t = Y T + ( X (2) t -X (2) T
) for all T ≤ t < 2T , and at the second maintenance time 2T :

Y 2T = Y T + (1 -ρ) ( X (2) 2T -X (2) T
)

.

More generally, we get:

Y t = Y nT + ( X (n+1) t -X (n+1) nT
) for all nT ≤ t < (n + 1) T, (3)

Y (n+1)T = Y nT + (1 -ρ) ( X (n+1) (n+1)T -X (n+1) nT ) (4)
where

X (n+1) t -X (n+1) nT is gamma distributed Γ(A (t) -A (nT ) , b) for all t ∈ [nT, (n + 1) T ]. Hence Y nT = (1 -ρ) n ∑ i=1 ( X (i) iT -X (i) (i-1)T ) (5) is gamma distributed Γ ( A (nT ) , b 1-ρ
) .

Except for the case ρ → 0 + , if t mod T ̸ = 0 , Y t is the sum of two independent and gamma distributed random variables (r.v.s) with different scale parameters, and it is not gamma distributed. Its expectation and variance are given by:

E (Y t ) = A (t) -ρ A (nT ) b , var (Y t ) = A (t) -ρ (2 -ρ) A (nT ) b 2 , ( 6
)
for nT ≤ t < (n + 1) T .

Remark 1. It is easy to check that E (Y t ) and var (Y t ) are decreasing with respect to ρ.

Next result provides some more insight than the previous remark into the impact of the maintenance efficiency ρ on the deterioration level of the maintained system. To state it, two different efficiency parameters ρ 1 and ρ 2 (ρ 1 , ρ 2 ∈ (0, 1)) are envisioned. The resulting ARD1 processes are denoted by

( Y (i) t
)

t≥0
for i = 1, 2, respectively.

Proposition 1. We have:

1. Y nT decreases with respect to ρ in the sense of the likelihood order: If

ρ 1 < ρ 2 , then Y (2) nT ≺ lr Y (1)
nT ; 2. Y t decreases with respect to ρ in the sense of both increasing convex and concave orders:

If ρ 1 < ρ 2 , then Y (2) t ≺ icx Y (1) t and Y (2) t ≺ icv Y (1) t .
Proof. From (5), we know that

Y (i) nT ∼ Γ ( A (nT ) , b 1-ρ i ) for i = 1, 2. If ρ 1 < ρ 2 , we have b 1-ρ 1 ≤ b 1-ρ 2 and the first result follows from point 1 of Lemma 1, which entails that Y (2) nT ≺ icx Y (1) nT and Y (2) nT ≺ icv Y (1) nT . Let nT ≤ t < (n + 1) T . Based on (4), we have Y (2) t -Y (2) nT = Y (1) t -Y (1) nT = X (n+1) t -X (n+1) nT so that Y (2) t -Y (2) nT ≺ icx Y (1) t -Y (1) nT . We derive that Y (2) t = Y (2) nT + ( Y (2) t -Y (2) nT ) ≺ icx Y (1) t = Y (1) nT + ( Y (1) t -Y (1) nT )
because the icx order is stable through convolution (Shaked and Shanthikumar, 2007, Thm 4.A.8 (d) page 186). The proof is similar for ≺ icv .

Remark 2. Note that ≺ lr is not stable under convolution in a general setting so that the second point of the previous proposition would not be valid for this order (except if

X (n+1) t -X (n+1) nT
has a log-concave density (Shaked and Shanthikumar, 2007, Thm 1.C.9. page 46), namely if A (t) -A (nT ) ≥ 1 for t ∈ [nT, (n + 1)T )). A counter-example showing that the second point does not hold for the likelihood order is provided in Remark 4 later on.

Second model: Arithmetic Reduction of (virtual) Age of order 1 (ARA1)

The ARA1 model is based on the notion of virtual age as introduced by [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] in the context of recurrent events, which we first recall. Let us consider a system with successive lifetimes U 1 , U 2 , . . . , U n , . . . and instantaneous repairs at times

T n = ∑ n i=1 U i , n = 1, 2, . . . (with T 0 = 0). Let F be the survival function of U 1 .
Assume that there exists a sequence (V n ) n≥1 of non negative random variables such that after the n -th maintenance action at time T n , the next lifetime U n+1 has the same conditional distribution given V n = y as the remaining lifetime of a new system at time y:

P(U n+1 > t|V n = y) = P(U 1 -y > t|U 1 > y) = F (t + y) F (y) (7)
for all t, y > 0. Then, V n is called the virtual age of the system at time T n . After a maintenance action at time T n , the next lifetime U n+1 has the same distribution as if the calendar age of the system were equal to V n . Between repairs, the virtual age evolves with speed 1, just as the calendar time, so that for

T n ≤ t < T n+1 , the virtual age is V (t) = V n + (t -T n ).
In the context of the present paper, we use a similar notion of virtual age for deteriorating systems. To be more specific, considering a system with intrinsic deterioration modeled by (X t ) t≥0 and instantaneous repairs at times

T n , n = 1, 2, . . . (with T 0 = 0), we say that V (t) = V n +(t-T n ) stands
for the virtual age of the system at time t ∈ [T n , T n+1 ) if, given V n , the deterioration level of the maintained system is conditionally identically distributed as (X Vn+t-Tn ) Tn≤t<T n+1 (with independence between the deterioration and the virtual age, details further). The idea is just the same as for recurrent events: After a maintenance action at time T n , the system behaves just as if its calendar age were equal to V n at time T n and between maintenance actions, the virtual age evolves with speed 1. We now come to the specific virtual age model developed in this paper, which is called Arithmetic Reduction of Age of order 1 (ARA1) model after [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF]. It is based on the Kijima II imperfect repair model [START_REF] Kijima | Some results for repairable systems with general repair[END_REF] and each (periodic) repair removes the ρ% of the age accumulated by the system since the last maintenance action (or from the origin).

Let (Z t ) t≥0 be the process that describes the degradation level of the maintained system under this model of repair. At the first maintenance time T , the virtual age of the system is reduced by ρT units and it becomes V (T ) = (1-ρ)T . Recalling that (X (i) t ) t≥0 models the deterioration between the i-th and i + 1-th repairs, the degradation level on [0, T ] is given by

Z t = X (1) t for t < T, Z T = X (1)
(1-ρ)T . It means that, at time T , the system goes back into its past: The system is rejuvenated (from the age T to the age (1 -ρ)T ) and the deterioration level is reduced (from X

(1)

T to X (1) (1-ρ)T ). In case of A(•) convex, the deterioration rate is also reduced (from A ′ (T ) to A ′ ((1 -ρ)T )).
For T ≤ t < 2T , the system age is

V (t) = V (T ) + (t -T ) = t -ρT .
The corresponding deterioration level is identically distributed as X t-ρT . It is equal to the deterioration level at time T plus the increment of deterioration on (T, t], which leads to

Z t = Z T + ( X (2) t-ρT -X (2) (1-ρ)T ) ,
where X

(2)

t-ρT -X (2) (1-ρ)T is independent on Z T = X (1)
(1-ρ)T . At time 2T -(just before the repair), the age of the system is V (2T -) = 2T -ρT which is reduced by ρT at time 2T . The age hence is V (2T ) = 2(1ρ)T at time 2T . The corresponding deterioration level is given by:

Z 2T = Z T + ( X (2) 2(1-ρ)T -X (2) (1-ρ)T
) .

More generally, for nT ≤ t < (n + 1) T , the virtual age at time t is V (t) = t -ρnT (just as for an ARA1 model for recurrent events, see [START_REF] Doyen | Classes of imperfect repair models based on reduction of failure intensity or virtual age[END_REF])) and the system degradation is given by

Z t = Z nT + ( X (n+1) t-ρnT -X (n+1) (1-ρ)nT ) , ( 8
)
Z (n+1)T = Z nT + ( X (n+1) (1-ρ)(n+1)T -X (n+1) (1-ρ)nT ) (9)
where

X (n+1) t-ρnT -X (n+1) (1-ρ)nT is gamma distributed Γ (A (t -ρnT ) -A ((1 -ρ) nT ) , b) for all t ∈ [nT, (n + 1) T ].
Hence

Z nT = n ∑ i=1 ( X (i) (1-ρ)iT -X (i) (1-ρ)(i-1)T ) ( 10
)
and it is gamma distributed Γ

(A ((1 -ρ)nT ) , b).
Here, Z t is the sum of two independent gamma distributed r.v.s which share the same scale parameter b and it is gamma distributed Γ(A(t -ρnT ), b).

Also:

E (Z t ) = A (t -ρnT ) b , var (Z t ) = A (t -ρnT ) b 2 (11)
for all nT ≤ t < (n + 1)T .

Remark 3. Here again, E (Z t ) and var (Z t ) are decreasing with respect to ρ.

Just as for the ARD1 model, we next set

( Z (i) t ) t≥0
to be the ARA1 process with repair efficiency ρ i for i = 1, 2.

Proposition 2. Z t decreases with respect to ρ for the likelihood ratio order (and hence also for both increasing convex and concave orders): If

ρ 1 < ρ 2 , then Z (2) t ≺ lr Z (1) t .
Proof. Let ρ 1 < ρ 2 and nT ≤ t < (n + 1) T . For i = 1, 2, from (9), we have:

Z (i) t ∼ Γ (A (t -ρ i nT ) , b) , with A (t -ρ 2 nT ) ≤ A (t -ρ 1 nT ). We hence derive from Lemma 1 that Z (2) t ≺ lr Z (1) t .
Remark 4. We can see from Propositions 1 and 2 that as expected, the more efficient the maintenance action is (namely the larger ρ is), the smaller the deterioration level is for both ARD1 and ARA1 models. Based on (Shaked and Shanthikumar, 2007, Thm 1.C.5. p 44), the previous result implies for instance that, for an ARA1 model, we have

[Z (2) t |Z (2) t > h] ≤ sto [Z (1) t |Z (1) t > h] ( 12
)
for all h > 0. Imagine that h is an alert threshold in an application context and that the crossing of h triggers a signal, then, the previous relation means that, given that the signal has already been triggered, the deterioration level is stochastically all the smaller as the efficiency of the maintenance action is higher. Now considering A(t) = t 2 , b = 1, h = 0.5, t = 1.9, ρ 1 = 0.9 < ρ 2 = 0.95, we get that E(Y

(2) t |Y (2) t > h) ≃ 1.16 > E(Y (1) t |Y
(1) t > h) ≃ 1.08, which shows that ( 12) is not valid any more for the ARD1 model, in concordance with Remark 2. The stronger likelihood ratio result obtained for the ARA1 model may hence lead to different consequences from those for the ARD1 model in an application context (and in particular, conditional expectations are not necessarily ranked in an intuitive way).

Comparison of the moments

We now come to the main object of the paper, which is the comparison between the two models of imperfect repairs. Note that, in an application context, there is no reason why the estimated repair efficiency should be the same when the impact of the maintenance is modeled by an ARD1 or ARA1 model. Our point hence is to compare Y

(1) t and Z

(2) t , with efficiency ρ 1 and ρ 2 , respectively. In this section, we focus on the comparison of their respective means and variances.

Proposition 3. Let us consider the two following assertions:

E ( Z (2) t ) ≥ E ( Y (1) t
) for all t, T > 0 (13)

and A ((1 -ρ 2 )t) ≥ (1 -ρ 1 )A(t) for all t > 0. ( 14)

Then:

1. Assertion (13) implies assertion (14). 2. If A(•) is concave, then the converse is also true (namely (14) =⇒ ( 13)). 3. As a special case, if A(•) is concave and ρ 2 ≤ ρ 1 , then (13) is true.

All the previous results are valid with reversed inequalities and concave substituted by convex.

Proof. Based on ( 6) and ( 11), assertion ( 13) is equivalent to

A(t -ρ 2 nT ) ≥ A(t) -ρ 1 A(nT ) ( 15
)
for all T > 0, all n ∈ N and all nT ≤ t < (n + 1) T . Taking t = nT in (15), it implies ( 14) for all t = nN with T > 0 and n ∈ N, and hence for all t > 0. This shows the first point.

For the second point, assume A(•) to be concave. Then A (t -ρ 2 nT ) -A(t) is non decreasing with respect to t and for nT ≤ t < (n + 1)T , we get that

A(t -ρ 2 nT ) -A(t) ≥ A((1 -ρ 2 )nT ) -A(nT ). ( 16
)
Assuming ( 14) to be true, we easily derive that

A(t -ρ 2 nT ) -A(t) ≥ (1 -ρ 1 )A(nT ) -A(nT ) = -ρ 1 A (nT )
so that ( 15) is true. This implies (13) and the second point is proved.

In the specific case where A(•) is concave and ρ 2 ≤ ρ 1 , we have

A((1 -ρ 2 )nT ) = A ((1 -ρ 2 )nT + ρ 2 0) ≥ (1 -ρ 2 )A(nT ) (17) ≥ (1 -ρ 1 )A(nT )
for all n and T (using A(0) = 0 in the first line). This implies that ( 14) is true. Point three now is a direct consequence of point two.

The reasoning is similar for reversed inequalities and it is omitted.

Remark 5. Note that Condition ( 14) is less restrictive than ρ 2 ≤ ρ 1 in the concave case. (The same with a reverse inequality in the convex case). For instance, consider

A (t) = t β with 0 < β < 1. Then Condition (14) means that (1 -ρ 2 ) β ≥ 1 -ρ 1 , which is less restrictive than ρ 1 ≥ ρ 2 .
In the following example, we look at the comparison of the expectations when the conditions in Proposition 3 are not fulfilled.

Example 1. As a first case, we take A (t) = e t -1 (convex function), b = 1, T = 1 ρ 1 = 0.95, ρ 2 = 0.5. As a second case, we take A (t) = 1 -e -t (concave function), b = 1, T = 1 ρ 1 = 0.5, ρ 2 = 0.95. Neither condition (14) nor the reversed inequality is true on the whole real line. The corresponding expectations are plotted in Figures 1a (first case) and 1b (second case). We can see that the respective means of Y

(1) t and Z

(2) t are not ordered in the same way on the whole real line.

Proposition 4. Let us consider the two following assertions:

V ar

( Z (2) t ) ≥ V ar ( Y (1) t
) for all t, T > 0 (18)

and A ((1 -ρ 2 )t) ≥ (1 -ρ 1 ) 2 A(t) for all t > 0. ( 19
)
Then:

1. Assertion (18) implies assertion (19).

1 2 3 4 t Y (1) t Z (2) t (a) A (t) = e t -1, ρ 1 = 0.95, ρ 2 = 0.5 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 t Y (1) t Z (2) t (b) A (t) = 1 -e -t , ρ 1 = 0.5, ρ 2 = 0.95 Figure 1: Expectations of Y (1) t and Z (2) t with b = 1, T = 1 2. If A(•) is concave, then the converse is also true (namely (19) =⇒ (18)). 3. As a special case, if A(•) is concave and 1 -ρ 2 ≥ (1 -ρ 1 ) 2 , then (18)
is true.

All the previous results are valid with reversed inequalities, and concave substituted by convex.

Proof. Based on ( 6) and (11), inequality ( 18) is equivalent to

A(t -ρ 2 nT ) ≥ A(t) -ρ 1 (2 -ρ 1 ) A(nT ) (20)
for all T > 0, all n ∈ N, all nT ≤ t < (n + 1) T . Considering t = nT , it implies (19), so that point one is true. Now assume A(•) to be concave. Based on (16), we have:

A(t -ρ 2 nT ) -A(t) + ρ 1 (2 -ρ 1 )A(nT ) ≥ A((1 -ρ 2 )nT ) -(1 -ρ 1 ) 2 A(nT ).
If ( 19) is true, then ( 20) is consequently true and ( 18) too. This shows the second point.

In the specific case where

A(•) is concave and (1 -ρ 2 ) ≥ (1 -ρ 1 ) 2 , we have A((1 -ρ 2 )nT ) ≥ (1 -ρ 2 )A(nT ) ≥ (1 -ρ 1 ) 2 A(nT )
(based on (17) for the first inequality). Hence (20) is true, so that ( 18) is true too. This provides point three.

The reasoning is similar for reversed inequalities and it is omitted.

0 1 2 3 4 t Y (1) t Z (2) t (a) A (t) = e t -1, ρ 1 = 0.95, ρ 2 = 0.5 0 1 2 3 4 5 6 0 0.2 0.4 0.6 0.8 t Y (1) t Z (2) t (b) A (t) = 1 -e -t , ρ 1 = 0.5, ρ 2 = 0.95

Example 2. We now consider the same data as for Example 1, where neither condition (19) nor the reversed inequality is true on the whole real line. The corresponding variances are plotted in Figures 2a (first case) and 2b (second case). We can see that the respective variances of Y

(1) t and Z

(2) t are not ordered in the same way on the whole real line.

Finally, we easily derive from Propositions 3 and 4 the following corollary, where the expectation and variance of the two processes are compared assuming a power law shape function for the gamma process.

Corollary 1. We consider A(t) = αt β with α, β > 0. Then we get

1. If β ≤ (≥)1, then E ( Z (2) t ) ≥ (≤) E ( Y (1) t ) ∀t ⇔ (1 -ρ 2 ) β ≥ (≤) (1 -ρ 1 ) . 2. If β ≤ (≥)1, then V ar ( Z (2) t ) ≥ (≤)V ar ( Y (1) t ) ∀t ⇔ (1 -ρ 2 ) β ≥ (≤) (1 -ρ 1 ) 2 .

Stochastic comparison between Y

(1) t and Z

(2) t

We now come to the stochastic comparison between Y

(1) t and Z

(2) t , as given by ( 4) and ( 9), with ρ substituted by ρ i , i = 1, 2, respectively.

Proposition 5. If A ((1 -ρ 2 )nT ) ≥ (≤) (1 -ρ 1 )A(nT ), ( 21
)
then Y

(1)

nT ≺ icx (≻ icv ) Z (2)
nT .

Proof. From ( 5) and ( 10), we know that Y

(1)

nT and Z

(2)

nT are gamma distributed, with distributions Γ

( A (nT ) , b 1-ρ 1 )
and Γ(A ((1 -ρ 2 )nT ) , b), respectively. The results can hence be obtained by application of Lemma 1.

We next focus on the comparison between

Y (1) t -Y (1) nT and Z (2) t -Z (2) nT . Proposition 6. If A (t) is convex (concave), then for all nT ≤ t < (n + 1) T : • Z (2) t -Z (2) nT ≺ lr (≻ lr ) Y (1) t -Y (1) nT , • Y (1) t -Y (1) nT ≺ lc (≻ lc ) Z (2) t -Z (2) nT . Proof. For nt ≤ t < (n + 1) T , let us set Y = Y (1) t -Y (1) nT = X (n+1) t -X (n+1) nT (see (4)) and Z = Z (2) t -Z (2) nT = X (n+1) t-ρ 2 nT -X (n+1) (1-ρ 2 )nT
(see ( 9)). Then Y and Z are gamma distributed, with distributions Γ(A(t)-A(nT )) and Γ(A (t -ρ 2 nT ) -A ((1 -ρ 2 ) nT )), respectively. Now, the likelihood ratio comparison result is a direct consequence of Lemma 1, because the convexity (concavity) of A(•) entails that

A(t) -A(nT ) ≥ (≤)A (t -ρ 2 nT ) -A ((1 -ρ 2 ) nT ) .
As for the log-concave order, using the notations of Definition 2, we have ( log

( f Y f Z )) ′ (y) = [A (t) -A (nT )] -[A (t -ρ 2 nT ) -A ((1 -ρ 2 ) nT )] y .
This function decreases (increases) when A (t) is convex (concave), which provides the result.

(f Y (1) t /f Z (2) t ) for α = 1.25 Figure 3: Plots of log(f Y (1) t /f Z (2) t ) for A(t) = t α
Hence, if A (t) is convex (concave), the increment between times nT and t (with nT ≤ t < (n + 1) T ) is smaller (larger) for the ARA1 model than for the ARD1 model in the sense of the likelihood ratio order (and consequently also for the usual stochastic and increasing convex/concave orders), but it has a larger (smaller) variability for the ARA1 model than for the ARD1 model in the sense of the log-concave order.

Based on the previous results, if A (•) is concave, we have Z

(2)

t -Z (2) nT ≺ lc Y (1) t -Y (1) nT but Y (1) nT ≺ lc Z (2)
nT (for nT ≤ t < (n + 1) T ). There consequently is no real hope that Y

(1) t

≺ lc Z (2) t . When A (•) is convex, we have both Y (1) t -Y (1) nT ≺ lc Z (2) t -Z (2) nT and Y (1) nT ≺ lc Z (2)
nT , and Y

(1) t ≺ lc Z

(2) t might be valid. However, remembering that the log-concave order is not closed under convolution, see [START_REF] Whitt | Uniform conditional variability ordering[END_REF], the question deserves to be further studied, which is done in the following example.

Example 3. The function log

( f Y (1) t /f Z (2) t )
is plotted in Figures 3a and3b for ρ 1 = 0.5, ρ 2 = 0.4, T = 1, A (t) = t α , n = 10, b = 1 at time t = 10.2, with α = 0.75 and α = 1.25, respectively. We observe that we do not have Y

(1) t ≺ lc Z

(2) t neither for α = 0.75 (as expected) nor for α = 1.25 for which the question was open. As a conclusion, in a general setting, Y

(1) t and Z

(2) t are not comparable with respect to the log-concave order.

Next result provides some conditions under which Y

(1) t and Z

(2) t are comparable with respect to either the increasing convex or concave order.

Theorem 1. If A (•) is concave and A ((1 -ρ 2 )t) ≥ (1 -ρ 1 )A(t) for all t > 0 ( 22
)
(condition of Proposition 3, known to be true if

ρ 1 ≥ ρ 2 ), then Y (1) t ≺ icx Z (2) t
for all t ≥ 0 (and all T ≥ 0).

If A (t) is convex with a reversed inequality in (22), then Z (2) t ≺ icv Y
(1) t for all t ≥ 0 (and all T ≥ 0).

Proof. Writing Y (1) t = ( Y (1) t -Y (1) nT ) + Y (1)
nT (the same for Z

(2) t ), the results are direct consequences from Propositions 5 and 6, based on the fact that both increasing convex and concave orders are stable under convolution (Shaked and Shanthikumar, 2007, Thm 4.A.8 (d)

page 186).

Example 4 ( Increasing convex order). We consider t = 10.5, T = 1, b = 1, ρ 1 = ρ 2 = 0.75 with A(t) = t 0.7 as a first case (concave case with condition (22) fulfilled) and A(t) = t 1.1 as a second case (convex case with reversed condition (22) fulfilled). The difference (see ( 1)). We observe that it changes sign in the convex case, which shows that Z

∫ +∞ x FZ (2) t (u) du - ∫ +∞ x FY (
(2) t and Y

(1) t

are not comparable with respect to the increasing convex order.

Example 5 ( Increasing concave order). We consider t = 10.5, T = 1, ρ 1 = 0.8, ρ 2 = 0.78, b = 1 with A(t) = t 0.9 (concave case with condition (22) fulfilled) and A(t) = t 1.1 (convex case with reversed condition (22) true). The difference ) and 5b (second case). We observe that, as expected, Z

∫ x 0 F Z (2) t (u)du - ∫ x 0 F Y (1) t (u)du is plotted in Figures 5a (first
(2) t ≺ icv Y

(1) t (see ( 2)) in the second case whereas Z

(2) t and Y

(1) t

are not comparable with respect to the increasing concave order in the first case.

Finally, we end this section by considering the case of a homogeneous gamma process (A(t) = at).

Corollary 2. Assume that A(t)

= at for all t ≥ 0, where a > 0. We have the following results:

(b) A(t) = t 1.1 Figure 4: Plots of ∫ ∞ x FZ (2) t (u)du - ∫ ∞ x FY (1) t (u)du for A(t) = t α x 0 2 4 6 8 -0.6 -0.4 -0.2 0 0.1 (a) A(t) = t 0.9 x 0 2 4 6 8 0 0.2 0.4 (b) A(t) = t 1.1 Figure 5: Plots of ∫ x 0 F Z (2) t (u)du - ∫ x 0 F Y (1) t (u)du for A(t) = t α 1. If ρ 1 ≥ ρ 2 , then Y (1) t ≺ icx Z (2) t
and hence E(Y

(1)

t ) ≤ E(Z (2) t ) for all t ≥ 0; 2. If ρ 1 ≤ ρ 2 , then Z (2) t ≺ icv Y (1) t and hence E(Z (2) t ) ≤ E(Y (1) t ) for all t ≥ 0; 3. If ρ 1 = ρ 2 , then Y (1) t ≺ cx Z (2) t and Z (2) t ≺ cv Y (1) t and hence E(Z (2) t ) = E(Y (1) t ) for all t ≥ 0; 4. V ar ( Z (2) t ) ≥ V ar ( Y (1) t ) for all t > 0 if and only if 1-ρ 2 ≥ (1-ρ 1 ) 2 .
Proof. Points 1, 2 and 4 are direct consequences of Theorem 1 and Corollary 1. As for point 3, assume that

ρ 1 = ρ 2 . This entails that E ( Y (1) t ) = E ( Z (2) t
)

. Based on the first two points, the result follows from (Shaked and Shanthikumar, 2007, Theorem 4.A.35 page 197).

Mostly equivalent imperfect repair models

In an applied context, parameters ρ 1 and ρ 2 for ARD1 and ARA1 models, respectively, will be estimated from feedback data, which will be typically gathered at maintenance times iT , i ≥ 1. As a consequence, we can expect that the estimated parameters ρ1 and ρ2 should be such that the corresponding expected deterioration levels should be very similar at maintenance times, namely such that

E ( Y (1) iT ) = (1 -ρ1 ) A (iT ) b ≃ E ( Z (2) iT ) = A ((1 -ρ2 ) iT ) b .
There hence is a specific interest for the applications to compare the ARD1 and ARA1 models under the condition

(1 -ρ 1 ) A (iT ) = A ((1 -ρ 2 ) iT ) for i ≥ 1, (23)
on (ρ 1 , ρ 2 ), which will lead to mostly equivalent deterioration levels (at least at maintenance times). However, the previous requirement (23) does not seem to have a solution for a general shape function A (•). We hence restrict the study to the power law case

A (t) = αt β (with α, β > 0), for which (23) is just equivalent to 1 -ρ 1 = (1 -ρ 2 ) β . ( 24
)
This section is hence devoted to this specific power law case with the previous relationship between ρ 1 and ρ 2 , which ensures that

E ( Y (1) iT ) = E ( Z (2) iT
) for all i ≥ 1.

Remark 6. This specific "equivalent" case has a similar spirit to that detailed in (Doyen and Gaudoin, 2004, Property 4), where the authors match the minimal wear intensities of two imperfect repair models for recurrent events, based on the reduction of either virtual age or failure intensity.

In case of a homogeneous gamma process (β = 1), the equivalent case corresponds to identical repair efficiencies for both ARD1 and ARA1 models (ρ 1 = ρ 2 ), which has already been studied in Corollary 2 (point 3). We now investigate the case of a general β.

In the equivalent case, there is equality in ( 14), ( 21) and ( 22). Based on the fact that A (•) is concave (convex) when β ≤ (≥) 1, we directly get the following results from Proposition 5, Theorem 1 and Corollary 1.

Corollary 3. Assume that A (t) = αt β (with α, β > 0) and that (ρ 1 , ρ 2 ) fulfills (24) (equivalent case). Then:

1. Z (2) nT ≺ cv Y (1) nT and Y (1) nT ≺ cx Z (2) nT (which both entail that var ( Z (2) nT ) ≥ var ( Y (1) nT )
, see (Shaked and Shanthikumar, 2007, (3.A.4) page 110)) for all n ≥ 1.

2. If β ≤ (≥) 1, then Y (1) t ≺ icx (≻ icv ) Z (2) t (which entails that E ( Y (1) t ) ≤ (≥) E ( Z (2) t
)

) for all t ≥ 0.

3.

If β ≤ 1, then V ar ( Z (2) t ) ≥ V ar ( Y (1) t
) for all t ≥ 0.

Remark 7. Based on the previous result, we can see that even if the two imperfect repair models provide similar expected deterioration levels at maintenance times, there are differences in both their location and spread between the repairs. For instance, considering β ≤ 1, a possible by-product of Y

(1) t

≺ icx Z (2) t is that E((Y (1) t -h) + ) ≤ E((Z (2)
t -h) + ) for any h > 0, see, e.g., (Shaked and Shanthikumar, 2007, (4)

p 182). If h is a critical deterioration level in an application context, (Y

(1) t -h) + and (Z (2) t -h) + correspond to the hazardous part of deterioration (beyond the critical level) and this means that the expected "risk" is lower for the ARD1 model than for the ARA1 one.

Note that Corollary 3 does not provide any insight for the comparison of the variances at time t when β > 1, which hence deserves a different analysis on which we now focus.

Proposition 7. Let β > 1 and let

g (x) = x β -(x -ρ 2 ) β - ( 1 -(1 -ρ 2 ) 2β ) (25) for x ∈ [1, 2). • If g (2 -) ≤ 0, then V ar ( Y (1) t ) ≤ V ar ( Z (2) t ) ( 26
)
for all t ≥ 0.

• If g (2 -) > 0, there exists one single x * ∈ (1, 2) such that g (x * ) = 0. Also:

-

Inequality (26) is true for all t ≥ t * = ⌈ 1 x * -1 ⌉
T , where ⌈•⌉ stands for the ceiling function;

-For each n < 1 x * -1 , inequality (26) is true for all t ∈ [nT, x * nT ), with a reversed inequality for t ∈ [x * nT, (n + 1) T ).
Proof. For nT ≤ t < (n + 1)T , it is easy to check that ( 26) is equivalent to

t β - ( 1 -(1 -ρ 2 ) 2β ) (nT ) β ≤ (t -ρ 2 nT ) β ,
which can also be written as g ( t nT

) ≤ 0 where 1 ≤ t nT ≤ 1 + 1 n and where g is defined by ( 25

) for x ∈ ∪ n≥1 [1, 1 + 1 n ) = [1, 2
). As β > 1, the function g increases from g (1) to g (2 -) and based on the fact that var

( Z (2) nT ) ≥ var ( Y (1) nT
) (first point of Corollary 3), we have g (1) ≤ 0. As for the sign of

g ( 2 -) = 2 β -(2 -ρ 2 ) β - ( 1 -(1 -ρ 2 ) 2β
) , there are two possibilities, which lead to the following cases:

• If g (2 -) ≤ 0, then g (x) ≤ 0 for all x ∈ [1, 2) and ( 26) is true for all t ≥ 0.

• If g (2 -) > 0, there exists one single x * ∈ (1, 2) such that g (x * ) = 0, with g (x) < 0 for all x ∈ [1, x * ) and g (x) > 0 for all x ∈ (x * , 2).

- 26) is true for all t ∈ [nT, (n + 1) T ). This inequality is hence true for all t ≥ n * T with n * = ⌈ 1 x * -1 ⌉ .

If 1 + 1 n ≤ x * (namely n ≥ 1 x * -1 ), then g (x) < 0 for all x ∈ [1, 1 + 1 n ) and (
-If 1 + 1 n > x * (namely n < 1 x * -1 ), then g (x) < 0 for all x ∈ [1, x * ) and ( 26) is true for all t ∈ [nT, x * nT ), with a reversed inequality for t ∈ [x * nT, (n + 1) T ).

Remark 8. In the case where g (2 -) > 0, note that the inequality n < 1

x * -1 is always valid for n = 1 so that for small n, the difference V ar

( Y (1) t ) - V ar ( Z (2) t
) will always cross 0 (from -to +) on [nT, (n + 1) T ) (and will remain negative for larger n).

Maintenance strategies

The reward function

This section is devoted to the analysis of maintenance strategies considering the two types of repair. In all the section, the system is assumed to provide a reward which decreases when the deterioration level of the system increases. Based on classical functions used in the insurance literature [START_REF] Rolski | Stochastic Processes for Insurance and Finance[END_REF]), we assume that the reward function is given by

g(x) = (b 1 -k 1 e α 1 x ) 1 {0≤x≤c} + (b 2 -k 2 e α 2 x ) 1 {c<x} , ( 27
)
with b 1 , b 2 , α 1 , α 2 , k 1 , k 2 , c > 0 and x ≥ 0, where g(x) stands for the unitary reward per unit time when the degradation level of the system is x. The function g is supposed to be continuous and positive on [0, c),

which implies that b 2 -k 2 e α 2 c = b 1 -k 1 e α 1 c > 0. ( 28
)
Also, we assume that α 1 ≤ α 2 and k 1 ≤ k 2 so that level c appears as a critical level, from which the system becomes less performing. With the previous assumptions, it is easy to check that g is a concave function and that g(x) > 0 if and only if x < L = ln(b 2 /k 2 ) α 2

. Level L hence appears as a critical threshold.

An example of reward function is plotted in Figure 6 (28). With this dataset L ≃ 10.0144. Two different maintenance strategies are envisioned in the two following subsections. In each case, the system is put into operation at time t = 0 and it degrades according to a non-homogeneous gamma process with parameters A(t) and β. For each maintenance strategy, the two imperfect repair models are envisioned (ARA1 or ARD1) and the comparison between the two types of repair is performed through their corresponding expected reward (profit) rates per unit time on a long time run.

(n, T ) policy

Starting from n ∈ N * and T > 0, the (n, T ) maintenance scheme is developed as follows:

• Imperfect repairs based on either one of the two models (ARD1 or ARA1) are performed at times T, 2T, 3T, . . .

• Each imperfect repair costs C r monetary units (m.u.),

• The profit per unit time is given by the reward function g from ( 27),

• The system is replaced by a new one at the time of the n-th imperfect repair (nT ) with a cost of C m.u..

Based on the renewal reward theorem, see, e.g., [START_REF] Tijms | A First Course in Stochastic Models[END_REF], the long time reward rate per unit time for this policy is given by:

R ARD (n, T ) = ∫ nT 0 E ( g(Y (1) s ) ) ds -(n -1)C r -C nT ( 29
)
when ARD1 repairs are considered and

R ARA (n, T ) = ∫ nT 0 E ( g(Z (2) s ) ) ds -(n -1)C r -C nT ( 30
)
for ARA1 repairs, where Y

(1) s and Z

(2) s are is given by ( 4) and ( 9), respectively, with ρ substituted by ρ i , i = 1, 2 .

Remark 9. Based on the fact that the function -g(•) is increasing and convex, the previously obtained theoretical results allow to derive several observations on the reward rates:

• From Propositions 1 and 2, we get that E

( -g(Y (1) t )
) and E

( -g(Z (2) t )
)

decrease with ρ 1 and ρ 2 for all t > 0, respectively. This entails that the objective profit rates R ARD (n, T ) and R ARA (n, T ) increase with the effectiveness of the repair (ρ 1 and ρ 2 , respectively) for both ARD1 and ARA1 models.

• If the shape function A(•) is concave and A((1 -ρ 2 )t) ≥ (1 -ρ 1 )A(t) for all t > 0 (Condition (22)), then Theorem 1 entails that

E ( g(Z (2) t ) ) ≤ E ( g(Y (1) t )
)

, ∀t > 0, from which we derive that the objective profit functions for the two repair models are comparable:

R ARA (n, T ) ≤ R ARD (n, T ) for all n, T.

We now come to some numerical illustrations.

Example 6. The parameters of the gamma process are A(t) = t 0.5 + t 0.75 (concave function) and b = 1; those for the reward reward function g are α Figure 7 shows the operating profit rate given in ( 29) for both ARD1 and ARA1 models using ρ 1 = ρ 2 = 0.5 (for which inequality ( 22) is true). The computations have been made using 8 points for T from 1 to 6 and 10 points for n from 1 to 10 with 5000 simulations in each point. The Simpson method is applied for the integrals in ( 29), with 20 points from 0 to nT .

1 = 0.1, α 2 = 0.25, k 1 = 1, k 2 = 1, b 1 = 11, c = 4,
Under an ARD1 model, the optimal profit rate is obtained for (n, T ) = (7, 2.43) with a profit rate of R ARD (7, 2.43) = 6.85 m.u. per unit time. Under an ARA1 model, the optimal profit rate is obtained for (n, T ) = (4, 3.14) with a profit rate of R ARA (4, 3.14) = 5.29 m.u. per unit time. Figure 8a shows the difference of the profit rates under the two repair models, that is, R ARD (n, T ) -R ARA (n, T ) for all n and T . As expected from the previous theoretical results, R ARA (n, T ) ≤ R ARD (n, T ). Also, we can observe that the difference between the two rates increases as n and T increases.

Example 7. Keeping the same parameters as in the previous example except from the repair efficiency for the ARD1 model which becomes ρ 1 = 0.31, Figure 8b shows the difference between the profit rates under the two repair models. Here A(t) = t 0.5 + t 0.75 is concave but condition ( 22) of Theorem 1 is not valid any more. We observe that there is no dominance of the profit rate of one model over the other.

Although the (n, T ) policy allows us to compare the two models of repair, under this maintenance policy, the system is (imperfectly) repaired even when the system is so degraded that the reward has become negative. We now suggest a more realistic condition-based maintenance strategy, where the maintenance action depends on the degradation level of the system.

(M,T) policy

Let T > 0 and let M ∈ [0, L) be a preventive maintenance threshold. (We recall that L is the critical threshold defined in Subsection 7.1 from where the reward becomes negative). The (M, T ) condition-based maintenance scheme is developed as follows:

• The system is inspected at times T, 2T, 3T, . . . and the system degradation level is checked.

• By an inspection:

-If the degradation level does not exceed the threshold M , an imperfect repair based on either one of the two models (ARD1 or ARA1) is performed with a cost of C r m.u.;

-If the degradation level is between levels M and L, a preventive replacement is performed and the system is instantaneously replaced by a new one with a cost of C p m.u.;

-If the degradation level exceeds L, a instantaneous corrective replacement takes place with a cost of C c m.u..

• The profit per unit time is given by the reward function g from ( 27).

The successive (corrective or preventive) replacements of the system appear as the points of a renewal process, and the long time profit rate per unit time is given by

C ARD (T, M ) (31) = E ( ∫ R 0 g ( Y (1) s ) ds -C r ([R/T ] -1) -C p 1 { M ≤Y (1) R <L } -C c 1 { L≤Y (1) R } ) E(R)
for the ARD1 model, with a similar expression for the ARA1 model (C ARA (T, M )), where R stands for the time to a system replacement and g denotes the reward function given by ( 6). Due to the complexity of the (M, T ) policy, there is no hope here to find analytical conditions that could ensure the dominance of one function C ARD (T, M ) or C ARA (T, M ) over the other.

Their comparison is hence made on a numerical example. Figure 9a shows the profit rate for the maintained system under an ARD1 repair. The optimal maintenance strategy is obtained for (M, T ) = (9.21, 3.05) with a profit rate of C ARD (3.05, 9.21) = 673.94 m.u. per unit time. Figure 9b shows the profit rate for the ARA1 repairs. The optimal maintenance strategy is obtained for (M, T ) = (10.24, 3.05) with a profit rate of C ARA (3.05, 10.24) = 684.34 m.u. per unit time. These figures have been computed considering a grid of 10 points for T from 1.14 to 4 and a grid of 13 points for M from 1 to L and 10000 simulations for each pair of points.

Figure 10 shows the difference between the profit rates C ARD (T, M ) -C ARA (T, M ) for this dataset. Although conditions of Theorem 1 are fulfilled, we can see that the sign of the difference changes, so that there is no dominance of the profit rate of one model on the other.

Conclusions and perspective

Two imperfect repair models for a degrading system are compared in this paper. The comparison is performed in terms of location and spread of the two resulting stochastic processes. Results are provided in terms of moments and likelihood ratio ordering, as well as in terms of (increasing) convex/concave ordering, which is not so common in the reliability literature. Two maintenance strategies are also developed, which are assessed through a reward function, which takes into account the effective deterioration level of the system (the lower the deterioration level, the higher the reward), contrary to classical objective functions from the literature.

The paper is developed under a periodic imperfect repair scheme, for sake of simplification. It is however easy to check that all the results of the paper would remain valid under a deterministic non periodic repair scheme (with the same maintenance times for both imperfect repair models), with very slight modifications. Even more, considering random maintenance times (independent on the deterioration level and identically distributed for both imperfect repair models), most results would also remain valid, such as the likelihood ratio and (increasing) convex/concave comparison results, based on the closure under mixture property of these stochastic orders (Shaked and Shanthikumar, 2007, Thm 1.C.15. p 48, Thm 4.A.8. p 185).

Note also that if the paper focuses on some specific stochastic orders, other ones could also be considered such as Laplace transform or Excess Wealth orders for instance. Other questions of interest concern the comparison of remaining lifetimes, considering the system as failed (or too degraded) when its deterioration level is beyond a fixed failure (critical) threshold. From a theoretical point of view, this seems a difficult issue in a general setting. One could then look at partial results in specific situations.

Figure

  Figure 2: Variances of Y (1) t and Z (2) t with b = 1, T = 1

  case

Figure 6 :

 6 Figure 6: Reward function for α1 = 0.1, b1 = 11, α2 = 0.25, k1 = 1, k2 = 1, c = 4

  which implies that b 2 = 12.2265. The (n, T ) policy is considered with C r = 2 m.u. as cost of imperfect repair and C = 25 m.u. as replacement cost.

Figure 7 :

 7 Figure 7: Operational profit rate, (n, T ) policy

  Case ρ 1 = 0.31 and ρ 2 = 0.5

Figure 8 :

 8 Figure 8: Differences between the operational profit rates, (n, T ) policy

Example 8 .

 8 The parameters of the gamma process are A(t) = 1.3t and b = 0.8. For the reward function g, they areα 1 = 0.4, α 2 = 0.5, b 1 = 800, k 1 = 1.05, k 2 = 1.07, c = 8,which implies b 2 = 832.6609 and L = 13.3139. The repair efficiencies of the ARD1/ARA1 repairs are ρ 1 = ρ 2 = 0.9. Their common cost is C r = 200 m.u.. The cost of a preventive replacement is C p = 1000 m.u. whereas it is C c = 1300 m.u. for a corrective one.

  ρ 2 = 0.9, ARA1 model

Figure 9 :

 9 Figure 9: Operational profit rates, (M, T ) policy

Acknowledgements

Both authors warmly thank the three referees and the Editor for their constructive remarks and comments, which have led to a much clearer and better structured paper.

Another point of interest would be to try and compare the two types of imperfect repairs dropping the gamma-process assumption. Based on the fact that, under technical conditions, normal random variables are comparable with respect to several stochastic orders (see, e.g., [START_REF] Müller | Comparison methods for stochastic models and risks[END_REF]), one can wonder whether it could be possible to get some similar results as in the paper for Wiener processes with drift. (Not all however, because comparison results between normal random variables in the likelihood ratio ordering sense require that the random variables share the same variance, which cannot be the case in our context. The same for the logconcave order, which requires that the two normal random variables share the same mean, see [START_REF] Whitt | Uniform conditional variability ordering[END_REF]). Other deteriorations processes might also be envisioned, such as inverse Gaussian or inverse Gamma processes.

Finally, other maintenance strategies could be envisioned, based on either one of the two types of imperfect repair.


Hafiz Tiomoko Ali

Romain Couillet
email: romain.couillet@centralesupelec.fr

Random matrix improved community detection in heterogeneous networks

This article proposes a new spectral method for community detection in large dense networks following the degree-corrected stochastic block model. We theoretically support and analyze an approach based on a novel "α-regularization" of the modularity matrix. We provide a consistent estimator for the choice of α inducing the most favorable community detection in worst case scenarios. We further prove that spectral clustering ought to be performed on a 1 -α regularization of the dominant eigenvectors (rather than on the eigenvectors themselves) to compensate for biases due to degree heterogeneity. Although focused on dense graph models, our clustering method is shown to be very promising on real world networks with competitive performances versus state-of-the-art spectral techniques developed for sparse homogeneous networks.

I. INTRODUCTION

One of the prominent features of real world networks is their community structure which subdivides the network graphs in disjoint clusters: nodes are densely connected inside each cluster and sparsely connected across clusters. Extracting these clusters from the observation of the graph and mapping the nodes to their respective clusters allows for a better understanding and analysis of large complex networks and is one of the challenging tasks in modern network mining. There are two major methods in the literature to perform this task: statistical inference methods which rely on some knowledge of the parameters of an underlying graph generative model and spectral algorithms which classify vertices using the eigenvectors corresponding to outlying eigenvalues of some matrix representation of the network [START_REF] Fortunato | Community detection in graphs[END_REF].

Large real world networks are usually sparse in the sense that the average degree of the nodes remains constant with respect to the network size n as n → ∞ and have heterogeneous degree distributions, often modelled as power laws [START_REF] Fortunato | Community detection in graphs[END_REF]. Standard spectral algorithms based on variants of the adjacency matrix have been proved suboptimal in sparse network community detection where they are supplanted by more powerful operators arising from statistical physics, notably the recently proposed non-backtracking [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] or Bethe Hessian (BH) matrix approaches [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. The latter have however been developed under a stochastic block model (SBM) assumption which does not allow for degree heterogeneity inside clusters and tend to fail in the presence of highly varying node degrees in the network. To allow for degree heterogeneity modeling, in this article we instead consider the so-called degree-corrected stochastic block model (DCSBM), initially proposed in [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF]. Denoting G a K-class graph of n vertices with communities

Modularity

Bethe Hessian

Algo with α = 1 Algo with αopt Fig. 1. Two dominant eigenvectors (x-y axes) for n = 2000,

K = 3 classes C 1 , C 2 and C 3 of sizes |C 1 | = |C 2 | = n 4 , |C 3 | = n 2 , intrinsic probabilities taking two values q 1 = 0.1, q 2 = 0.5, matrix of weights C = 1 3 1 T 3 + 100 √ n I 3 .
C 1 , . . . , C K and q i , 1 ≤ i ≤ n, an intrinsic probability for node i to connect to any other network node, the DCSBM assumes an adjacency matrix A ∈ {0, 1} n×n with A ij independent Bernoulli random variables with parameter P ij = q i q j C ab , for i ∈ C a and j ∈ C b , where C ab is a community-related correction factor (in general C aa > C ab for all a = b).

To illustrate the aforementioned limitations of spectral methods under the DCSBM model, the top two graphs of Figure 1 provide 2D representations of dominant eigenvector 1 versus eigenvector 2 for the standard modularity matrix and the sparsity-improved BH matrix, when half the nodes connect with low probability q 1 and half the nodes with high probability q 2 . For both methods, this erroneously induces the detection of extra communities and even a confusion of genuine communities in the BH approach. The bottom graphs of Figure 1 are the results of the "α-regularized" modularity approach proposed in this article, for an arbitrary and an optimized values of the parameter α, which both restore (with different levels of accuracy) the genuine community structure.

In addition to studying realistic degree corrected graph models, our investigation is motivated by the lack of unification of the spectral methods for community detection in the literature, based on various normalizations of the adjacency matrix: the adjacency matrix itself A (A ij = 1 if nodes i and j are connected, else 0), the modularity matrix B = A-dd T d T 1n , with d = A1 n the vector of degrees, the unnormalized Laplacian L = D -A, with D = diag(d), the normalized Laplacian L norm = I n -D -1/2 AD -1/2 , the random walk matrix D -1 A, etc. [START_REF] Fortunato | Community detection in graphs[END_REF]. In this article, we theoretically support the importance to study a generalized version L α ∝ D -α BD -α , with α ∈ [0, 1]. In particular, L 0 ∝ B and L 1/2 ∝ I n -L norm .

Our key objective is thus to understand the eigenstructure of L α and its consequence to community detection.

Our first main result is to show that L α is asymptotically well approximated by a spiked random matrix Lα , more amenable to analysis than L α itself. Spiked random matrices [START_REF] Benaych | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF] are low rank perturbations of standard random matrices and their spectrum is generally composed of one or many "bulks" of eigenvalues and, whenever a phase transition is met, of additional eigenvalues which isolate from these bulks. In this case, the eigenvectors associated with the isolated eigenvalues are correlated to the low rank matrix eigenvectors. In our context, the approximation Lα of L α tells us that: (i) the phase transition, which corresponds to the community detectability threshold, depends heavily on α and (ii) the dominant eigenvectors of L α are biased by the degree distribution (unless α = 0) and thus need to be properly normalized. A deeper analysis reveals the existence of an optimal choice α opt of α for which the detectability threshold is smallest, thus allowing for the weakest detectable community structure. We finally provide a consistent estimator αopt of α opt .

Due to space limitations, the proofs of the main nonstandard results are only sketched in the present manuscript but are available in full in an extended version.

Notations: Vectors (matrices) are denoted by lowercase (uppercase) boldface letters. {v a } n a=1 is the column vector v with (scalar or vector) entries v a and {V ab } n a,b=1 is the matrix V with (scalar or matrix) entries V ab . The operator

D(v) = D ({v a } n a=1
) is the diagonal matrix having (scalar or vector) v 1 , . . . , v n down its diagonal. The vector 1 n ∈ R n stands for the column vector filled with ones. The Dirac measure at x is δ x . The vector j a is the canonical vector of class C a defined by (j a ) i = δ i∈Ca and J

= [j 1 , . . . , j K ] ∈ {0, 1} n×K . The set C + is {z ∈ C, [z] > 0}.

II. MAIN RESULTS

A. Model and assumptions

We consider an n-node graph with

K classes C 1 , . . . , C K of sizes |C k | = n k .
To enforce a variation of the node degrees inside classes, we define a probability q i for node i to connect to any node in the graph. We define C ∈ R K×K a matrix of class weights C ab , independent of the q i 's, affecting the connection probability between nodes in C a and nodes in C b . The adjacency matrix A of the graph has then independent entries (up to symmetry) which are Bernoulli random variables with probability P ij = q i q j C ab ∈ (0, 1) when i ∈ C a and j ∈ C b and we set A ii = 0 for all i. In the dense regime under consideration, P ij = O(1). For convenience of exposition and without loss of generality, we assume that nodes indices are sorted by clusters and discard the nodes having no neighbor.

We shall perform spectral clustering on the family of matrices

L α ≡ n 2α-1 2 D -α A - dd T d T 1 n D -α (1)
where

D = D(d) with d = A1 n .
Clustering is asymptotically trivial when the weights C ab differ by O(1), as a vanishing error classification rate is easily guaranteed [START_REF] Gulikers | A spectral method for community detection in moderately-sparse degree-corrected stochastic block models[END_REF]. We shall instead consider the regime where the clustering performance is not asymptotically perfect. This regime is ensured by the following growth rate conditions. Assumption 1. As n → ∞, for all a,b ∈ {1, . . . , K} and i ∈ {1, . . . , n}:

1)

C ab = 1 + M ab √ n , where M ab = O(1); we shall denote M = {M ab } K a,b=1 .
2) q i are i.i.d. random variables with probability measure µ having compact support in (0, 1). We shall denote

m µ = tµ(dt). 3) ni n → c i > 0 and we will denote c = {c k } K k=1 . B.

Preliminary results

First note that, under Assumption 1, since

C ab → 1 as n → ∞, di √ d T 1n
is a uniformly consistent estimator of q i for 1 ≤ i ≤ n, i.e., max

1≤i≤n q i - d i d T 1 n a.s. -→ 0. (2)
Also observe that we can write

1 √ n A = 1 √ n qq T A d, √ n + 1 n q (a) q T (b) M ab K a,b=1 A d,1 + 1 √ n X Ar,1
where q (i) = [q n1+...+ni-1+1 , . . . , q n1+...+ni ] T ∈ R ni (n 0 = 0) and X = {X ij } n i,j=1 has independent (up to symmetry) entries of zero mean and variances σ 2 ij = q i q j (1 -q i q j ) + O(n -1 2 ). To study L α , which is not a standard random matrix model due to the dependency between D, dd T , and A, we shall perform a Taylor expansion of these matrices to retrieve an approximation Lα of L α which is asymptotically consistent in operator norm. To this end, note that A d, √ n , A d,1 , and A r,1 have spectral norms respectively of order O( √ n), O(1), and O(1), 1 so that 1

√ n A -dd T d T 1n
has all terms in O(1). Denoting D q = D(q), we next find

D = q T 1 n D q I n + √ n D -1 q q T 1 n (D {A d,1 1 n } + D {A a,1 1 n })
and the first order Taylor expansion of D -α is easily obtained as the last two terms in bracket are O(n -1 2 ). Putting all things together, we get the following approximation of L α .

Theorem 1. Let Assumption 1 hold and let L α be given by [START_REF] Fortunato | Community detection in graphs[END_REF]. Then, as n → ∞, L α -Lα → 0 in operator norm, almost surely, where

Lα = 1 m 2α µ 1 √ n D -α q XD -α q + UΛU T U = D 1-α q J √ n D -α q X1n q T 1n Λ = I K -1 K c T M I K -c1 T K -1 K -1 T K 0 . Since m 2α µ Lα is a perturbation of 1 √ n D -α q XD -α q
(full rank having zero mean and variance O(1) entries) by UΛU T (of maximum rank K), it follows a spiked random matrix model [START_REF] Chapon | The outliers among the singular values of large rectangular random matrices with additive fixed rank deformation[END_REF]. As illustrated in Figure 2, the spectrum of m 2α µ Lα (or equivalently of m 2α µ L α ) is asymptotically given by a compact spectrum with a density (red curve in Figure 2) and (sometimes) by additional isolated eigenvalues, hereafter called spikes.

Because of the (asymptotic) spiked model nature of L α , these spikes are only present beyond a phase transition threshold which depends on the norm of UΛU T (and thus indirectly of M) and on the spread of the compact component of the spectrum of L α (which essentially depends on α). In the following, we shall first study the conditions for this phase transition phenomenon and shall find a resulting "best choice" for α, which ensures the most favorable phase transition for worst case matrices M.

Beyond the phase transition, the dominant eigenvectors u i of L α shall correlate to the space spanned by the columns of U, so in particular to the D 1-α q j a , 1 ≤ a ≤ K, and not, as one would expect, to the canonical vectors j a (determining the classes) themselves. The dominant eigenvectors of L α are then biased by the degrees carried in D q and we thus propose to perform spectral clustering on the normalized vectors D α-1 u i rather than on the vectors u i themselves (recall that, up to a constant, D is a consistent estimator for D q ). These intuitive discussions are made rigorous in the subsequent section.

C. Eigenstructure of L α

We first aim at defining the aforementioned transition point beyond which eigenvalues isolate from the main spectrum of L α and thus non-trivial clustering ought to be achievable. To this end, we first identify the support S α of the limiting spectral distribution of m 2α µ L α . The latter is defined through its Stieltjes transform z → e α 2 (z) as follows. 2 Lemma 1 (Limiting spectral distribution). For z ∈ C + , the system

e α 1 (z) = q 1-2α -z -e α 1 (z)q 1-2α + e α 2 (z)q 2(1-α) µ(dq) (3) e α 2 (z) = q 2(1-α) -z -e α 1 (z)q 1-2α + e α 2 (z)q 2(1-α) µ(dq) (4)
2 Recall that the Stieltjes transform of a measure ν is defined, for z ∈ C \ supp(ν), as (t -z) -1 ν(dt). 2 δq 1 + 1 2 δq 2 , q 1 = 0.4, q 2 = 0.9, M defined by

M ii = 12, M ij = -4, i = j.
λ i ( M) > τ α ≡ -lim x↓S α 1 e α 2 (x)
with e α 2 defined in Lemma 1. In this case,

λ i (L α ) = - 1 m 2α µ e α 2 (λ i ( M))
.

The value τ α defined in Theorem 2 corresponds to a community detectability threshold beyond which sufficiently large eigenvalues λ i ( M) induce isolated eigenvalues in the spectrum of L α .

When separability is guaranteed, we subsequently show that the properly normalized dominant eigenvectors of L α tend to be close to linear combinations of the j a 's (and thus are essentially noisy step functions), the closeness to j a being ensured by the amplitude of the λ i ( M). This result, which again follows from classical spiked random matrix considerations (see e.g., [START_REF] Benaych | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF], [START_REF] Chapon | The outliers among the singular values of large rectangular random matrices with additive fixed rank deformation[END_REF]), is given precisely as follows.

Theorem 3. Under Assumption 1, let λ i ( M) and λ i (L α ) be an eigenvalue pair as defined in Theorem 2. We further assume λ i ( M) of unit multiplicity and denote u i the eigenvector associated to λ i (L α ). Then, letting ūi = D α-1 ui D α-1 ui and

Π = K a=1 jaj T a na , as n → ∞, almost surely, • If λ i ( M) ≤ τ α , ūT i Πū i → 0. • If λ i ( M) > τ α , lim inf n ūT i Πū i > 0.
• For all > 0, there exists T α > 0 such that, if λ i ( M) > T α , then lim inf n ūT i Πū i > 1 -. From Theorems 2-3, it is clear that the smaller τ α the more likely the condition λ i ( M) > τ α is met. Spectral clustering with L αopt where α opt = argmin α∈[0,1] {τ α } is thus "optimal" in the sense that it allows for non-trivial clustering when λ i ( M) is only slightly larger than τ αopt .

However, knowing α opt means being capable of estimating e α 2 (x) for all α ∈ [0, 1], as per Lemma 2. This is in fact doable thanks to Equation (2) which ensures that all q i , and thus µ, can be consistently estimated from the degrees d i . We thus have the following result.

.

While all these results provide strong hints on the expected performance of spectral clustering from L α , it still remains that the actual content of the dominant eigenvectors u i is unknown for generic settings. In the longer version of this article, we explicitly retrieve the "noisy plateaus" structure of the D α-1 u i and establish the consequences to clustering.

III. NUMERICAL RESULTS

This section illustrates the performances of our proposed method as compared to state-of-the-art spectral clustering methods, such as the BH approach (appropriate for sparse networks generated from the SBM). We provide simulations both for synthetic data generated from the DCSBM and for real world benchmarks commonly considered in the literature. The performance evaluation is the overlap to ground truth communities, defined in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] as

Overlap ≡ 1 n n i=1 δ(g i gi ) -1 K 1 -1 K
where g i and gi are the true and estimated label of node i, respectively. As the last step of all the spectral algorithms, we have performed 100 trials of the k-means algorithm on the K-1 leading eigenvectors (or D α-1 -normalized eigenvectors) u 1 , . . . , u K-1 starting from K points in R K-1 randomly extracted from rows of the matrix [u 1 , . . . , u K-1 ]. The selected k-means is chosen as the one with maximal modularity

1 d T 1n n i=1 n j=1 A ij - didj d T 1n δ(g i , gj ) [8].

A. Synthetic graphs

In this section, we consider graphs generated from the DCSBM. We first consider a rather realistic graph scenario (quite sparse with power law degree distribution). Figure 3 presents the theoretical (asymptotic) ratio between the largest eigenvalue λ of m 2α µ L α and the right edge S α of the main eigenvalue spectrum support, which one expects large to achieve good clustering performance, versus the amplitude of the elements of M. As predicted by our theoretical analysis, in the worst case scenario for M, this ratio is greater than 1 (and thus an eigenvalue escapes S α ) only when α = α opt . Besides, for fixed amplitude of the eigenvalues of M, the ratio obtained with α opt is seen to be the largest as compared to other values of α in the range [0, 1]; this suggests (without any theoretical support) that better clustering is achieved using L αopt even in non-worst case scenarios.

Figure 4 subsequently shows the overlap performance under the setting of Figure 3. Although we consider in that example a quite sparse regime which, for n = 2000, is not very consistent with our assumptions, our proposed method outperforms the BH approach, especially around the worst case values for M. It is worth mentioning that the empirically observed phase transitions closely match the theoretical ones (drawn in circles and the same as in Figure 3) but for the case α = 1. This mismatch is due to C taking values too far from 1, especially for large ∆, thereby no longer conforming Assumption 1.

We finally present in Figure 5 an example where the BH algorithm fails due to strongly heterogeneous node degrees. Assuming nodes connect with either low q 1 = 0.1 or high q 2 > q 1 intrinsic probability, we observe a sudden drop of the BH overlap once q 2 -q 1 is too large. This phenomenon is consistent with the fact, observed earlier in Figure 1, that BH 4 δq 1 + 1 4 δq 2 with q 1 = 0.1 and q 2 ∈ [0.1, 0.9], M defined by M ii = 10, M ij = -10, i = j, c i = 1 3 .

creates artificial communities out of nodes with the same q i parameter. This is a practical demonstration of the need for a proper eigenvector normalization to avoid degree biases.

B. Real world benchmarks

We finally confront the overlap performance on real world benchmarks in Table I. The best overlap score for each benchmark is set in boldface and quasi-equal scores in italic. Our approach largely outperforms the BH method on some benchmarks and has competitive performances on others. However note that, for so small network sizes, the performance achieved by L αopt may be quite unsatisfactory.

IV. CONCLUSION

We have introduced a new approach to community detection in large dense heterogeneous graphs using tools from random matrix theory. A salient feature of the approach lies in its automatically choosing the regularization of the modularity matrix which offers optimal performances. This work complements the recent breakthroughs in community detection, currently focusing on sparse but homogeneous graphs. A natural next step is to confront both approaches on a common ground, which is a challenging task due to the inadequacy (down to the mathematical tools) between sparse and dense network considerations. Studying the BH matrix in the dense heterogeneous regime is a natural first step into this investigation.

Fig. 2 .

 2 Fig. 2. Eigenvalues of m 2 µ L 1 , K = 3, n = 2000, c 1 = 0.3, c 2 = 0.3, c 3 = 0.4, µ = 1 2 δq 1 + 1 2 δq 2 , q 1 = 0.4, q 2 = 0.9, M defined by M ii = 12, M ij = -4, i = j.

  admits a unique solution (e α 1 (z), e α 2 (z)) ∈ (C + ) 2 . The map z → e α 2 (z) is the Stieljies transform of the limiting spectral distribution of m 2α µ L α with compact support S α = [-S α , S α ]. Our next objective is to determine the existence and, if so, position of the isolated eigenvalues of m 2α µ L α . Following now popular spiked model tools, this entails the following result.

Theorem 2 .

 2 Let Assumption 1 hold and let λ i ( M) be a non zero eigenvalue of M ≡ D(c) -cc T M. Then, for α ∈ [0, 1], there exists a corresponding isolated eigenvalue λ i (L α ) if and only if3

Lemma 2 .

 2 Define μ = 1 n n i=1 δ qi with qi = di √ d T 1n and êα i (z), i ∈ {1, 2}, Ŝα as in Lemma 1 but for µ replaced by μ. Then, as n → ∞, αopt → α opt almost surely, where αopt = argmin α∈[0,1] {τ α } with τα ≡ -1 lim x↓ Ŝα êα 2

Fig. 3 .

 3 Fig. 3. Largest eigenvalue λ of m 2α µ Lα as a function of the largest eigenvalue of M, M = ∆I 3 , c i = 1 3 , for ∆ ∈ [10, 150], µ a power law with exponent 3 and support [0.05, 0.3], for α ∈ {0, 1 4 , 1 2 , 3 4 , 1, αopt} (indicated below the graph). Here, αopt = 0.28. Circles indicate phase transition.

Fig. 4 .Fig. 5 .

 45 Fig. 4. Overlap for n = 3000, K = 3, c i = 1 3 , µ a power law with exponent 3 and support [0.05, 0.3], M = ∆I 3 , for ∆ ∈ [10, 150]. Here αopt = 0.28.

The notation O(•) is with respect to the operator norm.

The limit lim x↓S α e α 2 (x) is well defined in [-∞, 0] as x → e α 2 (x) is a growing negative function on the right side of S α .


Pierre-Antoine Ganaye

Michaël Sdika

Hugues Benoit-Cattin

TOWARDS INTEGRATING SPATIAL LOCALIZATION IN CONVOLUTIONAL NEURAL NETWORKS FOR BRAIN IMAGE SEGMENTATION

Keywords: brain MRI segmentation, CNN, spatial context, landmarks, probability atlas

Semantic segmentation is an established while rapidly evolving field in medical imaging. In this paper we focus on the segmentation of brain Magnetic Resonance Images (MRI) into cerebral structures using convolutional neural networks (CNN). CNNs achieve good performance by finding effective high dimensional image features describing the patch content only. In this work, we propose different ways to introduce spatial constraints into the network to further reduce prediction inconsistencies.

A patch based CNN architecture was trained, making use of multiple scales to gather contextual information. Spatial constraints were introduced within the CNN through a distance to landmarks feature or through the integration of a probability atlas. We demonstrate experimentally that using spatial information helps to reduce segmentation inconsistencies.

INTRODUCTION

Multi-Atlas methods [START_REF] Heckemann | Automatic anatomical brain mri segmentation combining label propagation and decision fusion[END_REF][START_REF] Klein | Mindboggle : Automated brain labeling with multiple atlases[END_REF] are one of the main approaches for the segmentation of brain structures : the segmentation maps produced individually by a set of atlas are combined for better quality. Machine learning can be used to improve the segmentation mapping or the label fusion of multi-atlas methods [START_REF] Sdika | Enhancing atlas based segmentation with multiclass linear classifiers[END_REF][START_REF] Wang | Multi-atlas segmentation with joint label fusion and corrective learning-an open source implementation[END_REF]. It can also be used by itself for MR brain image segmentation [START_REF] Marijn | Parallel multi-dimensional lstm, with application to fast biomedical volumetric image segmentation[END_REF][START_REF] Guha | Error corrective boosting for learning fully convolutional networks with limited data[END_REF][START_REF] Moeskops | Automatic Segmentation of MR Brain Images With a Convolutional Neural Network[END_REF].

Following the democratization of deep neural networks, [START_REF] Lee | Towards a deep learning approach to brain parcellation[END_REF] proposed an approach where a CNN is trained to predict the class of the center pixel of a patch. The segmentation map is obtained by applying the trained model as a sliding window over the image. Other variants [START_REF] Moeskops | Automatic Segmentation of MR Brain Images With a Convolutional Neural Network[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF] integrate patches of different sizes and resolutions, across parallel branches, in order to bring more contextual information to the network.

Recently another contribution based on an encoderdecoder architecture, shown its efficiency for brain segmentation [START_REF] Guha | Error corrective boosting for learning fully convolutional networks with limited data[END_REF] by firstly pre-training a CNN on a dataset annotated with FreeSurfer and optimizing a custom loss function during fine-tuning. This method produced state of the art segmentation maps. However, this type of architecture requires storing large features maps, which leads to memory consumptions issues. In this paper we chose to work on patch-based segmentation models, which provide interesting modeling properties to partially constrained problems.

Multi-atlas segmentation approaches [START_REF] Heckemann | Automatic anatomical brain mri segmentation combining label propagation and decision fusion[END_REF][START_REF] Klein | Mindboggle : Automated brain labeling with multiple atlases[END_REF][START_REF] Sdika | Enhancing atlas based segmentation with multiclass linear classifiers[END_REF][START_REF] Wang | Multi-atlas segmentation with joint label fusion and corrective learning-an open source implementation[END_REF] based on diffeomorphic registration methods are known to preserve the topology of the structures. Unlike atlas based segmentation methods, patch based CNN do not use the spatial position of the patch within the image volume. We propose and evaluate several ways to incorporate this prior knowledge to existing by-patch architectures, in order to yield better constraints learned from spatial features.

The integration of spatial information in a brain segmentation model was explored in [START_REF] Anbeek | Probabilistic segmentation of brain tissue in mr imaging[END_REF], where a nearest neighbor classifier combines voxel level intensities and spatial knowledge. Patch-based segmentation CNN using spatial features have also been proposed in the literature, [START_REF] Ghafoorian | Deep multi-scale locationaware 3d convolutional neural networks for automated detection of lacunes of presumed vascular origin[END_REF] used a combination of spatial coordinates and landmarks, which requires a pre-segmentation, [START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF] proposed to use a vector composed of distance to centroids, which requires to iterate a number of times to improve the estimates. In this paper, unlike [START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF] and [START_REF] Ghafoorian | Deep multi-scale locationaware 3d convolutional neural networks for automated detection of lacunes of presumed vascular origin[END_REF], we propose an architecture to integrate spatial context knowledge in any patch based CNN segmentation method, without requiring any kind of initialization nor iteration and is robust to degenerate cases such as nested structures.

We show the benefit of our approach, which does not require additional annotation, by applying it to a multiresolution CNN, on the multi-atlas brain segmentation challenge of MICCAI 2012.

In section 2. we first present the multi-resolution CNN and show how we introduce spatial context and information a priori. In section 3. we present the dataset and the evaluation metrics, followed by the implementation details. Finally in section 4. we illustrate and comment the obtained results.

METHODS

Our segmentation model is composed of a multi-resolution CNN (BaseNet), which is extended by tree additional branches (3dBranch, DistBranch, ProbBranch), each bringing more information and spatial context. Note that each of the additional branch can be used separately, as an individual component.

Multi-Resolution CNN

In Fig 1, we describe the main 2D multi-resolution classifier, named BaseNet in the remainder of the paper. It is composed of 3 branches inspired from [START_REF] Moeskops | Automatic Segmentation of MR Brain Images With a Convolutional Neural Network[END_REF][START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF]. The three branches are similar, having as input 25x25 patches but with three different scales. The two coarser patches are created by subsampling the original ones (of size 51 2 and 71 2 ), so that a larger context is accounted without increasing the number of parameters. The features produced by the three branches are concatenated and a 1x1 convolution is applied in order to reduce the feature space while keeping consistent information. Finally three consecutive fully connected layers are used to produce the scores together with the softmax function. To take advantage of the volumetric nature of the data, we try to test the importance of using a 3D branch (3dBranch). For doing this, a patch of size 15 3 is extracted and merged into the network with a dedicated branch. The size of the 3D patch is chosen to limit the brutal increase of parameters due to convolutions with 3D kernels while focusing on the fine scale information. The BaseNet architecture is used as an initial basis block throughout this article and will serve as reference to measure the performance of the proposed following branches.

Spatial information

Although the spatial position within the brain is a relevant information for brain structure segmentation, the classification produced by the BaseNet CNN relies on the patch content only. We propose to introduce a spatial representation of the patch position.

Let L be a set of landmarks defined uniformly along each axis of the input volume, where L ⊂ R 3 . For a patch centered at position x, we evaluate D ∈ R |L| , the euclidean distance vector of x to each landmarks in L. Unlike [START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF][START_REF] Ghafoorian | Deep multi-scale locationaware 3d convolutional neural networks for automated detection of lacunes of presumed vascular origin[END_REF], our landmarks are evenly distributed along each axis of the input image, they are not related to any region of the brain, thus do not require any pre-segmentation. As the landmarks are on a regular grid, D can be represented as a 3D image, which enables the use of convolutional layers. This distance image D serves as input to the DistBranch composed of convolutions with different kernel sizes, inspired from the inception module in [START_REF] Szegedy | Inception-v4, inception-resnet and the impact of residual connections on learning[END_REF] and finally merged with the second fully connected layer. The use of a radial basis function kernel to normalize the input distance image D has also been evaluated :

rbf (D) = exp(-αD 2 ) α ∈ R +

A priori via probability atlas

We also considered to account for the voxel position into the BaseNet segmentation by merging the CNN probability output with a more standard single probabilistic atlas segmentation. For a given voxel, the class probability vector according to the probabilistic atlas passes through three fully connected layers of size each, where is the number of classes. The output is summed to the BaseNet output just before the softmax function. This change is described as Prob-Branch in Fig 3. EXPERIMENTS

Data

The dataset is composed of 1.5T MRI from the OASIS project, it was distributed during the multi-atlas segmentation challenge of MICCAI 2012. The images were manually segmented into = 135 classes (structures and background). The original training dataset (15 images) was split into two non-overlapping sets : training (10 images), validation (5 images). The test dataset (20 images) is used to assert the performance of the models on unseen data. The images were affine registered to a reference atlas with FSL Flirt [START_REF] Jenkinson | Improved optimization for the robust and accurate linear registration and motion correction of brain images[END_REF]. All the images were skull striped, using a brain extraction software. The mean and standard deviation were estimated on the training set, all the images were finally mean centered and reduced.

Evaluation

The different methods discussed in this work were evaluated using the Dice coefficient, the Hausdorff distance d H and the mean surface distance d msd :

dice(X, Y ) = 2|X ∩ Y | |X| + |Y | , d H (X, Y ) = max{ sup x∈X d(x, Y ), sup y∈Y d(y, X)}, d msd (X, Y ) = 1 2   x∈X d(x, Y ) |X | + y∈Y d(y, X) |Y |   ,
where X and Y are two segmentation maps of a same label. d is defined as the minimal distance between a point and a set.

Implementation details

The number of landmarks was fine-tuned experimentally on the validation dataset, by varying the parameter from 3 3 to 10 3 points. We noticed an increase in precision until 7 3 and kept this parameter as a good balance between performance and processing time. We evaluated several representations for the distance image D : a 1D vector, a set of 2D images or a 3D volume. The 2D approach showed to be a good balance between the moderate performance of the 1D model and the cost of the 3D model because of the 3D convolutions. In the radial basis function, the value of α was set to 0.01.

The cross-entropy is used as the cost function. The numerical optimization was performed with SGD, with an initial learning rate lr 0 = 1e -3 and a momentum of 0.9. As in [START_REF] Chen | Deeplab : Semantic image segmentation with deep convolutional nets, atrous convolution, and fully connected crfs[END_REF] the learning rate was updated at each epoch with the poly rate policy :

poly(iter) = lr 0 * 1 - iter max iter power
Where iter is the index of the epoch, max iter the maximum number of epochs and power = 0.9. Batch of size 256 gave the best results. ReLU is the default activation function. Our 'Full' model is composed of all the branches, uses auxiliary losses and data augmentation.

Regularization & Auxiliary Loss

To prevent overfitting of the model, we used l 2 regularization, defined by an additional penalty term on the objective function λ 2 w 2 . Dropout [START_REF] Srivastava | Dropout : a simple way to prevent neural networks from overfitting[END_REF] was used for regularization, by randomly setting units output to 0. To ease the training of the network, we used auxiliary cost functions as advised in [START_REF] Szegedy | Inception-v4, inception-resnet and the impact of residual connections on learning[END_REF], at two different locations (BaseNet, DistBranch). An auxiliary loss consists of a fully connected layer attached to a branch, followed by a softmax function. The global loss is composed by the network's loss and the auxiliary losses.

Class imbalance & Data Augmentation

Because the anatomical regions of the brain have varying volumes, sampling from the original distribution produces class imbalance. We tried to balance the classes by adjusting the weights accordingly in the cross-entropy cost function, unfortunately we did not see any improvement. In order to increase the variability of the images for better generalization, random data augmentations were applied on the largest 2D patch, by combining rescale with a factor in the range [0.9 ;1.1] and rotation in the range [-10 ;10] degrees.

RESULTS

To evaluate the optimal way of integrating the distances to landmarks into the BaseNet CNN, we tested the performance of four architectures, table 2 shows the results obtained. It is clear that using the radial basis function to normalize does not really help, consequently we chose not to include it. As a matter of fact, it seems the CNN is able to adjust the signal values with the help of linear operations, such as convolutions. This assumption is underpinned by the poor result of the model 'distance with rbf'. We decide ultimately to integrate 2D distances into the network with convolution layers and without any kind of normalization. Figure 2 shows an example of segmentation maps we produced with the tested models. A real performance gap can be noticed between BaseNet(e) and BaseNet+DistBranch (d), where the first detects background between the left and right lateral ventricles and the second is able to recover smooth structures.

In table 1, we can notice the impact of each branch in this incremental setup. Adding each of them successively brings better results. The 2D multi-resolution model (BaseNet) combined with the distance integration (DistBranch) shows a noticeable decrease in the average and standard deviation of the Hausdorff distance, thus reducing serious segmentation issues, with the help of better spatial constraints. The best model is finally a combination of all the proposed branches, leading to an average dice of 0.748. In comparison, [START_REF] De | Deep neural networks for anatomical brain segmentation[END_REF] which is the only to our knowledge to have used a by-patch segmentation approach for the original 135 classes problem, proposed a model composed of 30M parameters and reached an average dice of 0.725. Our model has 10 order of magnitude less parameters, with a better average dice. With this model, we would have been ranked 5th of the multi-atlas segmentation challenge at MICCAI 2012, with a segmentation time per image of approximately 9 minutes.

We briefly compare to a UNet [START_REF] Ronneberger | Convolutional Networks for Biomedical Image Segmentation[END_REF] like encoder-decoder architecture inspired from [START_REF] Guha | Error corrective boosting for learning fully convolutional networks with limited data[END_REF], with skip-connections and max unpooling. It was trained to segment slice by slice, optimized only with cross-entropy and dice loss, on the same dataset. It showed encouraging dice similarity, but poor Hausdorff performance, demonstrating that patch based segmentation is still a competitive task for brain segmentation.

CONCLUSION

In this study, we set a framework to integrate spatial constraints into any patch based classification network. We showed that integrating distances to landmarks into a 2D multi-resolution CNN can help reduce segmentation incoherence. Adding information from a probability atlas together with 3D patch helps to minimize segmentation errors further. Encoder-decoder segmentation architecture showed promising results in terms of speed and performance, future work could be done to introduce specific spatial constraints in such model to reduce segmentation inconsistencies.

ACKNOWLEDGEMENT

This work was funded by the CNRS PEPS "ReaDAP-TIM" and was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program "Investissements d'Avenir"(ANR-11-IDEX-0007) operated by the French National Research Agency (ANR). We would like to thank the E.E.A doctoral school for providing a financial support. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan X Pascal GPU used for this research. Also we would like to thank the IN2P3 computing center for sharing their resources.

Fig. 1 :

 1 Fig. 1: Architecture of the network, composed of a 2D multi-resolution CNN (BaseNet), 3dBranch for the 3D input, DistBranch for the distances to landmarks and ProbBranch to integrate knowledge from the probability atlas.

Fig. 2 :

 2 Fig. 2: Illustration of the segmentations. Coronal slice (a) and associated segmentation maps : ground truth (b), Full (c), Base-Net+DistBranch (d) and BaseNet (e). The segmentations were obtained on one patient of the test dataset.

Table 1 :

 1 [START_REF] Heckemann | Automatic anatomical brain mri segmentation combining label propagation and decision fusion[END_REF]. Distance and similarity metrics for each best performing models. MSD is the mean surface distance and N param is the number of parameters of the network. All the metrics are averaged over the test dataset. (average ± standard deviation)

	Model

Table 2 :

 2 Comparison of methods to integrate distance to landmarks, performance measured using the dice similarity and the Hausdorff distance. (average ± standard deviation)