Jean B Lasserre

A MODIFIED CHRISTOFFEL FUNCTION AND ITS ASYMPTOTIC PROPERTIES

We introduce a certain variant (or regularization) Λµ n of the standard Christoffel function Λ µ n associated with a measure µ on a compact set Ω ⊂ R d . Its reciprocal is now a sum-of-squares polynomial in the variables (x, ε), ε > 0. It shares the same dichotomy property of the standard Christoffel function, that is, the growth with n of its inverse is at most polynomial inside and exponential outside the support of the measure. Its distinguishing and crucial feature states that for fixed ε > 0, and under weak assumptions, limn→∞ ε -d Λµ n (ξ, ε) = f (ζ ε ) where f (assumed to be continuous) is the unknown density of µ w.r.t. Lebesgue measure on Ω, and ζ ε ∈ B∞(ξ, ε) (and so f (ζ ε ) ≈ f (ξ) when ε > 0 is small). This is in contrast with the standard Christoffel function where if limn→∞ n d Λ µ n (ξ) exists, it is of the form f (ξ)/ωE(ξ) where ωE is the density of the equilibrium measure of Ω, usually unknown. At last but not least, the additional computational burden (when compared to computing Λ µ n ) is just integrating symbolically the monomial basis (x α ) α∈N d n on the box {x : xξ ∞ < ε/2}, so that 1/ Λµ n is obtained as an explicit polynomial of (ξ, ε).

Introduction

The Christoffel function Λ µ n : R d → R + and the Christoffel-Darboux kernel are both associated with a measure µ whose support Ω ⊂ R d is here a compact set, and are indexed by a "degree" n. They originate in the theory of approximation and orthogonal polynomials. Among its interesting features:

• (1) The Christoffel function identifies the support of µ by a dichotomy of its behavior as n grows, depending on whether the point is inside or outside the support of the measure. Namely, inside the support the growth of its reciprocal is at most polynomial in the degree n, and at least exponential outside the support.

• (2) Under some regularity assumptions, the Christoffel function is also related to the density f of µ w.r.t. Lebesgue measure on Ω. Indeed then d+n d Λ µ n → f /ω E pointwise in Ω (and uniformly on compact subsets of Ω)

where ω E is a so-called equilibrium measure of Ω. However in general the equilibrium measure is unknown, except in a few special cases of sets Ω with special geometry; see e.g. [1, §9.7 ], [START_REF] Kroó | Christofffel functions and universality on the boundary of the ball[END_REF][START_REF] Kroó | Christofffel functions and universality in the bulk for multivariate polynomials[END_REF], [START_REF] Xu | Asymptotics of the Christoffel functions on a simplex in R d[END_REF].

• (3) Its explicit expression as a rational function, is easy to compute, e.g. via the inverse of the moment matrix of µ (provided that its size is reasonable and compatible with state-of-the-art linear algebra softwares).

As advocated in [START_REF] Lasserre | The Christoffel-darboux Kernel for Data Analysis[END_REF], in addition to being a mathematical object interesting in its own, the Christoffel function (CF) turns out to also provide an efficient and easy-to-use tool to help solve some problems in data analysis, e.g. outlier detection, support inference, and density estimation. In such problems, the underlying measure of interest µ is usually not available and one has rather access to a finite sample of data points generated according to µ. Then in lieu of the unknown µ, one uses the empirical measure associated with the cloud of finitely many data points. Remarkably, and even though the geometry of the support of the empirical measure is quite trivial, its associated CF is still close (in a certain sense) to that of µ, and inherits interesting features of µ that can be exploited in data analysis; see e.g. [START_REF] Lasserre | Sorting out typicality via the inverse moment matrix SOS polynomial[END_REF][START_REF] Lasserre | On the Christoffel function and classification in data analysis[END_REF][START_REF] Lasserre | The Christoffel-darboux Kernel for Data Analysis[END_REF][START_REF] Lasserre | The empirical Christoffel function with applications in data analysis[END_REF].

Moreover, some recent works have revealed additional properties of the CF, as well as connections (in the author's opinion some even surprising) with seemingless disconnected topics, e.g., convex duality, certificates of positivity in real algebraic geometry, Pell's equation, equilibrium measure of a compact set; see [START_REF] Lasserre | A disintegration of the Christoffel function[END_REF][START_REF] Lasserre | Pell's equation, sum-of-squares, and equilibrium measures of compacts sets[END_REF].

For more details on the Christoffel-Darboux kernel and the Christoffel function, the interested reader is referred to e.g. [START_REF] Dunkl | Orthogonal Polynomials of Several Variables[END_REF][START_REF] Lasserre | The Christoffel-darboux Kernel for Data Analysis[END_REF][START_REF] Nevai | Géza Freud, orthogonal polynomials and Christoffel functions. A case study[END_REF][START_REF] Simon | The Christoffel Darboux kernel[END_REF][START_REF] Stahl | General Orthogonal Polynomials, Encyclopedia of Mathematics and its Applications[END_REF][START_REF] Xu | Asymptotics of the Christoffel functions on a simplex in R d[END_REF] and many references therein.

Contribution. Our motivation is essentially concerned with point (2) above. Indeed our main goal is to still recover the density f of µ asymptotically, but without the annoying factor ω E of the equilibrium measure because it is unknown in general. But simultaneously, we also want to maintain the support inference capability in (1) as well as an efficient computation (point (3) above).

To achieve these goals, we introduce the following variant (or regularization) of the Christoffel function, namely the function Λµ n : R d × R → R + defined by:

(1.1) Λµ n (ξ, ε) := { inf p∈R[x]n { p 2 dµ : B∞(ξ,ε) p(x) dx ε d = 1 } , for all (ξ, ε) ∈ R d × R + and n ∈ N, where B ∞ (ξ, ε) := {x : x -ξ ∞ < ε/2} (with Lebesgue volume ε d ).
As we will see, 1/ Λµ n is obtained as an explicit polynomial of (ξ, ε), and can be also seen as a a polynomial of ξ, parametrized by ε > 0 fixed. In particular, we are interested in the asymptotic behaviors of Λµ n (ξ, ε) with ε > 0 fixed, as well as Λµ n (ξ, 1/n), as n grows.

More precisely, let Ω ⊂ R d be compact with nonempty interior, and let µ be a Borel measure on Ω with density f with respect to (w.r.t.) Lebesgue measure on Ω. In particular, for every n ∈ N, the moment matrix M n (µ) is non singular. We prove that Λµ n in (1.1) has the following properties:

• 1/ Λµ
n is a sum of squares polynomial with explicit form:

(1.2) (ξ, ε) → Λµ n (ξ, ε) -1 = ṽn (ξ, ε) T M n (µ) -1 ṽn (ξ, ε)
, where, with v n (x) being the vector of the monomial basis (

x α ) α∈N d n , ṽn (ξ, ε) := B∞(ξ,ε) v n (x) dx ε d .
Importantly, each entry ṽn,α , α ∈ N d n , of the vector ṽn , is obtained explicitly as a polynomial in (ξ, ε).

• With ε > 0 fixed, and assuming that f is bounded and f > 0 on Ω,

(1.3) lim n→∞ ε d Λµ n (ξ, ε) -1 = B∞(ξ,ε) 1 f dx ε d = ε d 1 B∞(ξ,ε) ε d f 2 L 2 (µ) , for all ξ with B ∞ (ξ, ε) ⊂ Ω. Moreover, if f is continuous then (1.4) lim n→∞ ε -d Λµ n (ξ, ε) = f (ζ ε ) , for some ζ ε ∈ B ∞ (ξ, ε), and so f (ζ ε ) ≈ f (ξ) if ε is small enough.
• Finally, and as for the standard CF Λ µ n , the growth with n of Λµ n (ξ, ε) -1 (ε > 0 fixed) and Λµ n (ξ, 1/n) -1 , is at most polynomial if ξ ∈ int(Ω) and at least exponential if ξ ∈ Ω.

So the extended CF Λµ

n essentially shares same important features as the standard CF with the advantage that the limit of ε -d Λµ n (ξ, ε) (ε > 0 fixed) as n increases, is close to f (ξ) where f is the density of µ, a highly desirable feature. Moreover, the additional price to pay to obtain its explicit form (1.2) as a function of (ξ, ε), is rather negligible as it only requires the closed form expression of B∞(ξ,ε) x α dx, for all α ∈ N d n . Interpretation. For fixed ε > 0, the function ξ → Λµ n (ξ, ε) has a simple interpretation. Let f be the density of µ with respect to Lebesgue measure on Ω and assume that f ≥ γ > 0 on Ω. Let ε ξ be the linear functional h → ε ξ (h) := h dλ where λ is the uniform probability mesure on B ∞ (ξ, ε).

If ξ ∈ int(Ω) is such that B ∞ (ξ, ε) ⊂ Ω, then ε ξ ∈ L 2 (µ) and is the function x → 1 B∞(ξ,ε) (x)/τ ε f ∈ L 2 (µ) where τ ε = vol(B ∞ (ξ, ε)). Next, let ε ξ,n ∈ R[x] n be the orthogonal projection of ε ξ on the finite-dimensional subspace R[x] n of L 2 (µ). Then Λµ n (ξ, ε) = 1/ ε ξ,n 2 L 2 (µ) with limit (1.3) as n → ∞.
Λµ n (ξ, ε) can be seen as an ε-regularization of Λ µ n (ξ), where instead of working with the point evaluation at ξ (i.e., the Dirac measure δ ξ which is not in L 2 (µ)), one rather works with the element ε ξ of some neighborhood N (δ ξ , η) of δ ξ (in the weak-topology of the space of finite Borel signed measures) because it can also be considered as an element of L 2 (µ). In doing so with ε > 0 fixed, lim n→∞ ε -d Λµ n (ξ, ε) exists (which is not the case for Λ µ n (ξ)) and is close to f (ξ) if ε is small and f is continuous, a highly desirable feature.

|α| = i α i ) with cardinal s(n) = d+n d . Let v n (x) = (x α
) α∈N d n be the vector of monomials up to degree n, and let Σ[x] n ⊂ R[x] 2n be the convex cone of polynomials of total degree at most 2n which are sum-of-squares (in short SOS). For every p

∈ R[x] n write x → p(x) = p, v n (x) , ∀x ∈ R d ,
where p ∈ R s(n) is the vector of coefficients of p in the monomial basis (x α ) α∈N d . For a real symmetric matrix A = A T the notation A 0 (resp. A 0) stands for A is positive semidefinite (p.s.d.) (resp. positive definite (p.d.)).

The support of a Borel measure µ on R d is the smallest closed set A such that µ(R d \ A) = 0, and such a set A is unique. With S ⊂ R d compact, denote by C (S) the Banach space of real continuous functions on S equipped with the sup-norm. Its topological dual C (S) * is the Banach space M (S) of finite signed Borel measures on S, equipped with the total-variation norm.

Moment matrix. Let µ be a finite Borel measure on R d with all moments µ = (µ α ) α∈N d assumed to be finite. The (degree-n) moment matrix M n (µ) associated with µ is the real symmetric matrix with rows and columns indexed by N d n (hence of size s(n)), and with entries

M n (µ)(α, β) := x α+β dµ = µ α+β , α, β ∈ N d n .
Obviously, M n (µ) 0 for all n since

p, M n (µ) p = p 2 dµ ≥ 0 , ∀p ∈ R[x] n .
Christoffel function. We here assume that M n (µ) 0 for all n ∈ N, and therefore the inverse M n (µ) -1 is well-defined for all n ∈ N. In particular this is true in our case of interest, i.e., when the support Ω ⊂ R d of µ is compact with nonempty interior and µ has a density w.r.t. Lebesgue measure on Ω.

The (degree-n) Christoffel function Λ µ n : R → R + , associated with µ, is defined by:

x → Λ µ n (x) -1 := v n (x) T M n (µ) -1 v n (x) , ∀x ∈ R d . Alternatively, if (P α ) α∈N d ⊂ R[x]
is a family of polynomials which are orthonormal with respect to µ, then

(2.1) Λ µ n (x) -1 = α∈N d n P α (x) 2 , ∀x ∈ R d .
The Christoffel function has also a variational formulation. Namely:

(2.2) Λ µ n (ξ) = inf p∈R[x]n { p 2 dµ : p(ξ) = 1 } , ∀ξ ∈ R d . Problem (2.
2) is a quadratic convex optimization problem that can be solved efficiently. Its unique optimal solution p * ∈ R[x] n reads:

x → p * (x) := α∈N d n P α (ξ)P α (x) α∈N d n P α (ξ) 2 = K µ n (ξ, x) K µ n (ξ, ξ) , x ∈ R d ,
where

K µ n : R d × R d → R, defined by: (2.3) (x, y) → K µ n (x, y) := α∈N n t P α (x) P α (y) , x, y ∈ R d ,
is the Christoffel-Darboux kernel associated with µ. In particular

Λ µ n (ξ) -1 = K µ n (ξ, ξ) , ∀ξ ∈ R d .

2.2.

A regularization and parametrization of the Christoffel function. Let µ be a finite Borel probability measure on a compact set Ω ⊂ R d , with density f w.r.t. Lebesgue measure on Ω, i.e., dµ = 1 Ω (x)f (x)dx and

Ω f (x) dx = 1. Let L 2 (Ω, µ) (in short, L 2 (µ)
) be the usual Hilbert space of square integrable functions w.r.t. µ. Next, given ε > 0 and ξ ∈ R d , let

B ∞ (ξ; ε) := {x ∈ R d : x -ξ ∞ ≤ ε/2} (2.4) τ ε := vol(B ∞ (ξ, ε)) = ε d (2.5) dφ ε ξ = 1 B∞(ξ,ε) (x) dx τ ε , ε > 0 ; φ 0 ξ := δ {ξ} . (2.6) Proposition 2.1. With φ ε ξ as (2.6), φ ε ξ ⇒ φ 0 ξ = δ {ξ} as ε ↓ 0, i.e., lim ε→0 h dφ ε ξ = h(ξ) = h dφ 0 ξ , ∀h ∈ C (R n ) .
Proof. Observe that as h is continuous,

h dφ ε ξ = h(ζ ε ) , for some ζ ε ∈ B ∞ (ξ, ε),
and therefore h(ζ ε ) → h(ξ) as ε ↓ 0.
In Proposition 2.4, the notation φ ε ξ ⇒ φ 0 ξ is standard and stands for the weak convergence of probability measures, i.e., for the weak-topology

σ(M (Ω), C (Ω)) of M (Ω). With ε ≥ 0, the extended Christoffel function Λµ n : R d × R + → R + is defined by: (2.7) (ξ, ε) → Λµ n (ξ, ε) := min p∈R[x]n { Ω p 2 dµ : p dφ ε ξ = 1 } , ξ ∈ R d ,
where we have included the case ε = 0 with φ 0 ξ = δ {ξ} . Notice that the standard Christoffel function Λ µ n satisfies

Λ µ n (ξ) = inf p∈R[x]n { p 2 dµ : p(ξ) = 1 } = inf p∈R[x]n { p 2 dµ : p dφ 0 ξ = 1 } = Λµ n (ξ, 0) , ξ ∈ R d , that is, Λµ n (ξ, 0) = Λ µ n (ξ), for all ξ ∈ R d .
An explicit form. Importantly, the extended Christoffel function Λµ n is obtained as an explicit rational function of ξ and ε. More precisely, 1/ Λµ n is obtained as an explicit sum-of-squares (SOS) polynomial of (ξ, ε). Lemma 2.2. Let Λµ n be as in (2.7) and let (P α ) α∈N d t be a family of polynomials orthonormal w.r.t. µ. Then the unique optimal solution p * n ∈ R[x] n of (2.7) satisfies:

x → p * n (x) = Λµ n (ξ, ε) v n (x) T M n (µ) -1 v n (y) dφ ε ξ (y) (2.8) = Λµ n (ξ, ε) K µ n (x, y) dφ ε ξ (y) . (2.9)
In addition,

(2.10) Λµ n (ξ, ε) -1 = ṽn (ξ, ε) T M n (µ) -1 ṽn (ξ, ε) , ξ ∈ R n , where ṽn ∈ R[x, ε] n is defined by: (2.11) ξ → ṽn (ξ, ε) := 1 τ ε B∞(ξ,ε) v n (y) dy , n ∈ N . Proof. Rewrite (2.7) Λµ n (ξ, ε) = min p∈R s(n) { p, M n (µ) p : p , v n (x) dφ ε ξ = 1 } , ξ ∈ R d ,
which is a convex quadratic optimization problem. Its unique optimal solution

p * n ∈ R s(n) satisfies 2 M n (µ) p * n = λ * v n (y) dφ ε ξ (y) ,
for some scalar λ * . Hence λ * = 2 Λµ n (ξ, ε), and

p * n = λ * 2 M n (µ) -1 v n (y) dφ ε ξ (y) ,
which in turn yields

x → p * n (x) = p * n , v n (x) = λ * 2 v n (x), M n (µ) -1 v n (y) dφ ε ξ (y) = Λµ n (ξ, ε) K µ n (x, y) dφ ε ξ (y) ,
which yields (2.8)-(2.9). Next, using the definition (2.3) of K µ n , and

P α P β dµ = δ α=β , for all α, β ∈ N d n , (2.12) Λµ n (ξ, ε) = (p * n ) 2 dµ = Λµ n (ξ, ε) 2 α∈N d n B∞(ξ,ε) P α (y) dy τ ε 2 .
To obtain (2.10), just use that with (P α (x)) α∈N d n = Dv n (x), by orthogonality of the (P α )'s w.r.t. µ, the resulting change of basis matrix

D ∈ R s(n)×s(n) satisfies D T D = M n (µ) -1 .
Finally it remains to prove that ṽn

∈ R[ξ, ε] n . With β ∈ N d n , (2.13)
B∞(ξ,ε)

y β dy τ ε = ε -d d i=1 (ξ i + ε/2) β i +1 -(ξ i -ε/2) β i +1 β i + 1
is indeed a polynomial in (ξ, ε). To see this, in each term of the product in (2.13), use the identity

a n+1 i -b n+1 i = (a i -b i ) ( n j=0 a n-j i b j i
) with a i = ξ i +ε/2 and b i = ξ i -ε/2, so that a i -b i = ε for ever i, which implies that the term ε -d is annihilated. This implies that ṽn (ξ, ε) = B∞(ξ,ε) v n (y) dy/τ ε is a polynomial vector in the variables (ξ, ε).

Computation.

• If the orthonormal polynomials (P α ) α∈N d are available then just use (2.12) to (i) compute the polynomials q α ∈ R[ξ, ε] n defined by:

x → q α (ξ, ε) := 1 τ ε B∞(ξ,ε) P α (y) dy , α ∈ N d n ,
and (ii), sum up

α∈N d n q α (ξ, ε) 2 . • On the other hand, if the moment matrix M n (µ) is available then (i) compute the polynomial vector ṽn (ξ, ε) ∈ R[ξ, ε] s(n) n in (2.11) and (ii), form the SOS polynomial ṽn (ξ, ε) T M n (µ) -1 ṽn (ξ, ε)/τ 2
ε to obtain the SOS polynomial Λµ n (ξ, ε) -1 . In both cases, computing the polynomials q α or ṽn is an easy task which can be done exactly and even symbolically in (ξ, ε), as it reduces to integrate a polynomial on the box B ∞ (ξ, ε) parametrized by ξ and ε.

Therefore it is worth emphasizing that in the end, one thus obtains the polynomial (ξ, ε) → Λµ n (ξ, ε) -1 as an explicit SOS polynomial in the variables (ξ, ε), via ṽn (ξ, ε) T M n (µ) -1 ṽn (ξ, ε), exactly as ξ → Λ µ n (ξ) -1 was obtained as an explicit SOS polynomial in ξ, via v n (ξ) T M n (µ) -1 v n (ξ).

Compared with computation of the standard Christoffel function Λ µ n , computing the extended Christoffel function Λµ n only requires an extra symbolic integration (e.g., of the vector v n (x) = (x α ) α∈N d n on the box B ∞ (ξ, ε)), an easy task once (as long as the dimension d is not too large). In addition, from inspection of (2.13) we may conclude that

ṽn (ξ, ε) = v n (ξ) if n ≤ 1, v n (ξ) + O(ε) if n > 1,
which in turn yields:

(2.14) Λµ n (ξ, ε) -1 = Λ µ n (ξ) -1 + O(ε) , ∀ξ ∈ R d . Example 1. Let Ω = [-1, 1] and µ = dx/ √ 1 -x 2
, so that the Chebyshev polynomials of first kind (T j ) j∈N are orthogonal w.r.t. µ, and the family (P j ) j∈N with P 0 = T 0 / √ π and P j = T j 2/π, j ≥ 1, is orthonormal w.r.t. µ. Then Λ µ 0 (ξ) = 1/π = Λµ 0 (ξ, ε), and

Λ µ 1 (ξ) = 2 π 1 2 + ξ 2 = Λµ 1 (ξ, ε) Λ µ 2 (ξ) = 2 π 3 2 -3 ξ 2 + 4 ξ 4 Λµ 2 (ξ, ε) = 2 π 3 2 -3 ξ 2 + 4 ξ 4 + ε 2 2 - 2 ξ 2 ε 2 3 = Λ µ 2 (ξ) + ε 2 π (1 - 4 ξ 2 3
) .

2.4.

A L 2 (µ)-norm interpretation. Recall that dµ = f dx for some unknown density f : Ω → R + , and f ≥ γ for some γ > 0. In particular, x → f (x) -1 is an element of L 2 (µ).

Lemma 2.3. Assume that f ≥ γ on Ω, for some γ > 0. Let ξ ∈ Ω and ε > 0 be fixed.

(i) If B ∞ (ξ, ε) ⊂ Ω then the linear functional ε ξ : R[x] → R, h → ε ξ (h) := 1 τ ε B(ξ,ε) h(x) dx , ∀h ∈ R[x] ,
is an element of L 2 (µ), that we still denote ε ξ . Moreover, it is represented by the function

x → 1 B∞(ξ,ε) (x) τ ε f (x) ∈ L 2 (µ) , ∀x ∈ Ω , hence with norm (2.15) ε ξ L 2 (µ) = Ω 1 B∞(ξ,ε) (x) τ ε f (x) 2 f dx 1/2 = 1 B∞(ξ,ε) τ ε f L 2 (µ)
.

(ii) In addition, if f is continuous then

(2.16) lim ε→0 τ ε ε ξ 2 L 2 (µ) = 1 f (ξ)
, ∀ξ ∈ int(Ω) .

Proof. (i) The linear functional ε ξ is a Borel measure on B ∞ (ξ, ε) with density 1/τ ε w.r.t. Lebesgue measure on B ∞ (ξ, ε).

If B ∞ (ξ, ε) ⊂ Ω and f ≥ γ > 0 on Ω, then ε ξ (h) = 1 τ ε B∞(ξ,ε) h(x) dx = Ω h(x) 1 B∞(ξ,ε) (x) τ ε f (x) f (x) dx = Ω h(x) 1 B∞(ξ,ε) (x) τ ε f (x) dµ = h, 1 B∞(ξ,ε) τ ε f , ∀h ∈ L 2 (µ) ,
and so the linear functional ε ξ is represented by

1 B∞(ξ,ε) τε f ∈ L 2 (µ). It is simply the measure on Ω with density 1 B∞(ξ,ε) /τ ε f w.r.t. µ, hence an element of L 2 (µ) with norm 1 B∞(ξ,ε) /τ ε f L 2 (µ) .
(ii) In addition, if f is continuous then

ε ξ 2 L 2 = Ω 1 B∞(ξ,ε) (x) τ ε f (x) 2 dµ = Ω 1 B∞(ξ,ε) (x) τ ε f (x) 2 f dx = B∞(ξ,ε) 1 τ ε f (x) dx τ ε = 1 τ ε f (ζ ε ) ,
for some ζ ε ∈ B ∞ (ξ, ε), which yields (2.16) when ε ↓ 0 (as f is continuous).

2.5. Relating ε ξ with Λµ n , and asymptotic properties. We first start with relating Λµ n with Λ µ n for a fixed degree "n".

Proposition 2.4. Let ε > 0 fixed, and assume that f is continuous. Then for every n ∈ N:

(2.17) Λµ n (ξ, ε) ≥ Λ µ n (ζ n ) for some ζ n ∈ B ∞ (ξ, ε).
In particular, using ε := 1/n yields that for every n ∈ N,

(2.18) Λµ n (ξ, 1/n) ≥ Λ µ n (ζ n ) for some ζ n ∈ B ∞ (ξ, 1/n). Proof. With p * n ∈ R[x] n as in Lemma 2.2, 1 = B∞(ξ,ε) p * n dx ε d = p * n (ζ n ) for some ζ n ∈ B ∞ (ξ, ε),
and so as

p * n (ζ n ) = 1, Λµ n (ξ, ε) = (p * n ) 2 dµ ≥ inf p∈R[x]n { p 2 dµ : p(ζ n ) = 1} = Λ µ n (ζ n ) ,
which yields (2.17) and (2.18).

We next consider the asymptotic limit of Λµ n (ξ, ε) as n increases, when ξ ∈ int(Ω) and ε > 0 is fixed, as well as when := 1/n r with r > 0.

Theorem 2.5. Assume that f ≥ γ on Ω, for some γ > 0, and let Λµ n be as in (2.7). Then :

(i) With ε > 0 fixed and for all ξ ∈ Ω such that B ∞ (ξ, ε) ⊂ Ω:

(2.19) lim n→∞ Λµ n (ξ, ε) -1 = ε ξ 2 L 2 = B∞(ξ,ε) 1 τ ε f dx τ ε .
Equivalently,

(2.20) lim n→∞ ε d Λµ n (ξ, ε) -1 = B∞(ξ,ε) 1 f dx τ ε = 1/f L 1 (B∞(ξ,ε),dx/τε) ,
and in addition, if f is continuous on B ∞ (ξ, ε) then

(2.21) lim n→∞ ε -d Λµ n (ξ, ε) = f (ζ ε ) , for some ζ ε ∈ B ∞ (ξ, ε). In particular, lim n→∞ ε -d Λµ n (ξ, ε) ≈ f (ξ) when ε is small.
(ii) For all ξ ∈ int(Ω) and all n ∈ N,

n d Λµ n (ξ, 1/n) ≥ inf x {f (x) : x ∈ B ∞ (ξ, 1/n)} , so that if f is continuous on B ∞ (ξ, ε) then, (2.22) lim inf n→∞ n d Λµ n (ξ, 1/n) ≥ f (ξ) .
Moreover, if r < 1 then

(2.23) lim n→∞ n d Λµ n (ξ, 1/n r ) = +∞. Proof. (i) By Lemma 2.3, ε ξ ∈ L 2 (µ). To prove (2.19), notice that ε ξ 2 L 2 = sup h∈L 2 (µ) ε ξ (h) 2 h 2 dµ ,
and therefore,

1 ε ξ 2 L 2 = inf h∈L 2 (µ) h 2 dµ ε ξ (h) 2 = inf h∈L 2 (µ) { h 2 dµ : ε ξ (h) 2 = 1} = inf h∈L 2 (µ) { h 2 dµ : ε ξ (h) = 1} ≤ inf h∈R[x]n { h 2 dµ : ε ξ (h) = 1} , ∀n ∈ N ≤ Λµ n (ξ, ε) , ∀n ∈ N , (2.24)
and therefore, with ε > 0 fixed,

(2.25) lim inf n→∞ Λµ n (ξ, ε) ≥ 1 ε ξ 2 L 2 , for all ξ ∈ Ω s.t. B ∞ (ξ, ε) ⊂ Ω .
On the other hand, and still with ε > 0 fixed, as B ∞ (ξ, ε) ⊂ Ω, let h ε :=

1 B∞(ξ,ε) τεf ∈ L 2 (µ). As Ω is compact, R[x] is dense in L 2 (µ) and therefore there exists a sequence (h n ) n∈N ⊂ R[x] such that h n -h ε L 2 (µ) → 0 as n increases. Moreover, φ ε ξ (|h n -h|) = B∞(ξ,ε) |h n -h ε | dx τ ε = B∞(ξ,ε) |h n -h ε | τ ε f f dx = Ω 1 B∞(ξ,ε) |h n -h ε | τ ε f dµ ≤ µ(Ω) τ ε f * Ω (1 B∞(ξ,ε) (h n -h ε )) 2 dµ 1/2 ≤ µ(Ω) τ ε f * Ω (h n -h ε ) 2 dµ 1/2 = µ(Ω) τ ε f * h n -h ε L 2 (µ) ,
where

f * = min {f (x) : x ∈ Ω} (> 0 as f is continuous and f > 0 on Ω compact). Therefore, φ ε ξ (h n ) → φ ε ξ (h) as n increases. Then from Λµ n (ξ, ε) ≤ h 2 n dµ ( h n dφ ε ξ ) 2 , ∀n , and since lim n→∞ φ ε ξ (h n ) = φ ε ξ (h), lim n→∞ h 2 n dµ = h 2 dµ, we conclude that for all ξ ∈ Ω with B ∞ (ξ, ε) ⊂ Ω, lim sup n→∞ Λµ n (ξ, ε) ≤ lim n→∞ h 2 n dµ ( h n dφ ε ξ ) 2 = h 2 dµ ( h dφ ε ξ ) 2 = 1 ε ξ 2 L 2
, which, combined with (2.25), yields (2.19). Finally, as f > 0 on Ω and f is continuous, B∞(ξ,ε)

1 f dx τε = f (ζ ε ) -1 for some ζ ε ∈ B ∞ (ξ, ε), which yields (2.20). Moreover, f (ζ ε ) → f (ξ) as ε ↓ 0 and so f (ζ ε ) ≈ f (ξ) for > 0 fixed and small.
(ii) Finally let ξ ∈ int(Ω) and let n := n -r with r > 0 so that as n increases, B ∞ (ξ, 1/n r ) ⊂ Ω for all n sufficiently large, and τ εn = ε d n = n -rd . Recall that with n fixed sufficiently large,

εn ξ 2 L 2 = 1/(τ n f (ζ n )) = n rd /f (ζ n ) , for some ζ n ∈ B ∞ (ξ, 1/n r ).
Then for each n fixed, and with same arguments as in (2.24), one obtains

n -rd f (ζ n ) = 1 εn ξ 2 L 2 = inf h∈L 2 (µ) h 2 dµ εn ξ (h) 2 = inf h∈L 2 (µ) { h 2 dµ : εn ξ (h) = 1} ≤ inf h∈R[x]n { h 2 dµ : εn ξ (h) = 1} ≤ Λµ n (ξ, 1/n r ) , (2.26)
from which we deduce that with r = 1, n d Λµ n (ξ, 1/n) ≥ f (ζ εn ), for every n ∈ N, and thus

n d Λµ n (ξ, 1/n) ≥ inf{f (x) : x ∈ B ∞ (ξ, 1/n r )} , ∀n ∈ N . If f is continuous then f (ζ n ) → ξ as n → ∞,
and we obtain (2.22). Finally,

n d Λµ n (ξ, 1/n r ) ≥ n d(1-r) f (ζ εn )
, for every n ∈ N, which yields (2.23). So Theorem 2.5 states that for fixed ε > 0, asymptotically when n increases, Λµ n (ξ, ε) -1 converges to the L 2 (µ)-norm of the linear functional ε ξ ∈ L 2 (µ) which can be viewed as an ε-average (or ε-regularization) of the point evaluation δ ξ ∈ L 2 (µ). However, and as expected, ε ξ L 2 increases with 1/ε.

On the other hand, for the standard Christoffel function, Λ µ n (ξ) -1 is the L 2 (µ)-norm of the the point evaluation linear functional δ ξ viewed as an element of the finite dimensional subspace

(R[x] n , •, • L 2 ) of L 2 (µ). But then of course this norm is unbounded as n increases because δ ξ ∈ L 2 (µ). Remark 2.6. In Theorem 2.5(i) one restricts to ξ ∈ int(Ω) such that B ∞ (ξ, ε) ⊂ Ω. The reason is that if ξ ∈ int(Ω) but µ(B ∞ (ξ, ε)∩(R d \Ω)) > 0, then for all h ∈ L 2 (µ), B∞(ξ,ε) h dx ε d = Ω h 1 B∞(ξ,ε)∩Ω ε d f dµ = h , 1 B∞(ξ,ε)∩Ω ε d f ,
and therefore, (2.20) becomes

(2.27) lim n→∞ ε d Λµ n (ξ, ε) -1 = B∞(ξ,ε)∩Ω 1 f dx ε d , while if f is continuous then (2.21) becomes lim n→∞ ε -d Λµ n (ξ, ε) = f (ζ ε ) τ ε vol(B ∞ (ξ, ε) ∩ Ω) ,
for some ζ ε ∈ B ∞ (ξ, ε) ∩ Ω. So now the limit in (2.27) is the product of -the nice term f (ζ ε ) related to the density of f in B ∞ (ξ, ε), with -the multiplicative term ε d /vol(B ∞ (ξ, ε) ∩ Ω), which may be hard to compute for arbitrary Ω.

On the other hand, the limit in (2.27) still converges to f (ξ) when ε ↓ 0.

Comparing with asymptotics of the standard

Λ µ n . Recall that s(n) = d+n d = O(n d ).
Notice that under some restrictive conditions, if ξ ∈ int(Ω),

(2.28)

lim n→∞ s(n) Λ µ n (ξ) = f (ξ) ω E (x) ,
where ω E (x) is the (in general unknown) density of the equilibrium measure of Ω.

(ξ, ε) = f (ζ ε ) , for some ζ ε ∈ B ∞ (ξ, ε), when ε > 0 is fixed, B ∞ (ξ, ε) ⊂ Ω,
(ξ, 1/n) ≥ f (ξ) ; lim n→∞ n d Λµ n (ξ, 1/n r ) = +∞ if r < 1 ,
when ε = 1/n or ε = 1/n r is not fixed anymore. So (2.28) is more precise but requires strong assumptions on f , and is not very helpful for practical computation since the density ω E of the equilibrium measure of Ω is not known in general.

On the other hand, (2.21) is obtained with a rather weak assumption on f , and provides an information of the density f in the box B ∞ (ξ, ε), as well as an approximation of f (ξ) when ε is small. However it holds for points ξ ∈ int(Ω) whose • ∞ -distance to the boundary ∂Ω is at least ε/2 > 0.

Notice that as indicated in (2.23), varying ε with n as ε = 1/n r , with r < 1, is not a good idea if one would like to compare the limit of n d Λµ n (ξ, 1/n r ) with that of n d Λ µ n (ξ) in (2.28). To obtain a meaningful bounded limit would require at least ε = 1/n. Support inference. As already mentioned in the introduction, one striking property of the Christoffel function Λ µ n is the dichotomy of its asymptotic behavior with n, depending on whether ξ ∈ int(Ω) or ξ ∈ Ω. As we next show, Λµ n inherits the same highly desirable property.

1. Inside Ω. First with ε > 0 fixed, by Theorem 2.5(i), for every ξ such that B ∞ (ξ, ε) ⊂ Ω,

lim n→∞ ε -d Λµ n (ξ, ε) -1 → B∞(ξ,ε) (1/f ) dx τ ε ,
so that the sequence ( Λµ n (ξ, ε)) n∈N is bounded (as it converges). So clearly, ( Λµ n (ξ, ε)) -1 is bounded, uniformly in n. Next, when ε is allowed to vary with n, like ε = 1/n r with r > 0, then from (2.26) in the proof of Theorem 2.5,

Λµ n (ξ, 1/n r ) -1 ≤ n rd γ , ∀ξ ∈ int(Ω) ,
that is, the growth of Λµ n (ξ, 1/n r ) -1 is at most polynomial in n.

2. Outside Ω.

Lemma 2.7. Let ξ ∈ Ω and fix ε > 0. Then the growth of Λµ n (ξ, ε) -1 is at least exponential in n. Similarly, the growth of Λµ 2n (ξ, 1/n) -1 is at least exponential in n.

A proof is postponed to §4.2.

So, similarly as for the standard Christoffel function Λ µ n (ξ), we also obtain for Λµ n (ξ, ε) a dichotomy in its asymptotic behavior as n increases, depending on whether ξ ∈ Ω (exponential decrease to zero) or ξ ∈ Ω and at • ∞distance ε of the boundary ∂Ω (convergence to a nonzero value). A similar result holds for Λµ 2n (ξ, 1/n) when > 0 is allowed to vary with n, like e.g., ε = 1/n.

Conclusion

We have introduced an extended version Λµ n (x, ε) of the standard Christoffel function Λ µ n (x) for a measure µ on a compact set Ω ⊂ R d . One main motivation was to improve the asymptotic behavior of Λ µ n (x) as n increases and x ∈ int(Ω), from a computational point of view. Indeed, when s(n)Λ µ n (x) converges (under certain assumptions), the limit is of the form f (x)/ω E (x) where f (resp. ω E ) is the density of µ (resp. of the equilibrium measure of Ω) w.r.t. Lebesgue measure. As ω E is in general unknown, this result is of limited value from a computational viewpoint. In contrast, if f is continuous, then with ε > 0 fixed and as n grows, ε

-d Λµ n (x, ε) → f (ζ ε ) for some ζ ε ∈ B ∞ (x, ε). So if ε is small,
one obtains an approximation of f (x) which converges to f (x) if ε ↓ 0. In addition, the extended Christoffel function still inherits highly desirable properties of 1/Λ µ n , namely (i) the dichotomy of its symptotic behavior (at most polynomial growth versus at least exponential growth) depending on whether x is inside or outside Ω, and (ii) an efficient computation of its closed form expression as a polynomial of (x, ε). Next, for every δ ∈ (0, 1), the needle polynomial q ∈ R[t 2 ] n(δ) introduced in Króo and Lubinsky [START_REF] Kroó | Christofffel functions and universality in the bulk for multivariate polynomials[END_REF], satisfies: (4.4) q(0) = 1 ; |q( x 2 )| ≤ 2 • 2 -δ n(δ) ∀x ∈ B 2 (0, 1) \ B 2 (0, δ) , and in particular, q( x 2 ) 2 ≤ 4 • 2 -δ 2n(δ) , ∀x ∈ T (Ω) , whenever B 2 (0, δ)∩T (Ω) = ∅ (which happens whenever δ is sufficiency small as 0 ∈ T (Ω)). So with q as in (4.4), letting v = q 2 ∈ R[x] 2n(δ) , one obtains v(0) = 1 and

B∞(0,ε/r) v(x) dx τ ε/r ≥ C -1 η (1 + η) -n(δ) .
Therefore, we deduce that Λν 2n(δ) (0, ε/r) ≤ T (Ω) v 2 dν B∞(0,ε/r) v dx/τ ε/r

≤ C η (1 + η) n(δ) T (Ω) v 2 dν ,

2. Main result 2 . 1 .

 21 Notation and definitions. Let R[x] denote the ring of real polynomials in the variables x = (x 1 , . . . , x d ) and R[x] n ⊂ R[x] be its subset of polynomials of total degree at most n. Let N d n := {α ∈ N d : |α| ≤ n} (where

4 . Appendix 4 . 1 .Lemma 4 . 1 . 2

 441412 Affine invariance of the Christoffel function. Let Ω ⊂ R d be compact and let ξ ∈ R d be fixed. With r > 0, consider the linear transformation T : R d → R d defined by (4.1)x → T (x) := (x -ξ)/r , x ∈ R d .If µ is a Borel measure on Ω, let ν be the pushforward measure of µ by the mapping T , that is,ν(C) := µ(T -1 (C)) , ∀C ∈ B(Ω) .The support of ν is the set T (Ω), andT (Ω) p(y) dν(y) = Ω p(T (x)) dµ(x) , ∀p ∈ R[x] .Let T be as in (4.1) and with µ supported on Ω, let ν be its pushforward by T (hence supported on T (Ω)). ThenΛ ν n (T (x)) = Λ µ n (x) , ∀x ∈ R d (4.2) Λν n (0, ε/r) = Λµ n (ξ, ε) , ∀ξ ∈ R d . ) 2 dν(y) : p(T (x)) = 1 } = min p∈R[y]n { Ω p(T (y)) 2 dµ(y) : p(T (x)) = 1 } = Λ µ n (x) ,where we have used that {p • T : p ∈ R[x] n } generates R[x] n . Moreover, observe that B(0,ε/r) p(u) du (ε/r) d = B∞(T (ξ),ε) p(T (x)) dx ε d , and since T (ξ) = 0, we also obtain Λν n (0, ε/r) (T (x)) dµ(x) : B∞(ξ,ε) p(T (v)) dv ε d = 1 } = Λµ n (ξ, ε) .

4. 2 .

 2 Proof of Lemma 2.7.Proof. Let T : R d → R d be as in (4.1) with r > 0 such that T (Ω) ⊂ B 2 (0, 1) (where B 2 (ξ, τ ) := {x : xξ 2 < τ }), and let the measure ν on T (Ω) be the pushforward of µ by T . Observe that T (ξ) = 0 ∈ T (Ω) since ξ ∈ Ω. By the Bernstein-Markov property of Lebesgue measure on the box B ∞ (0, ε/r), for every η > 0 there exists a constant C η > 0, such that for all p ∈ R[x],sup{|p(x)| : x ∈ B ∞ (0, ε/r)} ≤ C η (1 + η) deg(p) B∞(0,ε) p(x)2 dx τ ε/r , see [8, p. 51]. Therefore if p(0) = 1 then, B∞(0,ε/r) p(x) 2 dx τ ε/r ≥ C -1 η (1 + η) -deg(p) .

This work is supported by the AI Interdisciplinary Institute ANITI funding through the french program "Investing for the Future PI3A" under the grant agreement number ANR-19-PI3A-0004. This research is also part of the programme DesCartes and is supported by the National Research Foundation, Prime Minister's Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme. 1

whenever B 2 (0, δ) ∩ T (Ω) = ∅, and therefore Λν 2n(δ) (0, ε/r) ≤ 16 C η (1 + η) n(δ) 2 -δ 4n(δ) . So recalling that ξ ∈ Ω (hence 0 ∈ T (Ω)), fix δ ∈ (0, 1) such that B 2 (0, δ) ∩ T (Ω) = ∅. Then by choosing η such that a := (1 + η)/2 4δ < 1, e.g. η < 4δ ln 2, one obtains:

which in turn shows that Λµ 2n (ξ, ε) decreases exponentially fast to zero with n, whenever ξ ∈ Ω.

Next, if ε = 1/n for every n, then agin let δ > 0 be such that B 2 (0, δ) ∩ T (Ω) = ∅. As we did above with ε fixed,

. Again, by choosing η such that a := (1 + η)/2 4δ < 1, e.g. η < 4δ ln 2, one obtains:

which shows that Λµ 2n (ξ, 1/n) decrease exponentially fast to zero whenever ξ ∈ Ω.


domain: phys spi

text:
R Guibert

A Bayle

F Plouraboué

Geolocalization of water-waves origin within water distribution networks using time reversal of first event detection

Keywords: water distribution network, water-hammer, high-frequency detector, geolocalization 2010 MSC: 00-01, 99-00

Drinking water distribution networks in urban areas are daily -based subjected to fast propagating pressure waves resulting from routine operations. These water-hammer waves leads to structural aging and facility damages, the origin of which is not easy to find , but are sometimes of high managerial interest. In this contribution, we demonstrate that using a reasonable number of high-frequency pressure detectors distributed within the network combined with a proper post-processing method permits a close geolocalization of the damaging wave origin. The method is first tested and validated on a real water distribution network having approximately 26 000 pipes, whereas considering a known, prescribed wave-origin, so that the sensitivity to sensor number (sensor spatial density), sensor location and signal-to-noise ratio on the geolocalization robustness is are analyzed in detail. It is then applied and illustrated over real sensor recordings the result of which are validated on the field from history matching. This paper thus presents the first field-scale geolocalization of water-hammer events origin test as well conditions for which, given sensor density and signal-to-noise ratio, the geolocalization success is to be expected.

Introduction

Water hammer waves within Water Distribution Networks (WDNs) are a common hindrance resulting from many possible routine operations, either organized, volountary or accidental. Not only these waves are responsible of for structure's wear, but they are also capable of facility damages. Intermittent water supply operations have indeed been correlated to incident damages in WDNs (Agathokleous and Christodoulou, networks within a few seconds before decaying in less than a minute. It results in a myriad of reverberation causing as much disorder and possible breakages as pre-existing weaknesses. Since the life-time of these events inside the network is short, and since they are not numerous within a day (e.g less than a few per day in the considered city of about half a million people), the dynamic of each water-hammer event can be considered separately. Being able to find the origin of a cascade of reverberating water-waves events is interesting from various viewpoints: patrimonial management, security, insurances, etc... Nevertheless, this topic has not been developed in real WDNs because it necessitates to overcome overcoming several challenges, not been fully addressed yet. Firstly, urban WDNs are, in many cases, not sufficiently and reliably detected and observed in details for a precise modeling of transient waves into them to be relevant. Secondly, localizing the origin of water-hammer waves necessitates a real-time high-frequency monitoring with widely distributed detectors, rarely deployed in WDNs. Thirdly, the computational cost of direct transient wave modeling in water distribution networks is very challenging either using the Method Of Characteristics [START_REF] Wang | Pipe burst risk state assessment and classification based on water hammer analysis for water supply networks[END_REF][START_REF] Nault | Generalized flexible method for simulating transient pipe network hydraulics[END_REF][START_REF] Meniconi | Leak detection in a real transmission main through transient tests: Deeds and misdeeds[END_REF][START_REF] Moosavian | Unified matrix frameworks for water hammer analysis in pipe networks[END_REF][START_REF] Riaño-Briceño | Distributed and vectorized method of characteristics for fast transient simulations in water distribution systems[END_REF], or finite volumes e.g [START_REF] Pal | An overview of the numerical approaches to water hammer modelling: The ongoing quest for practical and accurate numerical approaches[END_REF][START_REF] Zhang | Algorithm for detecting multiple partial blockages in liquid pipelines by using inverse transient analysis[END_REF]. Fourthly, an inverse method capable of identifying the origin of a reverberating water hammer wave over a large network, solving as many direct problems as tentative tries for possible origin, is obviously even more challenging. Last but not least, from a more fundamental viewpoint, even if the uniqueness of the wave origin from detecting reverberation waves has been established on discrete wave models on arbitrary networks [START_REF] Caputo | Inverse source problem in a forced network[END_REF][START_REF] Plouraboué | Source identification of propagating waves inside a network[END_REF], it is still a pending issue for continuum ones.

Hence, even if the idea of using the entire time-course of signal reverberation within the network at the sensor location so as to enlarge as much as possible data collection is appealing, in practice, this approach is still very difficult to develop at the present state of the art [START_REF] Che | Transient wave-based methods for anomaly detection in fluid pipes: A review[END_REF]. As part of the European "Surge-Net" project, [START_REF] Ferrante | Leak detection in branched pipe systems coupling wavelet analysis and a Lagrangian model[END_REF] carried out leakage location tests using echo analysis, by combining Lagrangian and wavelet transforms signal analysis, in the Lintrahen's (Scotland) WDN main trunk. [START_REF] Shucksmith | Onsite leak location in a pipe network by cepstrum analysis of pressure transients[END_REF] performed leakage tests on the Bradford's (Yorkshire, UK) WDN. The authors operated at the neighbourhood scale (about 100m of weakly branched pipe) using spectral analysis for the leak echo-localization (cepstrum analysis) and with a wide variety of pipe materials (PVC, cast iron, asbestos-cement). Recently, [START_REF] Meniconi | Anomaly pre-localization in distribution-transmission mains by pump trip: Preliminary field tests in the Milan pipe system[END_REF] also combined a wavelet transform analysis with a Lagrangian method to preliminary examine one of the main pressure transimission lines of the Milan's (Italy) WDN. Although achievable on a district scale or on the main supply pipes of a water distribution network, the implementation of such leak detection techniques is questionable for large-scale analysis. Indeed, techniques based on leak wave reflection suffer from a high sensitivity to propagation speed uncertainties and surrounding noise [START_REF] Wang | Pipe burst risk state assessment and classification based on water hammer analysis for water supply networks[END_REF]. Moreover, an idealized test case without leakage or blockage is required for allowing direct comparison with the field pressure signal. Recently, Meniconi and co-workers have performed such measurements on a laboratory network (Meniconi et al., 2022a,b). However, such benchmark case is often not operationally possible or, alternatively prone to modelling errors (e.g. crosssection pressure interactions, wave attenuation, incorrect network topology, etc...) as in other contexts such as quality event detection [START_REF] Kessler | Detecting accidental contaminations in municipal water networks[END_REF]. This is why an approximate alternative method avoiding the cost of a complete time-course of wave propagation evaluation is proposed and developped in this contribution. Rather than trying to exploit the entire signal complexity recorded at each sensor location, we take advantage of the first arrival time only.

This already permits to back-track the wave origin using a time-reversal method. Time reversal methods, either at the individual pipe level [START_REF] Grigoropoulos | Time reversal of waves in hydraulics: experimental and theoretical proof with applications[END_REF] or at the network one [START_REF] Shen | Locating the source of diffusion in complex networks by time-reversal backward spreading[END_REF] have already permit computationaly efficient source identification. This method has been used to develop noisetolerant pipeline defect detection [START_REF] Wang | Pipeline leak detection using the matched-field processing method[END_REF][START_REF] Wang | A review on applications of remote sensing and geographic information systems (gis) in 330 water resources and flood risk management[END_REF][START_REF] Meniconi | Leak detection in a real transmission main through transient tests: Deeds and misdeeds[END_REF] in the precise context of water-hammer waves. The aim of this contribution is This contribution aims to demonstrate that combining first-event detection with time-reversal is an interesting strategy to geolocalize water-hammer events in real WDN.

The paper is organized as follows. Section 2.1 describes the material related to the urban network and the pressure sensors used within the study. Section 2.2 details the algorithm of the first event detection, time reversal method, and candidate ordering. Section 3.1 discusses the localization of prescribed events in a real WDN so as to test and validate the method efficiency and its sensitivity to detectors number (more precisely the sensor "density", the sensor number divided by the total node number in the network). Section 3.2 finally describes the application of the proposed method to the field's data so as to geolocalize real events and discuss the obtained results.

Materials and Methods

Materials

The materials consist in a WDN data within which high-frequency pressure detectors are disposed of.

Water distribution network data

The entire network from the city of Toulouse (France) illustrated in figure 1a is composed of a total of 26 094 pipes. It has been formatted in Epanet format providing the pipe lengths, structural properties and diameters. The distinct pipe's composition within the network is detailed in table 1, showing a great majority of cast iron material. The heterogeneity of pipe diameters and lengths is illustrated in figure 1b and c showing a great diversity of diameter and length over more than a decade. each pipe. Local wave-hammer wave velocity indeed depends on the pipe diameter, length, thickness, Young modulus, and material density. It is estimated from known validated formulas (cf Appendix). Furthermore, the chosen location for high-frequency pressure sensor's (18 in total) is depicted in Figure 1a with green dots, each described by a label between 1 to 18, being spatially uniformly distributed over the urban network.

High-frequency pressure detectors

Pressure detectors record at 128 Hz frequency with a 5000P a accuracy. The recording mode is triggered by an awaking threshold which empirically sets from the base signal. This precautionary procedure avoids recording embarrassing irrelevant signal series in the detection event database. Since the event amplitude is expected within the range of [0.8, 2]10 5 Pa, the relative precision on the detected peak amplitude is of the order of 6.25%. Also, only the time-arrival of the first peak is critical for the chosen method. Each detector has been primarily submitted to a pressure calibration test using a prescribed static pressure before field deployment. Each sensor has its own embarked battery and RAM. Once triggered to awaken mode, the detectors acquire at full high-frequency (i.e. 128Hz) during 300 seconds. The resulting collected data are then more lately transmitted latelly with a standard GSM 3G protocol keeping each signal associated with each detector identifier.

Methods

The method decomposes into various steps : (i) first event detection into each detectors, (ii) backpropagation from detectors to potential source of detected event using the time-reversal method and (iii) calibration of time-reversal method on the network and detector set. The methods associated with these three steps are now detailed. An additional noise sensitivity check of the method has also been used, the details of which is also given in 2.2.3.

First event detection and time reversal method

First event detection is performed using offline change point detection method [START_REF] Truong | Selective review of offline change point detection methods[END_REF]. The detection criteria isare associated with a functional minimization associated with the local gradient of the noisy signal. More precisely we use a binary change point detection to perform a fast signal segmentation, coupled with a L 2 cost function that detects mean-shifts in the signal. The determination of the arrival times is illustrated for a real signal in figure 2. 1a) using high-frequency detectors. The illustrated event is discussed in subsection 3.2. For each signal related to the event, the normalized pressure is reported versus time (cf figure 7). The obtained arrival time is represented by the vertical dashed line.

The signal is captured almost simultaneously at detectors #13 and #10 (respectively at time 112.92 and 112.96 s), then at detector #12 (at time 114.29 s), and finally at detector #18 (115.00 s).

Time-reversal of first event detection within the network follows the approach proposed in Shen et al.

105 (2016). The method principle is detailled in Figure 3. As a prerequisite, each pipe is associated with a time propagation resulting from computing the ratio between pipe lengtht and velocity evaluated following the Appendix formula. Then the signal processing starts when an event has been detected, resulting in a first-time arrival (denoted t i , i = 1, 2, 3 in figure 3a) of the signal at various awaken detectors (denoted d i , i = 1, 2, 3). For each pair of detector/non-detector nodes in the network, a total-time is evaluated from 110 by computing the sum of each time propagation within each pipes along the shortest path between those within the network, as depicted in figure 3b. Performing this total time of propagation between each node and one detector results in the detector back-propagation cartography depicted in figure 3c for detector d 1 , figure 3d for detector d 2 and figure 3e for detector d 3 . Now using the first-time arrival component vector color as the detector with which it is associated in figure 3f). In [START_REF] Shen | Locating the source of diffusion in complex networks by time-reversal backward spreading[END_REF] the source node is the one with minimum variance back-propagating-time vector as illustrated in figure 3f. Note, however, that depending on the recorded time, several source nodes are possibly found with this method as depicted in figure 3e. Hence, in the case of noisy recordings one can infer that, in this case, the true source might not necessarily be the one having the minimum variance. Hence, in order to give more robustness to the method 120 (but obviously less sensitivity) we extend the search for the true source fromby considering the sorted list of back-propagating-time vector variances in order to find the "best" source candidates. This sorted list has to be closed, keeping with a "reasonable" number of possible candidates. This is what we call the method "calibration" as detailed in the next section.

Calibration of the time reversal method

In the original research [START_REF] Shen | Locating the source of diffusion in complex networks by time-reversal backward spreading[END_REF], the locabilitylocality condition is widely discussed. Two main aspects are highlighted: (i) all edges spreading times must be sufficiently different and (ii) nodes with a single neighbor disrupt this condition. Consequently, the success rate (i.e. exact source location) of the method depends on the topology of the considered network and the number of detectors. In our approach, the signals are real and uncertainties necessarily exist in the signal measurements. Moreover, it is currently not possible to deploy sensors on 20% of the nodes in the city networks, as suggested by [START_REF] Shen | Locating the source of diffusion in complex networks by time-reversal backward spreading[END_REF]; the current order of magnitude is less than one percent. In response to this, the method has been adapted to achieve interesting success rates with a limited number of sensors. The resulting variance of each source candidate areis sorted in order to produce a tentative hierarchy of the best source candidate. Nevertheless, this priority list is not always relevant: the source having the minimum variance, i.e. the first source in the sorted list of variances, is not always the true source. Hence, a list of potentially successful sources, i.e. "the best choice"

among the entire variance list has to be defined. This "best choice" needs a potentially successful source number to be defined. For this, a dedicated "calibration" of the method needs to be performed so as to find the most sensible potentially successful source number. Obviously, This potentially successful source number has to be chosen from a trade-off between accuracy and selectivity, i.e. maximizing the probability of finding the true source whilst, on the order hand, not increasing too much the number of potential sources. This "calibration" is highly sensitive to the specific network at hands, the detector number, and their positions, as well as to the wave velocity variability among the pipes. Hence, this calibration is empirically evaluated using the real network data and the exact detector positions and numbers. For this, we randomly chose sources in the network, perform the time-reversal approach in each case, and built the success-rate histogram of finding the true source at each rank of the variance list. This is what we called the "calibration" of the method. One calibration example is illustrated in figure 4a where all tested sources (used for building statistics) have been colored in red (they represent 5% of the total network). Figure 4b depicts the probability of each candidate to be being a true source versus its rank in the potentially successful source list. The closest to one, the most probable being the true source. Once adding together all potentially successful source probability in the list result in the probability of having the true source versus rank, i.e. versus the chosen potentially successful source number. The 0.9, 0.95, and 0.99 probabilities have been depicted with vertical dotted lines in figure 4b. They are "calibrated" in this case by a potentially successful source number of 11, 23, and 75 respectively.

Time reversal method with noisy data

In order tTo test the influence of noise in the first event detection time, some noise is added to the celerity wave in each pipe. This noise modelsis modeling the uncertainty associated with structural or geometrical parameters both influencing the wave velocity. We chose to impose a multiplicative noise, i.e. a noise being a small fraction of the local value of the wave celerity in each pipe. In this framework, first event detection times become random variables to be evaluated a large number of times in order to access the reliability and sensitivity of the source identification to noise. Obviously, for a given configuration (i.e. fixed network topology, wave celerity field, detector number, and locations) the calibration procedure is performed only once, so that for each set of random sources, the noise distribution is applied in a randomly and differently way. Increasing the noise amplitude permits to quantifying the degradation of the method performance, so as being to be able to extrapolate it for increasing noise and/or uncertainty.

Results

Influence of detector density and signal-to-noise ratio

As discussed in section 2.2.2, the success-rate, i.e. the capability of finding the true source among the potentially selected ones depends both upon the sensor spatial density and on the signal-to-noise ratio. We hereby analyze both from randomly chosen detector locations within the network. The detector density influence is first analyzed without noise, i.e. supposing that the wave velocity is perfectly evaluated within each pipe. In this case, 30 randomly chosen detector configurations are analyzed for each varying number of detectors between 10 to 50. This corresponds to a detector "density" (i.e. the percentage of the detector within the total node number of the network) between 0.05% to 0.22%. Note that these detector densities have been voluntarily chosen in a much modest range than the ones tested in [START_REF] Shen | Locating the source of diffusion in complex networks by time-reversal backward spreading[END_REF], in line with what is currently deployed and what will be deployed in the near future soon. This is motivated by the practical constraint that detectors deployment within the network is costly. Hence it is required to test how the method performs for weak detector density. For such low detector density, it is not expected that the method could provide a highly reliable success rate. Nevertheless, Figure 5a shows that the success rate associated with the best source in the variance list can reach 45.5% of exact detection for the highest detector density (with 50 detectors). This is already an interesting performance from the practical point of view considering the poor detector deployment effort. Moreover, when relaxing/extending the identification success to a one-neighbor distance within the network, i.e. considering one node distance to the best source of the variance list is a success, then the success rate can reach 68.7% as illustrated in Figure 5a. This means that the obtained geolocalization of the best source positions in the network is useful to locate the true source in its vicinity in case of an unsuccessful search. Hence, it is interesting to relax the selectivity of the method (which is the maximum for a single choice associated with the best variance in the potentially successful source list) in order to study how increasing the sensor density with a given calibration (from varying the list success rate from 0.9 to 0.99) impacts its selectivity. This calibration sensitivity to detector density are is is shown in figure 5b. Whereas increasing the detector density had a weak impact on the success rate of the best source candidate in the list, on the contrary, one can observe in figure 5b that it has a strong impact on the potentially successful source number which strongly decays as the detector number increases from 10 to 50. This is true for every level of calibration (i.e. for every prescribed success rate from 0.9 to 0.99). Hence, even if the method accuracy is poorly sensitive to sensor density, it is selectivity is strongly influenced by it. For this reason, we also evaluate the convex hull of potentially successful sources positions so as to provide a Region Of Interest (ROI) for the source location. This convex hull is not only useful from the practical viewpoint to delineate the research ROI. It is also usefull as a possibly expandable region in order to increase research success, in case of unsuccessful researches within it. This convex -hull will be exemplified in two practical cases in the next section.

We now turn to the analysis of noise impact on success rate. In order To analyze the expected precision of the geolocalization produced by the first event detection method in the presence of noise, we analyze the performance of the method with the same number of detectors as the one used on the field (18) for a large number of randomly chosen sources when considering a relative uncertainty between 0 to 5% for each local wave velocity in pipes. As expected, the larger the noise amplitude, the lower the success rate. Interestingly enough, the performance degradation scales linearly with the noise amplitude as found in figure 6a and6b.

Hence this permits us to predict the method performances in different contexts for which the uncertainty of the wave velocity are is larger.

Real events geolocalization

The method is now applied to two real field's data analyies denoted case 1 and case 2. For case 1, figure 7a shows the time variation of the high-frequency pressure signal from the awaken detectors (four of those, # 10, 12, 13, 18 whose location is given in figure 1) of one event. The geolocalization of the source is provided in figure 7b using various convex-hull of potentially successful sources positions associated with various calibrations as previously explained in section 3.1. As quantified in table 2, the ROI area provided by the convex hull represents a small fraction (a few percent) of the total network area. Obviously, depending on the calibration level, this area varies. This area varies depending on the calibration level. (yellow) 0.75 %, (orange) 0.5 %, ansd (red) 0.1 % of potentially successful source number in overall nodes. The four detectors that awaken during this event are depicted in green. (c) Zoom on ROI with the largest hull in dashed lines and nodes colorized versus their variance level from yellow to red for increasing probability of being the source. It decreases as the potentially successful source number used also levels -down. The proportion of pipe length within the convex -hull for the case 1 is also provided in table 2. It shows that (in the less selective choice) only 0.6% of the total pipe length are is worth of investigationing for finding the source. This

shows that the methods already permits a huge screening over the total pipe length inside the network, a result of high operational interest. This is the geolocalized illustration of the method's selectivity. Figure 7c also provides a zoom into the ROI region in order to better visualize the convex -hull as well as the true source location. In the considered event, history matching has indeed permitted us to find the true location of the responsible event associated with a pipe breakage. The black arrow of figure 7c indicates that the breakage location indeed lies inside the ROI prediction of the method. This a posteriori confirmation is a supplementary demonstration that the proposed method is relevant for the field's water-hammer source identification.

Another geolocalization example associated with case 2 is illustrated in figure 8 for which the event responsible for the water-hammer wave has been identified as the start-up of one boost pump inside a drinking water production plant. In this case height detectors have been awakened (# 1, 2, 3, 4, 5, 6, 7, 18) all showing a sudden uprising from base-state illustrated in figure 8a as opposed to the sudden depression found in figure 7a. This is indeed the expected behavior from an pressure injection event associated with a boost pump turn-on. The location of the pump is however found (figure 8c arrow) at the very edge of the ROI in this case. This might be related to the presence of high-diameter pipe connections of the drinking water production plant into the distribution system for which the associated wave velocity is perhaps roughly approximated. Concerning the quantitative figures obtained for the case 2 provided in table 2 a similar conclusion to case 1 can be raised. Also, in this case, the methods permits a huge screening over the total pipe length worth investigating from the overall network.

Conclusions

This work has investigated the geolocalization of water-hammer events sources within a water distribution network from the use of high-frequency pressure detectors distributed within it. Combining first event detection with a time-reversal method, the accuracy and selectivity of the method have been analyzed within real network configurations, considering the very low density of pressure detectors. In this context, we demonstrate from a dedicated calibration procedure the relevance of the proposed method to perform a very good screening of potentially successful sources. The effect of noise either associated with the detector signal or the network uncertainty have has also been analyzed. Performance degradation of the method have has also been quantified for a noise range between 1 to 5% of the base signal. Finally, the relevance of the proposed method is illustrated in two field cases for which history-matching analysis provides a true source location consistent with the method's prediction.

Hence the proposed method has demonstrated the proof-of-concept that water-hammer events geolocalization are is possible with a combination of rather sparing computational and technological tools. Considering the possible managerial interest of such localization in insurance issues and/or repairing investigations, this contribution can lead to significant operational consequences on in the field.

Obviously, Uncertainties on about the pipe's properties as well as network topological reliability could hamper a direct application of the method in some urban networks. Nevertheless, some more involved sensors and/or a more elaborated use of the signal beyond first-event detection might be interesting to develop in the future in order to partly overcome these uncertainties. These expressions can be found in [START_REF] Zhang | FSI Analysis of Liquid-Filled Pipes[END_REF]; [START_REF] Tijsseling | Water hammer with fluid-structure interaction in thick-walled pipes[END_REF].

Figure 1 :

 1 Figure 1: (a) Topography of the considered water distribution network. The network is composed of 23 784 nodes, 26 094 pipes for a total pipe distance close to 1 200 km. The positions of high-frequency pressure detectors (from 1 to 18), arranged within the network are illustrated with green dots. (b) Probability Density Functions (PDF) (with the corresponding histogram in zoom's inset) of dimensionless lengths (median 0.012 indicated with vertical dotted lines) normalized by maximal length. (c) same as (b) for dimensionless diameters (median 0.05) normalized with maximal diameter.

Figure 2 :

 2 Figure 2: Illustration of change point detection applied to a real pressure signal detected on the considered WDN (figure 1a)

Figure 3 :

 3 Figure 3: Time-reversal propagation method principle. (a) awaiken detectors d i , i = 1, 2, 3 at first passage time t i , i = 1, 2, 3 (b) Shortest-path between detector d 1 and one node. (c,d,e) Back-propagation cartography for detector d 1 -d 2 -d 3 . (f) For a given arrival time vector, the source node is the one having a back-propagation vector with minimum variance. (g) Depending on the first passage time t i , the resulting time arrival time vector can lead to several possible source nodes.

Figure 4 :

 4 Figure 4: Calibration of the method for the city network and the 18 real sensors (cf. figure 1a). 5% of the network nodes are chosen independently and randomly as sources. (a) network with random sources in red. (b) Calibration histogram and probability density function.

Figure 5 :

 5 Figure 5: Influence of detector number on the selectivity and success rate of the source geolocalisation. For each detector number, 30 random sets are calibrated, using 5% of temptative randomly chosen true sources among all nodes (1189). (a) Exact and one-neighbor detection versus function of detector number. In each case, the average behavior is plotted with bold continuous lines where respective colored zone displays ± the standard deviation around it. (b) Average values of the potentially successful source number versus detector number.

Figure 6 :

 6 Figure 6: Noise influence on the geolocalization (a) accuracy and (b) selectivity. Using the 18 detector locations provided in figure 1a, the calibration is performed by randomly adding a multiplicative noise to each pipe travel time, for each tested random source. Evolutions of (a) success rates and (b) potentially successful source number versus noise amplitude. Dashed lines display linear regressions.

Figure 7 :

 7 Figure 7: Case 1 illustrating real event detection from the field's data recordings. (a) Normalized pressure signals recorded by high-frequency detectors. (b) Geolocalization of the source with three convex hulls associated with three calibration levels :

  isothermal fluid bulk modulus, the water-hammer wave speed is given by c = c -• c p , C s = E/ρ s c p being the pulse wave velocity ratio, whilst D = ρ f /ρ s being the fluid to solid density one.

Table 1 :

 1 Distribution of pipe materials. Others materials are inox, fiber, PE, and PVD.

	Material Quantity (%) Length (%)
	Iron	87.66	89.28
	Steel	1.97	4.04
	PVC	1.01	1.22
	Others	2.15	1.46
	Unknown	7.21	4.00

Table 2 :

 2 Case # % of pot. suc. source number Hull area (%) Hull areas (normalized by the total area of the network) and subgraph pipe lengths (normalized by the total pipe length of the network) versus the node proportion kept in the potentially successful source list for the real events illustrated in figure 7 and 8.

	Pipe length (%)

2016).Since the water-hammer wave speed in WDNs is very fast (between 700 for 350 m/s for the softest plastic pipes to 1200 m/s for cast iron pipes depending on pipes' mechanical properties) having a rather slow attenuation, the wave generated at a given location can propagate over a large portion of the entire urban Fully documented templates are available in the elsarticle package on CTAN.Preprint submitted to Water researchDecember 13, 2022

This result in wave-propagation heterogeneity to be taken care-off from of considering the adapted velocity in

t i , i = 1, 2, 3, results in back-propagating-time vectors at each nodes (each time is colored with the same

Acknowledgments

This work was supported by the collaborative ANRT Grant CIFRE 2019/1453 co-funded by SETOM, a dedicated society of Veolia for the public drinking water service of Toulouse Métropole operating under the brand Eau de Toulouse Métropole. The authors gratefully acknowledge F. Nospelle, from SETOM, for sharing his network knowledge and for his contribution in gathering the on-field informations necessary to confirm the sources of the transient events herein analyszed.

Declaration of Competing Interest

The authors have no competing interests to declare.

Appendix. Local water-hammer wave velocity

Given the Young modulus E of the pipe, the Poisson coefficient ν s , the pipe's density ρ s , the dimensionless thickness α = e/R (e being the pipe thickness and R its radius), the acoustic water velocity c 0 , K the