Marius Tucsnak
email: marius.tucsnak@u-bordeaux.fr

Julie Valein
email: julie.valein@univ-lorraine.fr

Chi-Ting Wu
email: chi-ting.wu@univ-lorraine.fr

Finite dimensional approximations for a class of infinite dimensional time optimal control problems

Keywords: distributed parameter systems, optimal control, numerical approximation AMS subject classifications. 93C25, 93B07, 93C20, 35R30

In this work we study the numerical approximation of the solutions of a class of abstract parabolic time optimal control problems with unbounded control operator. Our main results assert that, provided that the target is a closed ball centered at the origin and of positive radius, the optimal time and the optimal controls of the approximate time optimal problems converge (in appropriate norms) to the optimal time and to the optimal controls of the original problem. In order to prove our main theorem, we provide a nonsmooth data error estimate for abstract parabolic systems.

Introduction

Time optimal control of infinite dimensional systems is a subject of growing interest, motivated by numerous applications in domains such as guidance of complex systems or temperature regulation in large buildings. In recent year, using new tools from infinite dimensional systems theory, the literature devoted to this topic grew in a considerable manner (see Arada and Raymond [START_REF] Arada | Time optimal problems with Dirichlet boundary controls[END_REF], Barbu [START_REF] Barbu | Analysis and control of nonlinear infinite-dimensional systems[END_REF], [START_REF] Barbu | Nonlinear differential equations of monotone types in Banach spaces[END_REF], Fattorini [START_REF] Fattorini | Time-optimal control of solutions of operational differential equations[END_REF], [START_REF] Fattorini | Infinite Dimensional Linear Control Systems. The Time Optimal and Norm Optimal Control Problems[END_REF], [START_REF] Fattorini | Infinite Dimensional Optimization and Control Theory[END_REF], Kunisch and Wang [START_REF] Kunisch | Time optimal control of the heat equation with pointwise control constraints[END_REF], Kunisch and Wachsmuth [START_REF] Kunisch | On time optimal control of the wave equation, its regularization and optimality system[END_REF], Li and Yong [START_REF] Li | Optimal control theory for infinite-dimensional systems[END_REF] and Tröltzsch [START_REF] Tröltzsch | On generalized bang-bang principles for two time-optimal heating problems with constraints on the control and the state[END_REF]). The specific case of time optimal control for systems governed by parabolic PDE's has numerous applications, from which we quote optimization of building thermal storage (see, for instance [START_REF] Henze | Evaluation of optimal control for active and passive building thermal storage[END_REF] and references therein).

The aim of this paper, containing results partially announced in [START_REF] Tucsnak | Numerical approximation of some time optimal control problems[END_REF], is to study the approximation of the solutions of time optimal (internally or boundary) control problems for a class of infinite dimensional linear systems by projecting the original problem on an appropriate family of finite dimensional spaces. This is a delicate question since, as shown in the above mentioned references, time optimal controls are usually highly oscillating functions (due to the bang-bang property). As far as we know, the only papers having already investigated this issue are Knowles [START_REF] Knowles | Finite element approximation of parabolic time optimal control problems[END_REF], Wang and Wang [START_REF] Wang | Finite element approximations of optimal controls for the heat equation with end-point state constraints[END_REF], Wang and Zheng [START_REF] Wang | An approach to the optimal time for a time optimal control problem of an internally controlled heat equation[END_REF] which investigated finite elements approximation for systems governed by the heat 1 equation with internal controls. To our knowledge, the similar results obtained in these papers are unknown to boundary time optimal control problems. A similar topic dealing with the homogenization of parabolic system, where we study the asymptotic behavior of solutions and (time-)optimal controls has been developed in recent years (see for example, Carja [START_REF] Cârjȃ | On variational perturbations of control problems: minimum-time problem and minimum-effort problem[END_REF], [START_REF] Cârjȃ | The time optimal problem for boundary-distributed control systems[END_REF], Castro and Zuazua [START_REF] Castro | Some topics on the control and homogenization of parabolic partial differential equations[END_REF] and Tebou [START_REF] Tebou | Uniform null controllability of the heat equation with rapidly oscillating periodic density[END_REF] ).

To state our results we need some notation. Let X and U be real Hilbert spaces, and let A 0 : DpA 0 q Ñ X be a strictly positive operator. It is known that ´A0 generates an exponentially stable analytic semigroup, denoted by T t . We denote by }.} (resp. x., .y) the norm (resp. the inner product) on X. For γ ą 0 we denote by X γ the space DpA γ 0 q, endowed with the graph norm. For γ ă 0, X γ stands for the dual of X ´γ with respect to the pivot space X. We also introduce an operator B P LpU, X ´αq with 0 ď α ă 1 2 , called control operator. In this paper we consider systems of the form 9 zptq `A0 zptq " Buptq pt ě 0q, (

zp0q " z 0 pz 0 P Xq,

where u P L 8 pr0, `8q; U q is the input function and z is the state trajectory. It is well known (this follows, for instance, by combining Theorem 4.4.3 and Proposition 5.1.3 in [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]) that if z 0 P X and u P L 8 pr0, `8q; U q, there exists a unique solution z P Cpr0, `8q; Xq of (1.1)-(1.2) and z satisfies zptq " T t z 0 `Φt u, where Φ t u "

ż t 0 T t´σ Bupσqdσ. (1.3)
Given ε ą 0, denote by Bp0, εq the closed ball centered in zero and of radius ε in X. We consider the time optimal control problem which states as follows pTPq Determine τ 0 ą 0 such that τ 0 " mintτ ě 0 | there exists u P L 8 pr0, `8q; U q s.t. }u} L 8 pr0,`8q;U q ď 1 and zpτ q P Bp0, εqu, and the corresponding optimal controls u 0 .

Denote

U ad " tu P L 8 pr0, `8q; U q | }u} L 8 pr0,`8q;U q ď 1u.

(1.4)

It is well-known that the above optimal time τ 0 and optimal control u 0 exist under additional assumptions (see, for instance, [START_REF] Micu | Time optimal boundary controls for the heat equation[END_REF]).

Let pV h q hą0 (resp. pU h q hą0 ) be a family of finite dimensional subspaces of X 1 2 (resp. U ), which are normed spaces when endowed with the restriction of the norm of X 1 2 (resp. U ). We denote P h (resp. Q h ) the orthogonal projector from X onto V h (resp. U onto U h ). For each h ą 0 we consider the following system:

9 z h ptq `Ah z h ptq " B h u h ptq pt ě 0q, (1.5) z h p0q " P h z 0 , (1.6)
where u h P L 8 pr0, `8q; U h q, pA h q hą0 is defined by

ă A h ϕ, ψ ą " ă A 1 2 0 ϕ, A 1 2 0 ψ ą pϕ, ψ P V h q, (1.7)
and B h P LpU, V h q is defined by: xB h u, ϕy " xu, B ˚ϕy U pϕ P V h , u P U q.

(1.8)

The above system is the Galerkin approximation of (1.1)-(1.2), and its solution z h can be written as

z h ptq " T t,h P h z 0 `Φt,h u h , where Φ t,h u h " ż t 0 T t´σ,h B h u h pσqdσ, (1.9)
where T t,h " expp´tA h q is the semigroup generated by A h . Denote by Bh p0, εq the closed ball centered in zero in V h with radius ε. For each h ą 0, we consider the time optimal control problem for the above system (1.5)-(1.6) which states as follows: pTP h q Determine τ h ą 0 such that τ h " mintτ ě 0 | there exists u h P L 8 pr0, `8q; U h q s.t. }u h } L 8 pr0,`8q;U h q ď 1 and z h pτ q P Bh p0, εqu, and the corresponding optimal controls u h in the admissible set U ad,h , where U ad,h is defined by U ad,h " tu P L 8 pr0, `8q; U h q | }u} L 8 pr0,`8q;U h q ď 1u.

(1.10)

The goal of this work is to study the convergence of τ h to τ 0 and of u h to u 0 when h Ñ 0. To this aim, we need appropriate assumptions on the approximation properties of the spaces pV h q hą0 and pU h q hą0 . More precisely, we assume that there exist θ ą 0, h 1 ą 0, C ą 0 such that for every h P p0, h 1 q and 0 ď γ ď 1 we have:

pC1q }x ´Ph x} ď Ch θγ }x} γ for every x P X γ , pC2q
A h is uniformly (with respect to h) analytic, in the sense of the Definition 3.2 below, pC3q }pI ´Ph qB} LpU,Xq ď Ch θp1´αq , pC4q }P h B} LpU,V h q ď Ch ´θα , pC5q lim hÑ0 }Q h u ´u} U " 0 for every u P U , where α represents the unboundedness degree of the operator B (introduced in the beginning of this section).

We are now in a position to state the main results of this paper:

Theorem 1.1.
With the notation and assumptions on operators A 0 , B, A h , B h , assume that pC1q-pC5q hold and that z 0 P X, }z 0 } ą ε. Then lim hÑ0 τ h " τ 0 .

Theorem 1.2. With the notation and assumptions on operators A 0 , B, A h , B h , assume that pC1q-pC5q hold and that the only η P X for which there exists an open non-empty interval I with B ˚Tt η " 0 for t P I is η " 0. Then we have

u h Ñ u 0 strongly in L 2 pr0, T s; U q,
where T " 2 lnp}z 0 }{εq λ 1

(λ 1 is the first eigenvalue of the operator A 0 ) and extending pu hq h and u 0 to time T by zero.

It is worth mentioning that similar convergence results have been obtained in [START_REF] Knowles | Finite element approximation of parabolic time optimal control problems[END_REF] and [START_REF] Wang | An approach to the optimal time for a time optimal control problem of an internally controlled heat equation[END_REF]. In [START_REF] Knowles | Finite element approximation of parabolic time optimal control problems[END_REF], the control space takes only in finite dimensional subspaces of U and in [START_REF] Wang | An approach to the optimal time for a time optimal control problem of an internally controlled heat equation[END_REF] the initial data z 0 is taken in X 1 can be seen as generalizations of those obtained in [START_REF] Wang | An approach to the optimal time for a time optimal control problem of an internally controlled heat equation[END_REF], in the sense that we consider a class of abstract problems including the PDE systems studied there. The novelty of our results is that we weaken the regularity assumptions on the initial data with z 0 P X and also that we weaken the assumption on control operator B with B P LpU, X ´αq, 0 ď α ă 1 2 . The outline of the remaining part of this paper is as follows. Section 2 contains some necessary background on time optimal control problems. In Section 3, we provide some error estimate results with smooth initial data which play an essential role to prove our Theorem 1.1. In Section 4, we provide the proof of our main theorems. In Section 5, we apply our abstract results to some equations. Throughout the paper, we denote by C a positive constant that may change from line to line.

2 Some background on time optimal control problem This section, in which we continue to use the notation and assumptions introduced in Section 1 for X, U , A 0 , B, is devoted to some background on the time optimal control problem pTPq. We first introduce some basic notions on controllability (see for example Ch6. and Ch11. in [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]).

For τ ą 0, we denote Ψ τ P LpX, L 2 pr0, τ s; U qq the observation operator of pA 0 , B ˚q on r0, τ s defined by pΨ τ z 0 qpσq " B ˚Tσ z 0 pz 0 P X, σ P r0, τ sq.

It is clear that (see, for instance, [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Ch. 4]),

Φ τ " R τ Ψ τ pτ ą 0q, (2.11)
where R τ is the reflection operator on L 2 pr0, τ s; U q defined by R τ uptq " upτ ´tq.

We then recall the following existence result, which can be proved by using standard techniques providing the existence of time optimal controls for linear systems. We refer to [START_REF] Fattorini | Infinite Dimensional Linear Control Systems. The Time Optimal and Norm Optimal Control Problems[END_REF] or Lions [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF] for a detailed description of the methodology or to Micu et al [START_REF] Micu | Time optimal boundary controls for the heat equation[END_REF]Proposition 2.6] for a slightly more general version of this result.

Proposition 2.1. For every z 0 P X and ε ą 0, the time optimal control problem pTPq admits at least one solution pτ ˚, u ˚q.

We need below a special version of the maximum principle. Although maximum principle is a classical tool in optimal control problems, we were unable to find the version we need (abstract setting with unbounded control operators) in the existing literature. Therefore we give the precise statement and a short proof below. Theorem 2.2. With the notation of Proposition 2.1, assume that τ ˚ą 0. Then the time optimal control u ˚satisfies the maximum principle, i.e., there exists ξ P X, ξ " 0 such that,

xu ˚pσq, pΨ τ ˚ξq pτ ˚´σqy U " max vPU, }v}ď1
xv, pΨ τ ˚ξq pτ ˚´σqy U , pσ P p0, τ ˚q a.e.q.

(2.12)

Moreover, ξ satisfies the transversality condition, i.e.:

xξ, z ´Tτ ˚z0 ´Φτ ˚u˚y ě 0 pz P Bp0, εqq.

(2.13)

Proof. Let pτ ˚, u ˚q be a solution of the time optimal control in Proposition 2.1 and let

z 1 " T τ ˚z0 `Φτ ˚u˚.
Moreover, for each z 0 P X and τ ą 0 denote

U τ :" tu| r0,τ s | u P U ad u and R τ z 0 " T τ z 0 `Φτ U τ .
We claim that z 1 P BR τ ˚z0 X BBp0, εq.

(2.14)

Indeed, it can be easily checked that we necessarily have }z 1 } " ε. If we admit, by contradiction, that z 1 is an interior point of R τ ˚z0 , this implies that there exists z 2 P R τ ˚z0 with }z 2 } ă ε. Consequently, there exists

u 2 P U τ ˚such that T τ ˚z0 `Φτ ˚u2 " z 2 .
Since the map t Þ Ñ }T t z 0 `Φt u 2 } is continuous and that τ ˚ą 0, we obtain that there exists τ P p0, τ

˚q such that }T τ z 0 `Φτ u 2 } " ε.
This contradicts the fact that pτ ˚, u ˚q is a solution of our time optimal control problem.We have thus proved our claim (2.14).

Thus we clearly have that the sets R τ ˚z0 and Bp0, εq are non-empty, convex and they have no common point. Moreover, since Bp0, εq is open, we can apply the geometric version of the Hahn Banach theorem (see, for instance [START_REF] Brezis | Analyse fonctionnelle[END_REF]Theorem 1.6]) to obtain that there exists a hyperplane separating R τ ˚z0 and Bp0, εq. This means that there exists α P R and ξ P X, ξ " 0 such that

xξ, ηy ď α pη P R τ ˚z0 q, (2.15)
xξ, ηy ě α pη P Bp0, εqq.

(

The two above inequalities and (2.14) imply that xξ, z 1 y " α.

(2.17)

From (2.15) and (2.17) we deduce that xξ, z 1 y ě xξ, ηy pη P R τ ˚z0 q, which combined with the duality (2.11), clearly implies the maximum principle (2.12). Moreover, we easily deduce the transversality condition (2.13) from (2.16) and (2.17).

Remark 2.3. To be more precise, we can deduce from the transversality condition (2.13) that:

ξ " ´k pT τ ˚z0 `Φτ ˚u˚q , for some positive constant k.

(2.18)

Corollary 2.4. With the notation of Theorem 2.2, assume that τ ˚ą 0 and that the assumption of Theorem 1.2 holds, i.e., the only η P X for which there exists an open non-empty interval I with B ˚Tt η " 0 for t P I is η " 0. Then the time optimal control u ˚is unique and it has the bang-bang property: }u ˚ptq} " 1 pt P p0, τ ˚q a.e.q.

Proof. We first remark that the following statement holds: If η P X such that there exists a subset e Ă r0, τ s of positive measure with B ˚Tt η " 0 for t P e, then η " 0. Indeed, the facts that the map t Þ Ñ B ˚Tt η is analytic (which is due to the analyticity of T t ) and that it vanishes in e imply that this map vanishes for t in some non-empty open interval I. Therefore, according to the hypothesis, we have η " 0. With the above property, it is clear that pΨ τ ˚ξq pτ ˚´tq ‰ 0 for almost every t P p0, τ ˚q. Indeed, if it is not the case, we deduce from the above property that ξ " 0, which contradicts Theorem 2.2. We can deduce from the maximum principle (2.12) that u ˚ptq " pΨ τ ˚ξqpτ ˚´tq }pΨ τ ˚ξqpτ ˚´tq} . Thus, the bang-bang property of the time-optimal control holds.

The uniqueness of the time-optimal control can be deduced from the strict convexity of U . This ends the proof.

Nonsmooth data error estimates for abstract parabolic equations

In this section, we gather, for easy reference, some results providing error estimates for the semidiscrete Galerkin approximation of some abstract parabolic systems. Our results are strongly inspired by Lasiecka and Triggiani [20,Ch.4] and Badra [3,Ch.5], where they deal with a more general framework. We continue to use the same notation and assumption introduced in Section 1 for X, U , A 0 , B, A h , B h and V h and we recall that T t (resp. T t,h ) is the semigroup generated by ´A0 (resp. ´Ah ). We then consider the following linear equation and its approximated scheme (where we do not approximate u by u h ):

9 zptq `A0 zptq " Buptq pt ě 0q, (3.1)
zp0q " z 0 P X, (3.2)

9 z h ptq `Ah z h ptq " B h uptq pt ě 0q, (3.3)
z h p0q " P h z 0 , (3.4)
with u P L 8 pr0, `8q; U q. We first recall a classical characterization of analytic semigroups: Proposition 3.1. T t is an analytic semigroup if and only if there exist C, δ ą 0 such that

Σ Ă ρpA 0 q and } pλI ´A0 q ´1 } LpXq ď C |λ| pλ P Σq, (3.5)
where Σ is defined by Σ :

" λ P C | |argpλq| ď π 2 `δ( .
It is well known (see, for instance, [24, Ch. 2.6]) that if T t is analytic then

}A γ 0 T t } LpXq ď C t γ pt ě 0, γ P p0, 1qq. (3.6)
We also need the concept of uniform analyticity of the family of semigroups pT t,h q hě0 . Definition 3.2. T t,h is an uniformly analytic semigroup if there exist C independant of h, δ ą 0 such that for every h ą 0 we have

Σ Ă ρpA h q and }pλI ´Ah q ´1} LpXq ď C |λ| pλ P Σq,
where Σ is defined in Proposition 3.1.

Similarly, the above definition implies the following property:

}A γ h T t,h } LpXq ď C t γ pt ą 0, γ P r0, 1sq. (3.7)
It is also clear that B h u " pP h Bq u pu P U q, and (see assumption pC4q) }B h } LpU,V h q ď Ch ´θα .

We recall in the following some well-known results:

Lemma 3.3.
Assume that T t is analytic. Then, for any 0 ď γ ď 1, there exists C ą 0 such that the following estimates hold:

1. }A γ 0 pλI ´A0 q ´1} LpXq ď C |λ| 1´γ pλ P Σq. 2. T t z " 1 2πi
ş Γ e λt pλI ´A0 q ´1 z dλ pz P Xq, where Γ is the path composed from the two rays ρe iψ and ρe ´iψ , 0 ă ρ ă 8 and π{2 ă ψ ă π{2 `δ, Γ being oriented so that λ increases along Γ ( λ denotes the imaginary part of λ).

Proof. The proof of the first estimate can be found in [24, Ch. 2.6] and the proof of the last assertion can be found in [24, Th. 7.7, Ch. 1].

Similar to Lemma 3.3, it is clear that we have the following lemma: Lemma 3.4. Assume that T t,h is uniformly analytic. Then, for any 0 ď γ ď 1, there exists C ą 0 such that the following statement hold: 1. }A γ h pλI ´Ah q ´1} LpXq ď C |λ| 1´γ pλ P Σq. 2. T t,h z " 1 2πi ş Γ e λt pλI ´Ah q ´1 z dλ pz P V h q.

We gather in the proposition below, with no claim of originality, several error estimates which play a central role in the proof of our main result. Some of the proofs are very similar to those in [START_REF] Badra | Stabilisation par feedback et approximation des équations de Navier-Stokes[END_REF]Ch.5], so that they are detailed. Proposition 3.5. Assume that the assumption pC1q ´pC4q hold. Then, there exists C ą 0, such that for every h ą 0, we have the following properties:

1. }A ´1 0 ´A´1 h P h } LpXq ď Ch θ . 2. }A ´1 0 pB ´Bh q} LpU,Xq ď Ch θp1´αq .
3. }pλI ´A0 q ´1 ´pλI ´Ah q ´1P h } LpXq ď Ch θ pλ P Σq. 4. }pλI ´A0 q ´1B ´pλI ´Ah q ´1B h } LpU,Xq ď Ch θp1´αq pλ P Σq. 5. }T t,h P h ´Tt } LpXq ď Ct ´1h θ . 6. }T t,h B h ´Tt B} LpU,Xq ď Ct ´1h θp1´αq . 7. }T .,h B h ´T. B} L 1 pr0,ts;LpU,Xqq ď Ch θp1´αq p| ln h| `t `1q .

Proof.

1. We skip this proof which is based on Proposition 5.1.1 and Proposition 5.1.3 in [START_REF] Lasiecka | Control Theory for Partial Differential Equations: Continuous and Approximation Theories[END_REF].

2. We skip this proof as it is similar to the proof of previous assertion by using assumption pC3q.

3. Notice at first that by simple calculation, it is clear that we have pλI ´A0 q ´1 " ´A´1 0 `λpλI ´A0 q ´1A ´1 0 and pλI ´Ah q ´1P h " ´A´1 h P h `λpλI ´Ah q ´1A ´1 h P h .

Denote M h pλq " pλI ´A0 q ´1 ´pλI ´Ah q ´1P h , we then have:

M h pλq " ´A´1 0 `A´1 h P h `λpλI ´A0 q ´1A ´1 0 ´λpλI ´Ah q ´1A ´1 h P h " ´`A ´1 0 ´A´1 h P h ˘`λpλI ´A0 q ´1 `A´1 0 ´A´1 h P h ˘`λpλI ´A0 q ´1A ´1 h P h ´λpλI ´Ah q ´1A ´1 h P h " ´`I ´λpλI ´A0 q ´1˘`A ´1 0 ´A´1 h P h ˘`λM h pλqA ´1 h P h .
This leads to M h pλq " ´`I ´λpλI ´A0 q ´1˘p A ´1 0 ´A´1 h P h qpI ´λA ´1 h P h q ´1. Notice that pI ´λA ´1 h P h q ´1 " I ´Ph ´Ah pλI ´Ah q ´1P h .

Indeed, we have pI ´λA ´1 h P h q `I ´Ph ´Ah pλI ´Ah q ´1P h " I ´Ph ´Ah pλI ´Ah q ´1P h ´λA ´1 h P h `λA ´1 h P h `λA ´1 h A h pλI ´Ah q ´1P h " I ´Ph `pλI ´Ah qpλI ´Ah q ´1P h " I.

In the same manner, `I ´Ph ´Ah pλI ´Ah q ´1P h ˘pI ´λA ´1 h P h q " I. Thus, (3.8) holds.

By analyticity assumption (3.5) and assertion 1. in Lemma 3.4, we know that there exists a constant C independent of λ and h such that, }λpλI ´A0 q ´1} LpXq ď C |λ| |λ| " C and }A h pλI ´Ah q ´1} ď C pλ P Σq, which leads to, using assertion 1. of Proposition 3.5

}M h pλq} LpXq ď }I ´λpλI ´A0 q ´1} LpXq Ch θ }I ´Ph ´Ah pλI ´Ah q ´1P h } LpXq ď Ch θ .
This ends the proof of assertion 3..

4.

Denote M h,B pλq " pλI ´A0 q ´1B ´pλI ´Ah q ´1B h . Notice that pλI ´A0 q ´1 " A 0 pλI ´A0 q ´1A ´1 0 .

Then, we have: M h,B pλq " pλI ´A0 q ´1pB ´Bh q ``pλI ´A0 q ´1 ´pλI ´Ah q ´1˘B h " A 0 pλI ´A0 q ´1A ´1 0 pB ´Bh q ``pλI ´A0 q ´1 ´pλI ´Ah q ´1˘B h .

Since A 0 pλI ´A0 q ´1 is a bounded operator (by assertion 1. in Lemma 3.3), it is clear that by using assertions 2. and 3. of Proposition 3.5 and assumption pC4q, we have:

}M h,B pλq} LpU,Xq ď C}A ´1 0 pB ´Bh q} LpU,Xq `}pλI ´A0 q ´1 ´pλI ´Ah q ´1} LpXq }B h } LpU,V h q ď Ch θp1´αq `Ch θ }B h } LpU,V h q ď Ch θp1´αq .

5.

It is clear that, by using assertion 2. in Lemma 3.3 and assertion 2. in Lemma 3.4, we have:

T t ´Tt,h P h " 1 2πi
ż Γ e λt `pλI ´A0 q ´1 ´pλI ´Ah q ´1P h ˘dλ.

Then, by using assertion 3. in Proposition 3.5, we have:

}T t ´Tt,h P h } LpXq " C ˆżΓ |e λt |dλ ˙hθ ď Ct ´1h θ .
This ends the proof.

6. We skip this proof as it is similar to the proof of the previous assertion by using assertion 4. in Proposition 3.5.

We calculate:

}T .,h B h ´T. B} L 1 pr0,ts;LpU,Xqq "

ż h θ 0 }T σ,h B h ´Tσ B} LpU,Xq dσ`ż t h θ }T σ,h B h ´Tσ B} LpU,Xq dσ " I 1 `I2 .
By using assertion 6. of Proposition 3.5, we have:

I 2 ď Ch θp1´αq ˆż t h θ σ ´1dσ ˙ď Ch θp1´αq ˆż 1 h θ σ ´1dσ `ż t 1 1dσ ˙ď Ch θp1´αq p| ln h| `t `1q .
Note that the assumption B P LpU, X ´αq means that A ´α 0 B is a bounded operator in LpU, Xq.

Then, in order to deal with I 1 , using also the analyticity of T (see (3.5)), for h small enough, we have, for every 0 ď α ă 1 2 ,

}T σ B} LpU,Xq ď }T σ A α 0 } LpXq }A ´α 0 B} LpU,Xq ď Cσ ´α.
Similarly, by using again the assumption B P LpU, X ´αq, it is clear that A ´α h B h is a bounded operator in LpU, Xq. Thus, we have

}T σ,h B h } LpU,Xq ď }T σ,h A α h } LpU,Xq }A ´α h B h } LpU,Xq ď Cσ ´α, using (3.7).
Finally, by the two previous inequalities, we have :

I 1 ď C ż h θ 0 σ ´αdσ ď Ch θp1´αq .
Thus, we have }T .,h B h ´T. B} L 1 pr0,ts;LpU,Xqq ď Ch θp1´αq p| ln h| `t `1q .

This ends the proof.

We will now present the main result of this section which gives us a nonsmooth initial data error estimate.

Theorem 3.6. Let z 0 P X, let z be the solution of (3.1),(3.2) and let z h be the solution of (3.3), (3.4). Moreover, assume that pC1q ´pC4q hold. Then, there exist C ą 0, h ą 0 such that for h P p0, hq we have the following error estimate: }zptq ´zh ptq} X ď Ch θ t ´1}z 0 } X `Ch θp1´αq p| ln h| `t `1q }u} L 8 pr0,`8q;U q , @t ą 0.

(3.9)

Proof. Denote K h ptq " T t B ´Tt,h B h . We recall that z, solution of system (1.1)-(1.2), satisfies (1.3) and z h , solution of (1.5)-(1.6), satisfies (1.9). It is clear that

}zptq ´zh ptq} X ď }T t z 0 ´Tt,h P h z 0 } X `› › › › ż t 0 K h pt ´σqupσqdσ › › › › X .
By using assertions 5. and 7. in Proposition 3.5, we have:

}zptq ´zh ptq} X ď Ch θ t ´1}z 0 } X `}K h ˚u} L 1 pr0,ts;Xq ď Ch θ t ´1}z 0 } X `}K h } L 1 pr0
,ts;LpU,Xqq }u} L 8 pr0,ts;U q ď Ch θ t ´1}z 0 } X `Ch θp1´αq p| ln h| `t `1q }u} L 8 pr0,ts;U q .

This ends the proof.

Proof of the main results

We denote τ 0 pz 0 q (resp.τ h pz 0,h q) the optimal time associated with the initial state z 0 P X (resp. z 0,h P V h ) and with final state in Bp0, q (resp. Ď B h p0, εq). We use also the notation zpt, z 0 , uq for the solution of the system (1.1)-(1.2) (resp. z h pt, z 0,h , u h q for the system (1.5)-(1.6)) associated with the initial state z 0 (resp. z 0,h ) and with the control u (resp. u h ) at time t.

We first recall a generalized Aubin-Lions theorem (Theorem 4.1) and a standard energy estimate (Lemma 4.2).

Theorem 4.1. Let T ą 0 and let Y 0 , Y 1 , Y 2 be Banach spaces such that Y 0 Ă Y 1 Ă Y 2 . Assume that Y i is reflexive for i " 0, 1, 2 and that the embedding of Y 0 into Y 1 is compact. Let 1 ă r ă 8.
Then, the space Integrating the above estimate from 0 to τ and using the fact that

W " L 8 pr0, T s; Y 0 q X W 1,r pr0, T s; Y 2 q is compactly embedded in Cpr0, T s; Y 1 q. Proof.
ˇˇxBuptq, A 1´2α 0 zptqy ´α,α ˇˇď }Buptq} ´α › › A 1´α 0 zptq › › , it follows that }A 1 2 ´α 0 zpτ q} 2 `ż τ 0 }A 1´α 0 zptq} 2 dt ď }A 1 2 ´α 0 z 0 } 2 `}Bu} 2 L 2 pr0,τ s;X ´αq .
This leads to:

}zpτ q} 2 1 2 ´α `ż τ 0 }zpsq} 2 1´α ds ď ż τ 0 }Bupsq} 2 ´αds `}z 0 } 2 1 2
´α.

(4.2)

Similarly, we have: and integrating by time from 0 to τ , we obtain the following inequality:

x
}A 1 2 ´α 0 zpτ q} 2 `ż τ 0 }A ´α 0 9 zptq} 2 dt ď }A 1 2 ´α 0 z 0 } 2 `}Bu} 2 L 2 pr0,τ s;X ´αq .
This leads to:

}zpτ q} 2 1 2 ´α `ż τ 0 } 9 zpsq} 2 ´αds ď ż τ 0 }Bupsq} 2 ´αds `}z 0 } 2 1 2
´α.

(

Then (4.1) is deduced from (4.2) and (4.3). This ends the proof.

We give here some convergence results which play an important rule in the proof of our main results.

Lemma 4.3. Let z 0 P X, T ą 0, ū P L 8 pr0, T s; U q and τ P r0, T s. Let pu n q n be a sequence of L 8 pr0, T s; U q and pτ n q n be a sequence of r0, T s such that u n P L 8 pr0, T s; U hn q for any n P N and

u n Ñ ū weakly* in L 8 pr0, T s; U q, τ n Ñ τ in r0, T s.
Then, the following convergence results hold 1. lim nÑ`8 }zp., z 0 , u n q ´zp., z 0 , ūq} Cpr0,T s;Xq " 0, 2. lim nÑ`8 }zpτ n , z 0 , u n q ´zpτ , z 0 , u n q} " 0.

Proof.

1. Denote ψptq " zpt, z 0 , ūq and ψ n ptq " zpt, z 0 , u n q. By the standard energy estimate (4.1), we know that pψ n q n is a bounded sequence in

W " Cpr0, T s; X 1 2 ´αq X L 2 pr0, T s; X 1´α q X W 1,2 pr0, T s; X ´αq. (4.4)
Using Theorem 4.1 with Y 0 " X 1 2 ´α, Y 1 " X and Y 2 " X ´α, we deduce that:

D r ψ P Cpr0, T s; Xq s.t. ψ n Ñ ψ strongly in Cpr0, T s; Xq and ψ n Ñ ψ weakly in W.

Now we prove that r ψ " ψ. We know that ψ n satisfies:

9 ψ n `A0 ψ n " Bu n , ψ n p0q " z 0 .
Then it is clear that 9 ψ n Ñ 9 r ψ weakly in L 2 pr0, T s; X ´αq and A 0 ψ n Ñ A 0 r ψ weakly in L 2 pr0, T s; X ´αq since ψ n Ñ ψ weakly in W . Moreover, u n á ū weakly* in L 8 pr0, T s; U q implies that u n á ū weakly in L 2 pr0, T s; U q. Thus, Bu n á B ū weakly in L 2 pr0, T s; Xq. Finally, ψp0q " z 0 . Indeed, we know that z 0 " ψ n p0q Ñ ψp0q, since ψ n Ñ ψ strongly in Cpr0, T s; Xq.

Consequently, ψ satisfies:

9 ψ `A0 ψ " B ū, ψp0q " z 0 ,
which implies that ψ " ψ. This leads to the first assertion.

We first notice that

}zpτ n , z 0 , u n q ´zpτ , z 0 , u n q} ď }T τn z 0 ´Tτ z 0 } `› › › › ż τn 0 T τn´σ Bu n pσqdσ ´ż τ 0 T τ ´σBu n pσqdσ › › › › ď }T τn z 0 ´Tτ z 0 } `› › › › › ż Ă τn 0 T τn´σ Bu n pσqdσ ´ż Ă τn 0 T τ ´σBu n pσqdσ › › › › › `› › › › ż τn Ă τn T τn´σ Bu n pσqdσ › › › › `› › › › ż τ Ă τn T τ ´σBu n pσqdσ › › › › " L 1 `L2 `L3 `L4 , where r τ n " mintτ n , τ u.
By using the continuity of the map t Þ Ñ T t z 0 , it is clear that L 1 converges to zero when n tends to `8.

Moreover, since pu n q n Ă L 8 pr0, T s, U q and using (3.6), we have

L 3 ď ż τn Ă τn }T τn´σ Bu n pσq} dσ ď ż τn Ă τn }T τn´σ A α 0 } LpXq › › A ´α 0 B › › LpU,Xq }u n pσq} U dσ ď C ż τn Ă τn
pτ n ´σq ´αdσ " C 1 ´α pτ n ´r τ n q 1´α Ñ 0, since 1 ´α ą 0. Similarly, L 4 converges to zero. Now we prove that L 2 converges to zero. We have

L 2 " › › › › › ż Ă τn 0 pT τn´σ ´Tτ´σ q B pu n pσq ´ūpσqq dσ › › › › › `› › › › › ż Ă τn 0 pT τn´σ ´Tτ´σ q B ūpσqdσ › › › › › " L 2,1 `L2,2 .
By using the continuity of the map t Þ Ñ T t z 0 , it is clear that L 2,2 converges to zero. Now we prove that L 2,1 tends to 0 when n tends to `8.

We have

L 2,1 ď › › › › › ż Ă τn 0 T τn´σ B pu n pσq ´ūpσqq dσ › › › › › `› › › › › ż Ă τn 0 T τ ´σB pu n pσq ´ūpσqq dσ › › › › › " }T τn´Ă τn }}zp r
τ n , z 0 , u n q ´zp r τ n , z 0 , ūq} `}T τ ´Ă τn }}zp r τ n , z 0 , u n q ´zp r τ n , z 0 , ūq}.

Since }T τn´Ă τn } and }T τ ´Ă τn } are bounded, by using the first assertion of this lemma, we can deduce that L 2,1 converge to zero when n tends to `8. This ends the proof of the second assertion.

We then give some properties which will be used to prove the convergence results:

Lemma 4.4.

Denote λ 1 (resp. λ 1,h ) the first eigenvalue of the operator A 0 (resp. A 0,h ). For every z 0 P X and z 0,h P V h the following properties hold 1. τ 0 pz 0 q ď lnp}z 0 }{εq λ 1 .

2. τ h pz 0,h q ď lnp}z 0,h }{εq λ 1,h .

λ 1 ď λ 1,h .

Proof. To prove assertion 1. we notice that }zpt, z 0 , 0q} ď e ´λ1 t }z 0 } and taking t " lnp}z 0 }{εq λ 1 we obtain that }zpt, z 0 , 0q} ď ε. This proves 1. By a similar argument, it is clear that assertion 2. holds.

We end by proving the last assertion. In fact, this inequality is easily deduced by the min-max formula:

λ 1 " min zPX 1 2 }A 1 2 0 z} 2 }z} 2 and λ 1,h " min zPV h }A 1 2 0 z} 2 }z} 2 .
This ends the proof.

Given ε ą 0 and z 0 P X such that }z 0 } ą ε, the lemma below shows that τ h has a strictly positive lower and upper bounds, which are independent of h. Lemma 4.5. Assume that pC1q holds. For every z 0 P X such that }z 0 } ą ε, there exist c, C ą 0, h ą 0 such that for any h P p0, hq, we have c ď τ h pP h z 0 q ď C, where C " 2 lnp}z 0 }{εq λ 1 .

Proof. We begin by proving that τ h pP h z 0 q is bounded from below. Assume by contradiction that lim hÑ0 τ h pP h z 0 q " 0.

We first notice that }z h pτ h , P h z 0 , u hq ´zh p0, P h z 0 , u hq} "

› › › › › T τ h ,h P h z 0 `ż τ h 0 T τ h ´σ,h Bu hpσqdσ ´Ph z 0 › › › › › ď › › › ´Tτ h ,h ´Id ¯Ph z 0 › › › `› › › › › ż τ h 0 T τ h ´σ,h Bu hpσqdσ › › › › › .
We can apply the Trotter-Kato theorem (see, for instance, [2, Theorem 3.6.1, Proposition 3.6.2]), to get that the first term of the right hand side of the previous inequality tends to zero when h tends to zero. Moreover, concerning the left hand side, using (3.7), we have

› › › › › ż τ h 0 T τ h ´σ,h Bu hpσqdσ › › › › › ď ż τ h 0 › › ›T τ h ´σ,h A α h › › › LpXq › › A ´α h B › › LpU,Xq }u h pσq} U dσ ď C ż τ h 0 pτ h ´σq ´αdσ " C 1 ´α pτ h q 1´α Ñ 0, since 1 ´α ą 0.
Consequently, we have lim hÑ0 }z h pτ h , P h z 0 , u hq ´zh p0, P h z 0 , u hq} " lim hÑ0 }z h pτ h , P h z 0 , u hq ´Ph z 0 } " 0.

Using the fact that }z h pτ h , P h z 0 , u hq} ď ε, it is clear that lim hÑ0 }P h z 0 } ď ε. However, note that lim hÑ0 }P h z 0 ´z0 } " 0, which comes from pC1q since X γ is dense in X for any γ ą 0 and P h is a projection. This leads to the contradiction with the fact that }z 0 } ą ε.

We prove now that τ h pP h z 0 q is bounded from above. This is obvious by using Lemma 4.4, since

τ h pP h z 0 q ď lnp}P h z 0 }{εq λ 1,h ď 2 lnp}z 0 }{εq λ 1 ă `8.
Proof of Theorem 1.1. It suffices to prove the following two inequalities: lim inf hÑ0 τ h pP h z 0 q ě τ 0 pz 0 q, (4.5) lim sup hÑ0 τ h pP h z 0 q ď τ 0 pz 0 q. (4.6)

We begin by proving (4.5). In the following, we denote simply τ 0 pz 0 q by τ 0 and τ h pP h z 0 q by τ h when there is no ambiguity. We notice at first the following property: @ T ą 0, @ u P U ad , τ 0 pz 0 q ď T `τ 0 pzpT, z 0 , uqq.

We deduce from Theorem 3.6 that }zpτ h , z 0 , u hq ´zh pτ h , P h z 0 , u hq} ď Ch θ τ h ´1}z 0 } X `Ch θp1´αq p| ln h| `1 `τ h q }u h} L 8 pr0,`8q;U h q , which leads to: }zpτ h , z 0 , u hq} ď ε `Ch θ τ h ´1}z 0 } X `Ch θp1´αq p| ln h| `1 `τ h q }u h} L 8 pr0,`8q;U h q ď ε `Ch θ τ h ´1 `Ch θp1´αq p| ln h| `1 `τ h q.

We denote r z 0 " zpτ h , z 0 , u hq. According to (4.7) with T " τ h pP h z 0 q, we have: τ 0 pz 0 q ď τ h pP h z 0 q `τ 0 p r z 0 q.

Since u h P L 8 pr0, `8q; U h q Ă L 8 pr0, `8q; U q and }u hptq} ď 1, u h is then an admissible control for the continuous system. Then, according to Lemma 4.4, we have:

τ 0 ď τ h `lnppε `Ch θ τ h ´1 `Ch θp1´αq p| ln h| `1 `τ h qq{εq λ 1 ď τ h `Ch θ λ 1 ε τ h ´1 `C h θp1´αq λ 1 ε p1 `τ h q `C h θp1´αq | ln h| λ 1 ε . ( 4.8)
Thus, (4.5) can be deduced by taking h to zero and by the fact that C ą lim hÑ0 τ h pP h z 0 q ą c ą 0 (see Lemma 4.5).

We now prove the second inequality (4.6). We have: }z h pτ 0 , P h z 0 , Q h u ˚q ´zpτ 0 , z 0 , u ˚q} ď }z h pτ 0 , P h z 0 , Q h u ˚q ´zh pτ 0 , P h z 0 , u ˚q} `}z h pτ 0 , P h z 0 , u ˚q ´zpτ 0 , z 0 , u ˚q} ď }z h pτ 0 , P h z 0 , Q h u ˚q ´zh pτ 0 , P h z 0 , u ˚q} `Ch θ τ 0 ´1 `Ch θp1´αq p| ln h| `1 `τ 0 q, using (3.9). Denote f phq " }z h pτ 0 , P h z 0 , Q h u ˚q ´zh pτ 0 , P h z 0 , u ˚q}. Notice that lim hÑ0 f phq " 0. Indeed,

lim hÑ0 f phq " lim hÑ0 }z h pτ 0 , P h z 0 , Q h u ˚q ´zh pτ 0 , P h z 0 , u ˚q} " lim hÑ0 }Φ τ ˚,h pu ˚´Q h u ˚q}.
Moreover, we have

}Φ τ ˚,h pu ˚´Q h u ˚q} ď } `Φτ ˚,h ´Φτ ˚˘pu ˚´Q h u ˚q} `}Φ τ ˚pu ˚´Q h u ˚q}.
Since B P LpU, X ´αq, with 0 ď α ă 1 2 , it is known that Φ τ ˚P LpL 2 pr0, τ ˚s; U q, Xq (see for instance [START_REF] Tucsnak | Observation and control for operator semigroups[END_REF]Proposition 5.1.3]). Combining this fact with assertion 7. in Proposition 3.5 and with assumption pC5q, it is clear that lim hÑ0 f phq " 0.

Thus, we have:

}z h pτ 0 , P h z 0 , Q h u ˚q} ď ε `f phq `Ch θ τ 0 ´1 `Ch θp1´αq p| ln h| `1 `τ 0 q.
By the similar argument as in (4.8), we have: τ h pP h z 0 q ď τ 0 `lnppε `f phq `Ch θ τ 0 ´1 `Ch θp1´αq p| ln h| `1 `τ 0 qq{εq λ 1,h ď τ 0 `f phq `Ch θ τ 0 ´1 `Cp1 `τ 0 qh θp1´αq `Ch θp1´αq | ln h| λ 1 ε .

This leads to inequality (4.6) by taking h to zero.

Proof of Theorem 1.2. Denote T " 2 lnp}z 0 }{εq λ 1

. It is clear that τ h ď T for all h ą 0 and τ 0 ď T (see Lemma 4.4 and Lemma 4.5). We extend pu hq h and u 0 to time T by zero.

Since }u h} L 8 pr0,T s;U h q ď 1, there exist a control ū P L 8 pr0, T s; U q and a subsequence ph n q n Ñ 0, such that: u hn Ñ ū weakly* in L 8 pr0, T s; U q. Now we prove that ū " u 0 .

We then consider the standard P 1 finite element method with regular triangulation T h of Ω (see for example, Thomée [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF], Raviart and Thomas [START_REF] Raviart | Introduction à l'analyse numérique des équations aux dérivées partielles[END_REF]) and we build the finite element space V h Ă H 1 0 pΩq of X, defined by V h " ϕ P Cp Ωq |ϕ |T P P 1 pT q for every T P T h , ϕ |BΩ " 0 ( where P 1 pT q is the set of affine functions on T . We also define U h a finite dimensional subspace of U by U h :" tχ ω ϕ h | ϕ h P V h u .

Denote P h the L 2 ´projection from L 2 pΩq to V h , and we consider the following space semi-discrete scheme formulated in the form of (1.5)-(1.6): 9 z h ptq `Ah z h ptq " B h u h ptq pt ě 0q, (5.4) z h p0q " P h z 0 pz 0 P Xq, (

where u h is taken in L 8 pr0, `8q; U h q and the operator A h and B h are defined by

ă A h ϕ, ψ ą " ż Ω ∇ϕ∇ψ dx pϕ, ψ P V h q, (5.6)
and B h P LpU, V h q is defined by : xB h u, ϕy " xu, B ˚ϕy U pϕ P V h , u P U q.

(5.7)

Given ε ą 0 with the same notation as in Section 1, we consider the time optimal control problem pTPq associated to (5.1)-(5.3) (resp. pTP h q associated to (5.4)-(5.5)) and denote pτ 0 , u 0 q the corresponding optimal time and time-optimal control (resp. pτ h , u hq).

We can now state the following convergence results:

Proposition 5.2. With the notation above, let, for every ε ą 0, z 0 P X such that }z 0 } ą ε. Then we have that τ h Ñ τ 0 , and u h Ñ u 0 in L 2 pr0, T s; L 2 pΩqq.

Proof. In order to apply Theorem 1.1 and Theorem 1.2, we verify that conditions pC1q ´pC5q are satisfied with α " 0 and θ " 2. Indeed, pC1q is a standard error estimate when applying the P 1 finite element scheme (see for example [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF]). pC2q is proved in Bramble et al [START_REF] Bramble | Some convergence estimates for semidiscrete Galerkin type approximations for parabolic equations[END_REF] and in Lasiecka [START_REF] Lasiecka | Convergence estimates for semidiscrete approximations of nonselfadjoint parabolic equations[END_REF]. Moreover, pC3q´pC5q clearly hold true because here the control operator B is bounded. Furthermore, it is known that the assumption made in Theorem 1.2 holds true (see for example [START_REF] Micu | Time optimal boundary controls for the heat equation[END_REF]). Thus, we can apply Theorem 1.

  See for example [26, Cor. 4, p.85].

	it follows that if z satisfies the equation (1.1)-(1.2), then
				1 2	d dt	}A	1 2 0	´α	zptq} 2 `}A 1´α 0	zptq} 2 " xBuptq, A 1´2α 0	zptqy ´α,α .
	Lemma 4.2. Assume that z 0 P X 1 2	´α. Then, there exists C ą 0 such that
	}zpτ q} 2 1 2	´α	`ż τ 0	`} 9 zpsq} 2 ´α `}zpsq} 2 1´α ˘ds ď C	ˆż τ 0	}Bupsq} 2 ´αds `}z 0 } 2 1 2	´α˙p τ ą 0q,	(4.1)
	where z is the solution of (1.1)-(1.2).
	Proof. Using the facts that
	x 9 zptq, A 1´2α 0	zptqy ´α,α "	1 2	d dt	}A	1 2 0	´α	zptq} 2 ,	xA 0 zptq, A 1´2α 0	zptqy ´α,α " }A 1´α 0	zptq} 2 ,

  1 and Theorem 1.2 to conclude. Example 5.3. Let Ω Ă R N be a bounded open set with a boundary BΩ of class C 2 .

				We consider the
	heat equation
	9 zpx, tq ´ zpx, tq " 0	px P Ω, t ě 0q,	(5.8)
	with the initial and boundary conditions
	Bz Bν	|

BΩ px, tq " upx, tq px P BΩ, t ě 0q, (5.9) zpx, 0q " z 0 pxq P L 2 pΩq px P Ωq.

(5.10)

and with bounded operator B, which can not apply to boundary control problem. Our results

The main step here is to prove the following convergence property: }z hn pτ hn , P hn z 0 , u hn q ´zpτ 0 , z 0 , ūq} Ñ 0.

(4.9) Indeed, since Bp0, εq is complete (notice that Ď B h p0, εq Ă Bp0, εq), (4.9) leads to zpτ 0 , z 0 , ūq P Bp0, εq. Then, by the uniqueness of the time optimal control and Theorem 1.1, we deduce that ū " u 0 .

We are then reduced to prove (4.9). We have:

}z hn pτ hn , P hn z 0 , u hn q ´zpτ 0 , z 0 , ūq} ď }z hn pτ hn , P hn z 0 , u hn q ´zpτ hn , z 0 , u hn q} `}zpτ hn , z 0 , u hn q ´zpτ 0 , z 0 , u hn q} `}zpτ 0 , z 0 , u hn q ´zpτ 0 , z 0 , ūq}. (4.10)

The first term of the right hand side of (4.10) converges to zero using the error estimate (3.9). Moreover, by using the first and the second assertions in Lemma 4.3, it is clear that the second and the third terms of the right hand of (4.10) converges to zero. We then deduce (4.9). Thus, we have:

We deduce immediately that:

At last, since both u hn and u 0 are bang-bang controls (see Corollary 2.4), we have lim nÑ`8

}u hn } L 2 pr0,T s;U q " }u 0 } L 2 pr0,T s;U q .

According to the Radon-Riesz property, this leads to the strong convergence in L 2 pr0, T s; U q. This ends the proof. where ω Ă Ω is an open and non-empty subset with its characteristic function χ ω . The system (5.1)-(5.3) can be formulated in the abstract form (1.1)-(1.2) by taking A 0 " ´ with Dirichlet boundary conditions, of domain DpA 0 q " H 1 0 pΩq X H 2 pΩq on X " L 2 pΩq. The control operator B P LpU, Xq (here α " 0) is defined by

Bϕ " χ ω φ, where U " L 2 pωq and where φ is the extension by 0 of ϕ on Ω outside ω.

It is known that the system (5.8)-(5.10) can be written in the form of (1.1)-(1.2) by taking X " L 2 pΩq, U " L 2 pBΩq, DpA 0 q " z P H 2 pΩq | Bz Bν | BΩ " 0 ( , A 0 " ´ and the control operator B P LpU, X ´αq (here α " 1 4 `ε, ε ą 0) defined by

where N is the Neumann map. This map is defined by N v " z, where z P L 2 pΩq is the unique solution of the nonhomogeneous elliptic equation

(5.11)

We refer to [START_REF] Lasiecka | Galerkin approximations of abstract parabolic boundary value problems with rough boundary data-l p theory[END_REF]Par. 5] which shows that the Galerkin semi-discrete approximation of system (5.8)-(5.10) satisfies conditions pC1q ´pC5q (and for other numerical approximation method).

We can obtain the following result:

Proposition 5.4. With the notation above, let, for every ε ą 0, z 0 P X such that }z 0 } ą ε. Then we have that τ h Ñ τ 0 , and

Proof. In order to apply Theorem 1.1 and Theorem 1.2, we verify that conditions pC1q ´pC5q are satisfied with α " 1 4 `ε, ε ą 0 and θ " 2 (see for example [START_REF] Lasiecka | Galerkin approximations of abstract parabolic boundary value problems with rough boundary data-l p theory[END_REF]Ch. 5]). Moreover, it is known that the assumption made in Theorem 1.2 holds true (see for example [START_REF] Lasiecka | Control Theory for Partial Differential Equations: Continuous and Approximation Theories[END_REF]Section 3.3,Ch. 5.3]). Thus, we can apply Theorem 1.1 and Theorem 1.2 to conclude.

We illustrate in the follows some numerical results for the 1d-heat equation with Neumann boundary condition. More precisely, we take Ω " p0, 1q, " 1 (the radius of the target ball), control constraints |upt, 0q| ď 1 and |upt, 1q| ď 1 and the initial data z 0 pxq " 5sinpπxq. We choose to discretize totally the system, with the implicit Euler method for time and finite difference scheme for space. We also use the optimization solver IPOPT in Matlab.

We obtain the following result: Moreover, we find that u h pt, 0q " u h pt, 1q " ´1 for any t ą 0 and for any number of discretization, that validates the fact that the time optimal control for the discrete system converges towards a bang-bang time optimal control. Remark 5.5. The error estimates of the time and control convergence remain to study.


Hiep H Nguyen
email: huu-hiep.nguyen@inria.fr

Abdessamad Imine
email: abdessamad.imine@loria.fr

Michaël Rusinowitch
email: michael.rusinowitch@inria.fr

Detecting Communities under Differential Privacy

Keywords: Differential privacy, Community detection, LouvainDP, ModDivisive

Complex networks usually expose community structure with groups of nodes sharing many links with the other nodes in the same group and relatively few with the nodes of the rest. This feature captures valuable information about the organization and even the evolution of the network. Over the last decade, a great number of algorithms for community detection have been proposed to deal with the increasingly complex networks. However, the problem of doing this in a private manner is rarely considered.

In this paper, we solve this problem under differential privacy, a prominent privacy concept for releasing private data. We analyze the major challenges behind the problem and propose several schemes to tackle them from two perspectives: input perturbation and algorithm perturbation. We choose Louvain method as the back-end community detection for input perturbation schemes and propose the method LouvainDP which runs Louvain algorithm on a noisy supergraph. For algorithm perturbation, we design ModDivisive using exponential mechanism with the modularity as the score. We have thoroughly evaluated our techniques on real graphs of different sizes and verified that ModDivisive steadily gives the best modularity and avg.F1Score on large graphs while LouvainDP outperforms the remaining input perturbation competitors in certain settings.

INTRODUCTION

Graphs represent a rich class of data observed in daily life where entities are described by nodes and their connections are characterized by edges. Apart from microscopic (node level) and macroscopic (graph level) configurations, many complex networks display a mesoscopic structure, i.e. they appear as a combination of components fairly independent of each other. These components are called communities, ACM ISBN 978-1-4503-2138-9. DOI: 10.1145/1235 modules or clusters and the problem of how to reveal them plays a significant role in understanding the organization and function of complex networks. Over the last decade, a great number of algorithms for community detection (CD) have been proposed to address the problem in a variety of settings, such as undirected/directed, unweighted/weighted networks and non-overlapping/overlapping communities (for a comprehensive survey, see [START_REF] Fortunato | Community detection in graphs[END_REF]).

These approaches, however, are adopted in a non-private manner, i.e. a data collector (such as Facebook) knows all the contributing users and their relationships before running CD algorithms. The output of such a CD, in the simplest form, is a clustering of nodes. Even in this case, i.e. where only a node clustering (not the whole graph) is revealed, contributing users privacy may still be put at risk. For example, thirteen nodes in Fig. 1 are clustered in four communities [{0, 1, 2}, {3, 4}, {5, 6, 11, 12}, {7, 8, 9, 10}]. Assuming that a new edge (0,3) is added, the new clustering of nodes may become [{0, 1, 2, 3}, {4, 5, 6, 11, 12}, {7, 8, 9, 10}]. By observing the clusterings before and after the addition of the edge (0,3), an attacker infers that there might have been new edge(s) between 3 and the nodes in the community {0,1,2}.

In this paper, we address the problem of CD from the perspective of differential privacy [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF]. This privacy model offers a formal definition of privacy with a lot of interesting properties: no computational/informational assumptions about attackers, data type-agnosticity, composability and so on [START_REF] Mcsherry | Privacy integrated queries: an extensible platform for privacy-preserving data analysis[END_REF]. By differential privacy, we want to ensure the existence of connections between users to be hidden in the output clustering while keeping the low distortion of clusters compared to the ones generated by the corresponding non-private algorithms.

As far as we know, the problem is quite new and only mentioned in the recent work [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF] where Mülle et al. use a randomized response technique [START_REF] Dwork | The algorithmic foundations of differential privacy[END_REF] to perturb the input graph so that it satisfies differential privacy before running conventional CD algorithms. This scheme (we call it EdgeFlip afterwards) is classified as input perturbation in differential privacy literature (the other two categories are algorithm perturbation and output perturbation). Similarly, TmF approach [START_REF] Nguyen | Differentially private publication of social graphs at linear cost[END_REF] can apply to the true graph to get noisy output graphs as in the work of Mülle et al. Earlier, 1k-Series [START_REF] Wang | Preserving differential privacy in degree-correlation based graph generation[END_REF], Density Explore Reconstruct (DER) [START_REF] Chen | Correlated network data publication via differential privacy[END_REF] and HRG-MCMC [START_REF] Xiao | Differentially private network data release via structural inference[END_REF] are the best known methods for graph structure release under differential privacy. These methods can be followed by any exact CD algorithm to get a noisy clustering satisfying differential privacy. We choose Louvain method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] as such a CD algorithm. However, as we will see in the exper-iments, the output clusterings by the aforementioned methods have very low modularity scores. This fact necessitates new methods for CD problem under differential privacy.

Our main contributions are the new schemes LouvainDP (input perturbation) and ModDivisive (algorithm perturbation) which perform much better than the state-of-theart. LouvainDP is a high-pass filtering method that randomly groups nodes into supernodes of equal size to build a weighted supergraph. LouvainDP is guaranteed to run in linear time. ModDivisive is a top-down approach which privately divides the node set into the k-ary tree guided by the modularity score at each level. The main technique used in ModDivisive is the Markov Chain Monte Carlo (MCMC) to realize the exponential mechanism [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF]. We show that ModDivisive's runtime is linear in the number of nodes, the height of the binary tree and the burn-in factor of MCMC. The linear complexity enables us to examine million-scale graphs in a few minutes. The experiments show the high modularity and low distortion of the output clusters by Lou-vainDP and ModDivisive.

Our contributions are summarized as follows:

• We analyze the major challenges of community detection under differential privacy ( Section 3). We explain why techniques borrowed from k-Means fail and how the difficulty of ǫ-DP recommender systems justifies a relaxation of ǫ.

• We design an input perturbation scheme LouvainDP (Section 4) that runs in linear time using the high-pass filtering technique from [START_REF] Cormode | Differentially private summaries for sparse data[END_REF] and Louvain method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF].

• We propose an algorithm perturbation scheme Mod-Divisive (Section 5) as a divisive approach by using the modularity-based score function in the exponential mechanism. We prove that modularity has small global sensitivity and ModDivisive also runs in linear time.

• We conduct a thorough evaluation on real graphs of different sizes and show the outperformance of LouvainDP and ModDivisive over the state-of-the-art (Section 6).

Table 1 summarizes the key notations used in this paper.

RELATED WORK

Community Detection in Graphs

There is a vast literature on community detection in graphs. For a recent comprehensive survey, we refer to [START_REF] Fortunato | Community detection in graphs[END_REF]. In this section, we discuss several classes of techniques.

Newman and Girvan [START_REF] Newman | Finding and evaluating community structure in networks[END_REF] propose modularity as a quality of network clustering. It is based on the idea that a random graph is not expected to have a modular structure, so the possible existence of clusters is revealed by the comparison between the actual density of edges in a subgraph and the density one would expect to have in the subgraph if the nodes of the graph were connected randomly (the null model). The modularity Q is defined as where nc is the number of clusters, lc is the total number of edges joining nodes in community c and dc is the sum of the degrees of the nodes of c.

Q = nc c=1 lc m - dc 2m 2 (1)
Symbol Definition G = (V, E G ) true graph with n = |V | and m = |E G | G ′ = (V, E G ′ ) neighboring graph of G G = (V, E G) sample noisy output graph G 1 = (V 1 , E
Many methods for optimizing the modularity have been proposed over the last ten years, such as agglomerative greedy [START_REF] Clauset | Finding community structure in very large networks[END_REF], simulated annealing [START_REF] Medus | Detection of community structures in networks via global optimization[END_REF], random walks [START_REF] Pons | Computing communities in large networks using random walks[END_REF], statistical mechanics [START_REF] Reichardt | Statistical mechanics of community detection[END_REF], label propagation [START_REF] Raghavan | Near linear time algorithm to detect community structures in large-scale networks[END_REF] or InfoMap [START_REF] Rosvall | Maps of random walks on complex networks reveal community structure[END_REF], just to name a few. The recent multilevel approach, also called Louvain method, by Blondel et al. [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] is a top performance scheme. It scales very well to graphs with hundreds of millions of nodes. This is the chosen method for the input perturbation schemes considered in this paper (see Sections 3.2 and 4 for more detail). By maximizing the modularity, Louvain method is based on edge counting metrics, so it fits well with the concept of edge differential privacy (Section 3.1). One of the most recent methods SCD [START_REF] Prat-Pérez | High quality, scalable and parallel community detection for large real graphs[END_REF] is not chosen because it is about maximizing Weighted Community Clustering (WCC) instead of the modularity. WCC is based on triangle counting which has high global sensitivity (up to O(n)) [START_REF] Zhang | Private release of graph statistics using ladder functions[END_REF]. Moreover, SCD pre-processes the graph by removing all edges that do not close any triangle. This means all 1-degree nodes are excluded and form singleton clusters. The number of output clusters is empirically up to O(n).

Graph Release via Differential Privacy

In principle, after releasing a graph satisfying ǫ-DP, we can do any mining operations on it, including community detection. The research community, therefore, expresses a strong interest in the problem of graph release via differential privacy. Differentially private algorithms relate the amount of noise to the computation sensitivity. Lower sensitivity implies smaller added noise. Because the edges in simple undirected graphs are usually assumed to be independent, the standard Laplace mechanism [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] is applicable (e.g. adding Laplace noise to each cell of the adjacency matrix). However, this approach severely deteriorates the graph structure.

The state-of-the-art [START_REF] Chen | Correlated network data publication via differential privacy[END_REF][START_REF] Xiao | Differentially private network data release via structural inference[END_REF] try to reduce the sensitivity of the graph in different ways. Density Explore Reconstruct (DER) [START_REF] Chen | Correlated network data publication via differential privacy[END_REF] employs a data-dependent quadtree to summarize the adjacency matrix into a counting tree and then reconstructs noisy sample graphs. DER is an instance of input perturbation. Xiao et al. [START_REF] Xiao | Differentially private network data release via structural inference[END_REF] propose to use Hierarchical Random Graph (HRG) [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF] to encode graph structural information in terms of edge probabilities. Their scheme HRG-MCMC is argued to be able to sample good HRG models which reflect the community structure in the original graph. HRG-MCMC is classified as an algorithm perturbation scheme. A common disadvantage of the state-of-the-art DER and HRG-MCMC is the scalability issue. Both of them incur quadratic complexity O(n 2 ), limiting themselves to medium-sized graphs.

Recently, Nguyen et al. [START_REF] Nguyen | Differentially private publication of social graphs at linear cost[END_REF] proposed TmF that utilizes a filtering technique to keep the runtime linear in the number of edges. TmF also proves the upper bound O(ln n) for the privacy budget ǫ. At the same time, Mülle et al. [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF] devised the scheme EdgeFlip using a randomized response technique. However, EdgeFlip costs O(n 2 ) and is runnable only on graphs of tens of thousands of nodes.

PRELIMINARIES

In this section, we review key concepts and mechanisms of differential privacy.

Differential Privacy

Essentially, ǫ-differential privacy (ǫ-DP) [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] is proposed to quantify the notion of indistinguishability of neighboring databases. In the context of graph release, two graphs G1 = (V1, E1) and G2 = (V2, E2) are neighbors if V1 = V2, E1 ⊂ E2 and |E2| = |E1| + 1. Note that this notion of neighborhood is called edge differential privacy in contrast with the notion of node differential privacy which allows the addition of one node and its adjacent edges. Our work follows the common use of edge differential privacy. The formal definition of ǫ-DP for graph data is as follows.

Definition 3.1. A mechanism A is ǫ-differentially private if for any two neighboring graphs G1 and G2, and for any output

O ∈ Range(A), P r[A(G1) ∈ O] ≤ e ǫ P r[A(G2) ∈ O]
Laplace mechanism [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF] and Exponential mechanism [START_REF] Mcsherry | Privacy integrated queries: an extensible platform for privacy-preserving data analysis[END_REF] are two standard techniques in differential privacy. The latter is a generalization of the former. Laplace mechanism is based on the concept of global sensitivity of a function f which is defined as ∆f = maxG 1 ,G 2 ||f (G1)-f (G2)||1 where the maximum is taken over all pairs of neighboring G1, G2. Given a function f and a privacy budget ǫ, the noise is drawn from a Laplace distribution Lap(λ) : p(x|λ) = 1 2λ e -|x|/λ where λ = ∆f /ǫ. Theorem 3.1. (Laplace mechanism [START_REF] Dwork | Calibrating noise to sensitivity in private data analysis[END_REF]) For any function

f : G → R d , the mechanism A A(G) = f (G) + Lap1( ∆f ǫ ), ..., Lap d ( ∆f ǫ ) (2)
satisfies ǫ-differential privacy, where Lapi( ∆f ǫ ) are i.i.d Laplace variables with scale parameter ∆f ǫ . 2

Geometric mechanism [START_REF] Ghosh | Universally utility-maximizing privacy mechanisms[END_REF] is a discrete variant of Laplace mechanism with integral output range Z and random noise ∆ generated from a two-sided geometric distribution Geom(α) : P r[∆ = δ|α] = 1-α 1+α α |δ| . To satisfy ǫ-DP, we set α = exp(-ǫ). We use geometric mechanism in our Lou-vainDP scheme.

For non-numeric data, the exponential mechanism is a better choice [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF]. Its main idea is based on sampling an output O from the output space O using a score function u. This function assigns exponentially higher probabilities to outputs of higher scores. Let the global sensitivity of u be ∆u = maxO,G 1 ,G 2 |u(G1, O)u(G2, O)|. [START_REF] Mcsherry | Mechanism design via differential privacy[END_REF]) Given a score function u : (G×O) → R for a graph G, the mechanism A that samples an output O with probability proportional to exp( ǫ.u(G,O) 2∆u ) satisfies ǫ-differential privacy. 2

Composability is a nice property of differential privacy which is not satisfied by other privacy models such as kanonymity.

Theorem 3.3. (Sequential and parallel compositions [START_REF] Mcsherry | Privacy integrated queries: an extensible platform for privacy-preserving data analysis[END_REF]) Let each Ai provide ǫi-differential privacy. A sequence of Ai(D) over the dataset D provides Σ n i=1 ǫi-differential privacy.

Let each Ai provide ǫi-differential privacy. Let Di be arbitrary disjoint subsets of the dataset D. The sequence of Ai(Di) provides max n i=1 ǫi-differential privacy. 2

Louvain Method

Since its introduction in 2008, Louvain method [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF] becomes one of the most cited methods for the community detection task. It optimizes the modularity by a bottom-up folding process. The algorithm is divided in passes each of which is composed of two phases that are repeated iteratively. Initially, each node is assigned to a different community. So, there will be as many communities as there are nodes in the first phase. Then, for each node i, the method considers the gain of modularity if we move i from its community to the community of a neighbor j (a local change). The node i is then placed in the community for which this gain is maximum and positive (if any), otherwise it stays in its original community. This process is applied repeatedly and sequentially for all nodes until no further improvement can be achieved and the first pass is then complete.

We demonstrate Louvain method in Fig. 1 by a graph of 13 nodes and 20 edges. If each node forms its own singleton community, the modularity Q will be -0.0825. In the first pass of Louvain method, each node moves to the best community selected from its neighbors' communities. We get the partition [{0, 1, 2}, {3, 4}, {5, 6, 11, 12}, {7, 8, 9, 10}] with modularity 0.46375. The second phase of first pass builds a weighted graph corresponding to the partition by aggregating communities. The second pass repeats the folding process on this weighted graphs to reach the final partition [{0, 1, 2}, {3, 4, 5, 6, 11, 12}, {7, 8, 9, 10}] with modularity 0.47.

This simple agglomerative algorithm has several advantages as stated in [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. First, its steps are intuitive and easy to implement, and the outcome is unsupervised. Second, the algorithm is extremely fast, i.e. computer simulations on large modular networks suggest that its complexity is linear on typical and sparse data. This is due to the fact that the possible gains in modularity are easy to compute and the number of communities decreases drastically after just a few passes so that most of the running time is concentrated on the first iterations. Third, the multi-level nature of the method produces a hierarchical structure of communities which allows multi-resolution analysis, i.e the user can zoom in the graph to observe its structure with the desired resolution. In addition, Louvain method is runnable on weighted graphs. This fact supports naturally our scheme LouvainDP as described in Section 4.1.

Challenges of Community Detection under Differential Privacy

In this section, we explain why community detection under differential privacy is challenging. We show how techniques borrowed from related problems fail. We also advocate the choice of ǫ as a function of the graph size n.

The problem of differentially private community detection is closely related to ǫ-DP k-Means clustering and recommender systems. The ǫ-DP k-Means is thoroughly discussed in [START_REF] Su | Differentially private k-means clustering[END_REF]. However, techniques from ǫ-DP k-Means are not suitable to ǫ-DP community detection. First, items in k-Means are in low-dimensional spaces and the number of clusters k is usually small. This contrast to the case of community detection where nodes lie in a n-dimensional space and the number of communities varies from tens to tens of thousands, not to say the communities may overlap or be nested (multi-scale). Second, items in ǫ-DP k-Means are normalized to [-1, 1] d while the same preprocessing seems invalid in ǫ-DP community detection. Moreover, the output of k-Means usually consists of equal-sized balls while this is not true for communities in graphs. Considering the graph as a high-dimensional dataset, we tried the private projection technique in [START_REF] Kenthapadi | Privacy via the johnson-lindenstrauss transform[END_REF] which is followed by spectral clustering, but the modularity scores of the output are not better than random clustering.

Recent papers on ǫ-DP recommender systems [START_REF] Guerraoui | D2p: distance-based differential privacy in recommenders[END_REF][START_REF] Banerjee | The price of privacy in untrusted recommender systems. Selected Topics in Signal Processing[END_REF] show that privately learning the clustering of items from user ratings is hard unless we relax the value of ǫ up to log n. Banerjee et al. [START_REF] Banerjee | The price of privacy in untrusted recommender systems. Selected Topics in Signal Processing[END_REF] model differentially private mechanisms as noisy channels and bound the mutual information between the generative sources and the privatized sketches. They show that in the information-rich regime (each user rates O(n) items), their Pairwise-Preference succeeds if the number of users is Ω(n log n/ǫ). Compared to ǫ-DP community detection where the number of users is n, we should have ǫ = Ω(log n). Similarly, in D2P scheme, Guerraoui et al. [START_REF] Guerraoui | D2p: distance-based differential privacy in recommenders[END_REF] draw a formula for ǫ as

ǫ (p,0,λ) D2P = ln(1 + (1 -p).NE p.|G λ | ) ( 3
)
where λ is the distance used to conceal the user profiles (λ = 0 reduces to the classic notion of differential privacy). NE is the number of items which is exactly n in community detection. |G λ | is the minimum size of user profiles at distance λ over all users (|G λ | = o(NE) except at unreasonably large λ). Clearly, at p = 0.5 (as used in [START_REF] Guerraoui | D2p: distance-based differential privacy in recommenders[END_REF]), we have ǫ (0.5,0,λ) D2P

≈ ln n. The randomized technique in D2P is very similar to EdgeFlip [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF] which is shown ineffective in ǫ-DP community detection for ǫ ∈ (0, 0.5 ln n) (Section 6). Note that ǫ-DP community detection is unique in the sense that the set of items and the set of users are the same. Graphs for community detection are more general than bipartite graphs in recommender systems. In addition, modularity Q (c.f. To further emphasize the difficulty of ǫ-DP community detection, we found that IDC scheme [START_REF] Gupta | Iterative constructions and private data release[END_REF] using Sparse Vector Technique [9, Section 3.6] is hardly feasible. As shown in Algorithm 1 of [START_REF] Gupta | Iterative constructions and private data release[END_REF], to publish a noisy graph that can approximately answer all cut queries with bounded error m 0.25 n/ǫ 0.5 , IDC must have B(α) "yes" queries among all k queries. B(α) may be as low as √ m but k = 2 2n . In the average case, IDC incurs exponential time to complete. The typical epsilon in the literature is 1.0 or less. However, this value is only applicable to graph metrics of low sensitivity O(1) such as the number of edges, the degree sequence. The global sensitivity of other metrics like the diameter, the number of triangles, 2K-series etc. is O(n), calling for smooth sensitivity analysis (e.g. [START_REF] Nissim | Smooth sensitivity and sampling in private data analysis[END_REF]). For counting queries, Laplace/Geometric mechanisms are straightforward on real/integral (metric) spaces. However, direct noise adding mechanisms on the space P of all ways to partition the nodeset V are non-trivial because |P| ≈ n n and P is nonmetric. Compared to ǫ-DP recommender systems discussed above which use ǫ = ln n for the super-exponential spaces of size O(2 n(n-1)/2 ), ǫ-DP community detection clearly needs lower privacy budget.

To conclude, ǫ-DP community detection is challenging and requires new techniques. In this paper, we evaluate the schemes for ǫ up to 0.5 ln n. At ǫ = 0.5 ln n, the multiplicative ratio (c.f. Definition 3.1) is e ǫ = e 0.5 ln n = √ n. We believe it is a reasonable threshold for privacy protection compared to ǫ = ln n (i.e. e ǫ = n) in ǫ-DP graph release and ǫ-DP recommender systems.

INPUT PERTURBATION

In this section, we propose the linear scheme LouvainDP that uses a filtering technique to build a noisy weighted supergraph and calls the exact Louvain method subsequently. Then we discuss several recent ǫ-DP schemes that can be classified as input perturbation. Fig. 2a sketches the basic steps of the input perturbation paradigm.

LouvainDP: Louvain Method on Noisy Supergraphs

The basic idea of LouvainDP is to create a noisy weighted supergraph G1 from G by grouping nodes into supernodes of equal size k. We then apply the filtering technique of Cormode et al. [START_REF] Cormode | Differentially private summaries for sparse data[END_REF] to ensure that there are only O(m) noisy weighted edges in G1. Finally, we run the exact community detection on G1.

In [START_REF] Cormode | Differentially private summaries for sparse data[END_REF], Cormode et al. propose several summarization techniques for sparse data under differential privacy. Let M be a contingency table having the domain size m0 and m1 non-zero entries (m1 ≪ m0 for sparse data), the conventional publication of a noisy table M ′ from M that satisfies ǫ-DP requires the addition of Laplace/geometric noise to m0 entries. The entries in M ′ could be filtered (e.g. removing negative ones) and/or sampled to get a noisy summary M ′′ . This direct approach would be infeasible for huge domain sizes m0. Techniques in [START_REF] Cormode | Differentially private summaries for sparse data[END_REF] avoid materializing the vast noisy data by computing the summary M ′′ directly from M using filtering and sampling techniques.

In our LouvainDP, the supergraph G1 is an instance of sparse data with the domain size m0 = n 1 (n 1 +1) 2 where n1 is the number of supernodes and m1 is the number of non-zero entries corresponding to non-zero superedges. We use the one-sided filtering [START_REF] Cormode | Differentially private summaries for sparse data[END_REF] to efficiently compute G1 with O(m) edges in linear time.

Algorithm

LouvainDP can run with either geometric or Laplace noise. We describe the version with geometric noise in Algorithm 1.

Algorithm 1 LouvainDP(G, s)

Input: undirected graph G, group size k, privacy budget ǫ Output: noisy partition C

1:

G 1 ← ∅, n 1 = ⌊ |V | k ⌋ -1, V 1 ← {0, 1, .., n 1 } 2: ǫ 2 = 0.01, ǫ 1 = ǫ -ǫ 2 , α = exp(-ǫ 1 ) 3: get a random permutation Vp of V 4: compute the mapping M : Vp → V 1 5: compute superedges of G 1 : E 1 = {e 1 (i, j)} where i, j ∈ V 1 6: m 1 = |E 1 | + Lap(1/ǫ 2 ), m 0 = n 1 (n 1 +1) 2 7: θ = ⌈log α (1+α)m 1 m 0 -m 1 ⌉ 8: s = (m 0 -m 1 ) α θ 1+α
9: for e 1 (i, j) in E 1 do 10: e 1 (i, j) = e 1 (i, j) + Geom(α)

11:

if e 1 (i, j) ≥ θ then

12:

add e 1 (i, j) to G 1 13: for s edges e 0 (i, j) / ∈ E 1 sampled uniformly at random do 14:

draw w from the distribution P r[X ≤ x] = 1 -α x-θ+1

15:

if w > 0 then

16:

add edge e 0 (i, j) with weight w to G 1 17: run Louvain method on G 1 to get C1 18: compute C from C1 using the mapping M 19: return C Given the group size k, LouvainDP starts with a supergraph G1 having ⌊ |V | k ⌋ nodes by randomly permuting the nodeset V and grouping every k consecutive nodes into a supernode (lines 1-4). The permutation prevents the possible bias of node ordering in G. The set of superedges E1 is easily computed from G. Note that m1 = |E1| ≤ m due to the fact that each edge of G appears in one and only one superedge. The domain size is m0 = n 1 (n 1 +1) 2 (i.e. we consider all selfloops in G1). The noisy number of non-zero superedges is m1 = |E1| + Lap(1/ǫ2). Then by one-sided filtering [START_REF] Cormode | Differentially private summaries for sparse data[END_REF], we estimate the threshold θ (line 7) and the number of passing zero superedges s (line 8). For each non-zero superedge, we add a geometric noise and add the superedge to G1 if the noisy value is not smaller than θ. For s zero superedges e1(i, j) / ∈ E1, we draw an integral weight w from the distribution P r[X ≤ x] = 1α x-θ+1 and add e1(i, j) with weight w to G1 if w > 0.

Complexity

LouvainDP runs in O(m) because the loops to compute superedges (Line 5) and to add geometric noises (lines 9-12) cost O(m). We have s = (m0 -m1) Line 7). So the processing of s zero-superedges costs O(m). Moreover, Louvain method (line 17) is empirically linear in m1 [START_REF] Blondel | Fast unfolding of communities in large networks[END_REF]. We come up with the following theorem.

α θ 1+α ≤ m 0 -m 1 1+α (1+α)m 1 m 0 -m 1 = m1 (see
Theorem 4.1. The number of edges in G1 is not larger than 2m. LouvainDP's runtime is O(m).

Privacy Analysis

In LouvainDP, we use a small privacy budget ǫ2 = 0.01 to compute the noisy number of non-zero superedges m1. The remaining privacy budget ǫ1 is used for the geometric mechanism Geom(α). Note that getting a random permutation Vp (line 3) costs no privacy budget. The number of nodes n is public and given the group size k, the number of supernodes n1 is also public. The high-pass filtering technique (Lines 6-16) inherits the privacy guarantee by [START_REF] Cormode | Differentially private summaries for sparse data[END_REF]. By setting ǫ1 = ǫ -ǫ2, LouvainDP satisfies ǫ-differential privacy (see the sequential composition (Theorem 3.3)).

Alternative Input Perturbation Schemes

1K-series [START_REF] Wang | Preserving differential privacy in degree-correlation based graph generation[END_REF], DER [START_REF] Chen | Correlated network data publication via differential privacy[END_REF], TmF [START_REF] Nguyen | Differentially private publication of social graphs at linear cost[END_REF] and EdgeFlip [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF] are the most recent differentially private schemes for graph release that can be classified as input perturbation. While 1K-series and TmF run in linear time, DER and EdgeFlip incur a quadratic complexity. DER and EdgeFlip are therefore tested only on two medium-sized graphs in Section 6.

The expected number of edges by EdgeFlip is

|E G| = (1 -s)m + n(n-1)
4 s (see [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF]) where s = 2 e ǫ +1 is the flipping probability. Substitute s into |E G|, we get

|E G| = m + ( n(n-1) 4 -m) 2 e ǫ +1
. The number of edges in the noisy graph G generated by EdgeFlip increases exponentially as ǫ decreases. To ensure the linear complexity for million-scale graphs, we propose a simple extension of EdgeFlip, called EdgeFlipShrink (see Appendix A.1).

ALGORITHM PERTURBATION

The schemes in the algorithm perturbation category privately sample a node clustering from the input graph without generating noisy sample graphs as in the input perturbation. This can be done via the exponential mechanism. We introduce our main scheme ModDivisive in Section 5.1 followed by HRG-Fixed, a variant of HRG-MCMC [START_REF] Xiao | Differentially private network data release via structural inference[END_REF] runnable on large graphs (Section 5.2). Fig. 2b sketches the basic steps of the algorithm perturbation paradigm.

ModDivisive: Top-down Exploration of Cohesive Groups

Overview

In contrast with the agglomerative approaches (e.g. Louvain method) in which small communities are iteratively merged if doing so increases the modularity, our ModDivisive scheme is a divisive algorithm in which communities at each level are iteratively split into smaller ones. Our goal is to heuristically detect cohesive groups of nodes in a private manner. There are several technical challenges in this process. The first one is how to efficiently find a good split of nodes that induces a high modularity and satisfies ǫ-DP at the same time. The second one is how to merge the small groups to larger ones. We cope with the first challenge by realizing an exponential mechanism via MCMC (Markov Chain Monte-Carlo) sampling with the modularity as the score function (see Theorem 3.2). The second challenge is solved by dynamic programming. We design ModDivisive as a k-ary tree (Fig. 3), i.e. each internal node has no more than k child nodes. The root node (level 0) contains all nodes in V and assigns arbitrarily each node to one of the k groups. Then we run the MCMC over the space of all partitions of V into no more than k nonempty subsets. The resulting subsets are initialized as the child nodes (level 1) of the root. The process is repeated for each child node at level 1 and stops at level maxL. Fig. 3 illustrates the idea with k = 3 for the graph in Fig. 1.

Algorithm

Algorithm 2 sketches the main steps in our scheme ModDivisive. It comprises two phases: differentially private sampling a k-ary tree of depth maxL which uses the privacy budget ǫ1 and finding the best cut across the tree to get a good clustering of nodes which consumes a budget maxL•ǫm.

The first phase (lines 1-14) begins with the creation of eA, the array of privacy budgets allocated to levels of the tree. We use the parameter λ ≥ 1 as the common ratio to form a geometric sequence. The rationale behind the common ratio is to give higher priority to the levels near the root which have larger node sets. By sequential composition (Theorem 3.3), we must have i eA[i] = ǫ1. All internal nodes at level i do the MCMC sampling on disjoint subsets of nodes, so the parallel composition holds (Theorem 3.3). Subsection 5.1.4 analyzes the privacy of ModDivisive in more detail. We use a queue to do a level-by-level exploration. Each dequeued node r's level will be checked. If its level is not larger than maxL, we will run ModMCMC (Algorithm 3) on it (line 9) to get a partition r.part of its nodeset r.S (Fig. 3). Each subset in r.part forms a child node and is pushed to the queue. The second phase (line 15) calls BestCut to find a highly modular partition across the tree.

Differentially Private Nodeset Partitioning Let P be the space of all ways P to partition a nodeset A to no more than k disjoint subsets, the direct application of exponential mechanism needs the enumeration of P. The probability of a partition P being sampled is exp( ǫp 2∆Q Q(P, G)) The space P is connected. It is straightforward to verify that the transitions performed in line 3 of ModMCMC are reversible and ergodic (i.e. any pair of nodeset partitions can be connected by a sequence of such transitions). Hence, ModMCMC has a unique stationary distribution in equilibrium. By empirical evaluation, we observe that ModMCMC converges after K|r.S| steps for K = 50 (see Section 6.3).

P ′ ∈P exp( ǫp 2∆Q Q(P ′ , G)) (4
Each node r in the tree is of type NodeSet. This structure consists of an array r.part where r.part[u] ∈ {0..k -1} is the group id of u. To make sure that ModMCMC runs in linear time, we must have a constant time computation of modularity Q(P ) (line 4 of Algorithm 3). This can be done with two helper arrays: the number of intra-edges r.lc[0..k -1] and the total degree of nodes r.dc[0..k -1] in each group. The modularity Q is computed in O(k) (Formula 1) using r.lc, r.dc. When moving node u from group i to group j, r.lc and r.dc are updated accordingly by checking the neighbors of u in G. The average degree d is a constant, so the complexity of ModMCMC is linear in the number of MCMC steps.

Algorithm 3 ModMCMC

Input: graph G, nodeset r.S, group size k, privacy budget ǫp Output: sampled partition r.part 1: initialize r.part with a random partition P 0 of k groups 2: for each step i in the Markov chain do 3:

pick a neighboring partition P ′ of P i-1 by randomly selecting node u ∈ r.S and moving u to another group.

4:

accept the transition and set P i = P ′ with probability

min(1, exp( ǫp 2∆Q Q(P ′ ,G)) exp( ǫp 2∆Q Q(P i-1 ,G))
) 5: // until equilibrium is reached 6: return a sampled partition r.part = P i Finding the Best Cut Given the output k-ary tree R with the root node root, our next step is to find the best cut across the tree. A cut C is a set of nodes in R that cover all nodes in V . As an example, a cut C in Fig. 3 returns the clustering [{0,1,2,3}, {4}, {5,6,11}, {8,9,10,7,12}]. Any cut has a modularity score. Our goal is to find the best cut, i.e. the cut with highest modularity, in a private manner.

We solve this problem by a dynamic programming technique. Remember that modularity is an additive quantity (c.f. Formula 1). By denoting opt(r) as the optimal modularity for the subtree rooted at node r, the optimal value is opt(root). The recurrence relation is straightforward

opt(r) = max{Q(r), t∈r.children opt(t)}
This idea is realized in three steps. The pseudo-code is provided in Appendix A.2. The first step uses a queue to fill a stack S. The stack ensures any internal node to be considered after its child nodes. The second step solves the recurrence relation. Because all modularity values are sensitive, we add Laplace noise Laplace(∆Q/ǫm). The global sensitivity ∆Q = 3/m (see Theorem 5.2), so we need only a small privacy budget for each level (ǫm = 0.01 is enough in our experiments). The noisy modularity modn is used to decide whether the optimal modularity at node r is by itself or by the sum over its children. The final step backtracks the best cut from the root node.

Complexity

ModDivisive creates a k-ary tree of height maxL. At each node r of the tree other than the leaf nodes, ModMCMC is run once. The run time of ModMCMC is O(K •|r.S|) thanks to the constant time for updating the modularity (line 4 of ModMCMC). Because the union of nodesets at one level is V , the total runtime is O(K • |V | • maxL). BestCut only incurs a sublinear runtime because the size of tree is always much smaller than |V |. The following theorem states this result Theorem 5.1. The time complexity of ModDivisive is linear in the number of nodes n, the maximum level maxL and the burn-in factor K.

Privacy Analysis

We show that ModMCMC satisfies differential privacy. The goal of MCMC is to draw a random sample from the desired distribution. Similarly, exponential mechanism is also a method to sample an output x ∈ X from the target distribution with probability proportional to exp( ǫu(x) 2∆u ) where u(x) is the score function (x with higher score has bigger chance to be sampled) and ∆u is its sensitivity. The idea of using MCMC to realize exponential mechanism is first proposed in [START_REF] Chaudhuri | A near-optimal algorithm for differentially-private principal components[END_REF] and applied to ǫ-DP graph release in [START_REF] Xiao | Differentially private network data release via structural inference[END_REF].

In our ModMCMC, the modularity Q(P, G) is used directly as the score function. We need to quantify the global sensitivity of Q. From Section 3.1, we have the following definition Definition 5.1. (Global Sensitivity ∆Q)

∆Q = max P,G,G ′ |Q(P, G) -Q(P, G ′ )| (5)
We prove that ∆Q = O(1/m) in the following theorem whose proof is provided in Appendix A.3.

A Variant of HRG-MCMC

Similar to DER and EdgeFlip, HRG-MCMC [START_REF] Xiao | Differentially private network data release via structural inference[END_REF] runs in quadratic time due to its costly MCMC steps. Each MCMC step in HRG-MCMC takes O(n) to update the tree. To make it runnable on million-scale graphs, we describe briefly a variant of HRG-MCMC called HRG-Fixed (see Fig. 4). Instead of exploring the whole space of HRG trees as in HRG-MCMC, HRG-Fixed selects a fixed binary tree beforehand. We choose a balanced tree in our HRG-Fixed. Then HRG-Fixed realizes the exponential mechanism by sampling a permutation of n leaf nodes (note that each leaf node represents a graph node). The next permutation is constructed from the current one by randomly choosing a pair of nodes and swap them. The bottom of Fig. 4 illustrates the swap of two nodes a and d. The log-likelihood L still has the sensitivity ∆L = 2 ln n as in HRG-MCMC [START_REF] Xiao | Differentially private network data release via structural inference[END_REF]. Each sampling (MCMC) operation is designed to run in O( d•log n). By running K • n MCMC steps, the total runtime of HRG-Fixed is O(K • m • log n), feasible on large graphs. The burn-in factor K of HRG-Fixed is set to 1,000 in our evaluation.

EXPERIMENTS AND RESULTS

In this section, our evaluation aims to compare the performance of the competitors by clustering quality and efficiency. The clustering quality is measured by the modularity Q and the average F1-score in which the modularity is the most important metric as we aim at highly modular clusterings. The efficiency is measured by the running time. All algorithms are implemented in Java and run on a desktop PC with Intel Core i7-4770@ 3.4Ghz, 16GB memory. Two medium-sized and two large real graphs are used in our experiments 1 . as20graph is the graph of routers comprising the internet. ca-AstroPh is a co-authorship network where two authors are connected if they publish at least one paper together. amazon is a product co-purchasing network where the graph contains an undirected edge from i to j if a product i is frequently co-purchased with product j. youtube is a video-sharing web site that includes a social network. Table 2 shows the characteristics of the graphs. The rows Com(munities) and Mod(ularity) are the output of Louvain method. The row 0.5ln n is the largest privacy budget ǫ considered in the experiments. The number of samples in each test case is 20.

The schemes are abbreviated as 1K-series (1K), EdgeFlip (EF), Top-m-Filter (TmF), DER, LouvainDP (LDP), Mod-Divisive (MD), HRG-MCMC and HRG-Fixed.

Quality Metrics

Apart from modularity Q, we use F1, the average F1-score, following the benchmarks in [START_REF] Yang | Defining and evaluating network communities based on ground-truth[END_REF]. The F1 score of a set A with respect to a set B is defined as the harmonic mean H of the precision and the recall of A against B. We define

prec(A, B) = |A∩B| |A| , recall(A, B) = |A∩B| |B| F1(A, B) = 2.prec(A, B).recall(A, B) prec(A, B) + recall(A, B)
Then the average F1 score of two sets of communities C and C ′ is defined as

F1(A, C) = max i F1(A, ci), ci ∈ C = {c1, .., cn} F1(C, C ′ ) = 1 2|C| c i ∈C F1(ci, C ′ ) + 1 2|C ′ | c i ∈C ′ F1(ci, C)
We choose the output clustering of Louvain method as the ground truth for two reasons. First, the evaluation on the real ground truth is already done in [START_REF] Prat-Pérez | High quality, scalable and parallel community detection for large real graphs[END_REF] and Louvain method is proven to provide high quality communities. Second, the real ground truth is a set of overlap communities whereas the schemes in this paper output only non-overlap communities. The chosen values of ǫ are {0.1ln n, 0.2ln n, 0.3ln n, 0.4ln n, 0.5ln n}.

Note that there is a well-known resolution limit of modularity [10, Section VI.C]. This phenomenon may occur in top-down community detection approaches in which the modularity of a clustering of large communities is not much different from that of a clustering of many smaller communities. Similarly, the F1-score may not be consistent with modularity. For example, on amazon graph, non-private Louvain method returns a clustering of 52,655 small communities with modularity 0.678 after first round. In the final round, it returns a clustering of 257 large communities with modularity 0.926. The F1-score between 52,655 small communities and 257 large communities (as ground truth) is 0.01. On youtube graph, the first round and the last round of non-private Louvain return 147,361 and 13,485 communities with modularity 0.644 and 0.710 respectively. The F1-score between the two clusterings is 0.14. These results confirm the inherent limitations of modularity and F1-score for community detection.

Performance of LouvainDP

We test LouvainDP for the group size k ∈ {4, 8, 16, 32, 64}. The results on youtube are displayed in Fig. 5. We observe a clear separation of two groups k = 4, 8 and k = 16, 32, 64. As k increases, the modularity increases faster but also saturates sooner. Similar separations appear in avg.F1Score and the number of communities. Non-trivial modularity scores in several settings of (k, ǫ) indicate that randomly grouping of nodes and high-pass filtering superedges do not destroy all community information of the original graph. At ǫ = 0.5 ln n and k = 4, 8, the total edge weight in G1 is very low ( < 0.05m), so many supernodes of G1 are disconnected and Louvain method outputs a large number of communities. The reason is that the threshold θ in LouvainDP is an integral value, so causing abnormal leaps in the total edge weight of G1. We pick k = 8, 64 for the comparative evaluation (Section 6.4).

Performance of ModDivisive

The effectiveness of ModDivisive is illustrated in Fig. 6 for graph youtube and λ = 2.0, K = 50. We select six pairs of (k, maxL) by the set {(2,10), [START_REF] Chaudhuri | A near-optimal algorithm for differentially-private principal components[END_REF][START_REF] Cormode | Differentially private summaries for sparse data[END_REF], [START_REF] Chen | Correlated network data publication via differential privacy[END_REF][START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF], [START_REF] Clauset | Hierarchical structure and the prediction of missing links in networks[END_REF][START_REF] Chen | Correlated network data publication via differential privacy[END_REF], [START_REF] Clauset | Finding community structure in very large networks[END_REF][START_REF] Chen | Correlated network data publication via differential privacy[END_REF], [START_REF] Fortunato | Community detection in graphs[END_REF][START_REF] Chaudhuri | A near-optimal algorithm for differentially-private principal components[END_REF]}. Modularity increases steadily with ǫ while it is not always the case for avg.F1Score.

We choose λ = 2.0 to obtain a good allocation of ǫ among the levels. Fig. 7a shows the modularity for different values of λ. Note that λ = 1.0 means ǫ is equally allocated to the maxL levels. By building a k-ary tree, we reduce considerably the size of the state space P for MCMC. As a result, we need only a small burn-in factor K. Looking at Fig. 7b, we see that larger K = 100 results in only tiny increase of modularity in comparison with that of K = 50.

Comparative Evaluation

In this section, we report a comparative evaluation of Lou-vainDP and ModDivisive against the competitors in Fig. 9. The dashed lines in subfigures 9i, 9j, 9k and 9l represent the ground-truth number of communities by Louvain method. ModDivisive performs best in most of the cases.

On as20graph and ca-AstroPh, HRG-MCMC outputs the whole nodeset V with the zero modularity while 1K-series, TmF, DER also give useless clusterings. EdgeFlip produces good quality metrics exclusively on ca-AstroPh while Lou-vainDP returns the highest modularity scores on as20graph. However, the inherent quadratic complexity of EdgeFlip On the two large graphs, ModDivisive dominates the other schemes by a large margin in modularity and avg.F1Score. LouvainDP is the second best in modularity at k = 64. Our proposed HRG-Fixed is consistent with ǫ and has good performance on youtube. Note that HRG-MCMC is infeasible on the two large graphs due to its quadratic complexity. Again, 1K-series, TmF and EdgeFlipShrink provide the worst quality scores with the exception of 1K-series's avg.F1Score on youtube.

The runtime of the linear schemes is reported in Fig. 8. EdgeFlipShrink, 1K-series, TmF and LouvainDP benefit greatly by running Louvain method on the noisy output graph G. ModDivisive and HRG-Fixed also finish their work quickly in O(K • n • maxL) and O(K • m • log n) respectively.

CONCLUSION

We have given a big picture of the problem ǫ-DP community detection within the two categories: input and algorithm perturbation. We analyzed the major challenges of community detection under differential privacy. We explained why techniques borrowed from k-Means fail and how the difficulty of ǫ-DP recommender systems enables a relaxation of privacy budget. We proposed LouvainDP and ModDivisive as the representatives of input and algorithm perturbations respectively. By conducting a comprehensive evaluation, we revealed the advantages of our methods. ModDivisive steadily gives the best modularity and avg.F1Score on large graphs while LouvainDP outperforms the remaining input perturbation competitors in certain settings. HRG-MCMC/HRG-Fixed give low modularity clusterings, indicating the limitation of the HRG model in divisive CD. The input perturbation schemes DER, EF, 1K-series and TmF hardly deliver any good node clustering except EF on the two medium-sized graphs.

For future work, we plan to develop an ǫ-DP agglomerative scheme based on Louvain method and extend our work for directed graphs and overlapping community detection under differential privacy. Instead of outputting G, EdgeFlipShrink computes Ĝ that has the expected number of edges m by shrinking E G. First, the algorithm compute the private number of edges m using a small budget ǫ2 (Lines 2-3 ). The new flipping probability s is updated (Line 5). The noisy expected number of edges in the original EdgeFlip is shown in Line 6. We obtain the shrinking factor p = m m 0 (Line 7). Using p, every 1-edge is sampled with probability 1-s 2 p instead of 1-s 2 as in [START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF]. The remaining 0-edges are randomly picked from EG as long as they do not exist in Ĝ (Lines [START_REF] Kenthapadi | Privacy via the johnson-lindenstrauss transform[END_REF][START_REF] Mcsherry | Mechanism design via differential privacy[END_REF][START_REF] Mcsherry | Privacy integrated queries: an extensible platform for privacy-preserving data analysis[END_REF][START_REF] Medus | Detection of community structures in networks via global optimization[END_REF][START_REF] Mülle | Privacy-integrated graph clustering through differential privacy[END_REF][START_REF] Newman | Finding and evaluating community structure in networks[END_REF].

The if Ĝ does not contain (i, j) then

18:

add edge (i, j) to Ĝ

19:

n 0 --

20: return Ĝ′

A.2 BestCut Algorithm

Algorithm 5 outlines the steps of BestCut (Section 5.1.2).

A.3 Proof of Theorem 5.2

Proof. Given the graph G and a partition P of a nodeset Vp ⊆ V (for any node of the k-ary tree other than the root node, its nodeset Vp is a strict subset of V ), the neighboring graph G ′ has E G ′ = EG ∪ e. We have two cases Case 1. The new edge e is an intra-edge within the community s. The modularity Q(P, G) is Because ∆Q is the absolute value of d1, we consider the most positive and the most negative values of d1. Remember that k c dc ≤ 2m, so the positive bound d1 <

4:

S.push(r)

5:

for child node r i in r.children do

6:

Q.enqueue(r i ) 7: dictionary sol ← ∅ 8: while S is not empty do 9:

Figure 1 :

 1 Figure 1: Louvain method Theorem 3.2. (Exponential mechanism [15]) Given a score function u : (G×O) → R for a graph G, the mechanism A that samples an output O with probability proportional to exp( ǫ.u(G,O) 2∆u ) satisfies ǫ-differential privacy. 2

Figure 2 :

 2 Figure 2: Two categories of ǫ-DP community detection Formula 1) is non-monotone, i.e. for two disjoint sets of nodes A and B, Q(A ∪ B) may be larger, smaller than or equal toQ(A) + Q(B).To further emphasize the difficulty of ǫ-DP community detection, we found that IDC scheme[START_REF] Gupta | Iterative constructions and private data release[END_REF] using Sparse Vector Technique [9, Section 3.6] is hardly feasible. As shown in Algorithm 1 of[START_REF] Gupta | Iterative constructions and private data release[END_REF], to publish a noisy graph that can approximately answer all cut queries with bounded error m 0.25 n/ǫ 0.5 , IDC must have B(α) "yes" queries among all k queries. B(α) may be as low as √ m but k = 2 2n . In the

Figure 3 :

 3 Figure 3: Example of ModDivisive with k = 3. A cut C is shown by the dot-dashed line

Theorem 5 . 2 .mFigure 4 :

 524 Figure 4: HRG-Fixed

Figure 5 :

 5 Figure 5: LouvainDP on youtube

Figure 6 :

 6 Figure 6: ModDivisive on youtube with λ = 2.0, K = 50

  K (at ǫ = 0.5ln n)

Figure 7 :

 7 Figure 7: ModDivisive: modularity vs. λ and K on amazon

Figure 8 :

 8 Figure 8: Runtime (in second)

Figure 9 :

 9 Figure 9: Quality metrics and the number of communities on four graphs

2 s - 1 ) -ǫ 2 5: s = 2 e ǫ +1 6 : 4 s 7 : p = m m 0 8 : 13 :

 2122647813 expected edges of Ĝ is E[ Ĝ] = E[ m] = m and the running time of EdgeFlipShrink is O(m). Algorithm 4 EdgeFlipShrink(G, s) Input: undirected graph G, flipping probability s Output: anonymized graph Ĝ 1: Ĝ ← ∅ 2: ǫ 2 = 0.01 3: m = m + Lap(1/ǫ 2 )4: ǫ = ln( m 0 = (1 -s) m + n(n-1) // process 1-edges 9: n 1 = 0 10: for edge (i, j) ∈ E G do 11:add edge (i, j) to Ĝ with prob. // process 0-edges 14: n 0 = m -n 1 15: while n 0 > 0 do 16: random pick an edge (i, j) / ∈ E G 17:

  ) 2 ). The difference d1 = Q(P, G ′ ) -Q(P, G)

r 3 m+1 < 3 m . Case 2 .

 32 ← S.pop(), r.modn = r.mod + Laplace(∆Q/ǫm) 10: if r is a leaf node then 11: sol.put(r.id, (val = r.modn, self = T rue)) sm = r i ∈r.children sol[r i .id].modn 14: if r.modn < sm then 15:sol.put(r.id, (val = sm, self = F alse)) sol.put(r.id, (val = r.modn, self = T rue))18: list C ← ∅, queue Q ← root 19: while Q is not empty do 20: r ← Q.dequeue()21:if sol[r.id].self == T rue then 22:C = C ∪ {r} 23: else 24:for child node r i in r.children do 25:Q.enqueue(r i ) return C use the constraints k c lc ≤ m and ds ≤ 2m, so d1 > 1 m+1 -m m(m+1) -2m+1 (m+1) 2 > -2 m+1 .As a result, ∆Q = |d1| < The new edge e is an inter-edge between the communities s and t. Similarly, we have Q(P, G) = k c ( lc m -) 2 . Again, we consider the most positive and the most negative values of d2, using the constraint k c dc ≤ 2m, the positive bound d2 < (2m+1)4m 2 4m 2 (m+1) 2 < 2 m+1 . For the negative bound, we use the constraintsk c lc ≤ m and ds + dt ≤ 2m, so d2 > -m m(m+1) -4m+2 4(m+1) 2> -2 m+1 . As a result, ∆Q = |d2| < 2 m+1 < 2 m . To recap, in both cases ∆Q < 3 m .

Table 1 :

 1 List of notations

)

  Algorithm 2 ModDivisive Input: graph G, group size k, privacy budget ǫ, max level maxL, ratio λ, BestCut privacy at each level ǫm Output: noisy partition C 1: compute the array eA[0..maxL -1] s.t. = ǫ 1 , eA[i] = eA[i + 1] * λ where ǫ 1 = ǫ -maxL.ǫm |A|, i) where S(|A|, i) is the Stirling number of the second kind, S(n, k) ≈ k n k! . This sum is exponential in |A|, so enumerating P is computationally infeasible. Fortunately, MCMC can help us simulate the exponential mechanism by a sequence of local transitions in P.

	2: initialize the root node with nodeset V 3: root = NodeSet(G, V, k)
	4: root.level = 0
	5: queue Q ← root 6: while Q is not empty do
	7:	r ← Q.dequeue()
	8:	if r.level < maxL then
	9: 10: 11: 12: 13: 14:	r.part = ModMCMC(G, r.S, k, eA[r.level]) for subset S i in r.part do P i = NodeSet(G, S i , k) P i .level = r.level + 1 r.child i ← P i Q.enqueue(P i )
	15: C ← BestCut(root, ǫm)
	16: return	C
	However, |P| =	k i=1 S(

i eA

[i]

Table 2 :

 2 Characteristics of the test graphs

		as20graph ca-AstroPh amazon	youtube
	Nodes	6,474	17,903 334,863 1,134,890
	Edges	12,572	196,972 925,872 2,987,624
	Com	30	37	257	13,485
	Mod	0.623	0.624	0.926	0.710
	0.5ln n	4.4	4.9	6.4	7.0

1

  m+1 + (2m+1)4m 2 4m 2 (m+1) 2 < 3 m+1 . For the negative bound, we Algorithm 5 BestCut Input: undirected graph G, root node root, privacy budget at each level ǫm Output: best cut C 1: stack S ← ∅, queue Q ← root 2: while Q is not empty do 3:

r ← Q.dequeue()

http://snap.stanford.edu/data/index.html

APPENDIX A. ALGORITHMS AND PROOFS

A.1 EdgeFlipShrink