Ismael Bailleul

James Norris

Diffusion in small time in incomplete sub-Riemannian manifolds

For incomplete sub-Riemannian manifolds, and for an associated second-order hypoelliptic operator, which need not be symmetric, we identify two alternative conditions for the validity of Gaussian-type upper bounds on heat kernels and transition probabilities, with optimal constant in the exponent. Under similar conditions, we obtain the small-time logarithmic asymptotics of the heat kernel, and show concentration of diffusion bridge measures near a path of minimal energy. The first condition requires that we consider points whose distance apart is no greater than the sum of their distances to infinity. The second condition requires only that the operator not be too asymmetric.

Introduction and summary of results

Let M be a connected C ∞ manifold, which is equipped with a C ∞ sub-Riemannian structure X 1 , . . . , X m and a positive C ∞ measure ν. Thus, X 1 , . . . , X m are C ∞ vector fields on M which, taken along with their commutator brackets of all orders, span the tangent space at every point, and ν has a positive C ∞ density with respect to Lebesgue measure in each coordinate chart. Consider the symmetric bilinear form a on T * M given by

a(x) = m =1 X (x) ⊗ X (x).
Let L be a second order differential operator on M with C ∞ coefficients, such that L1 = 0 and L has principal symbol a/2. In each coordinate chart, L takes the form

L = 1 2 d i,j=1 a ij (x) ∂ 2 ∂x i ∂x j + d i=1 b i (x) ∂ ∂x i (1)
on Ω x,y by µ x,y t = µ t,x,y • σ -1 t .

We focus mainly on two problems, each associated with a choice of the endpoints x and y, and with the limit t → 0. The first is to give conditions for the validity of Varadhan's asymptotics for the heat kernel t log p(t, x, y) → -d(x, y) 2 /2

(2

)
where d is the sub-Riemannian distance. The second is to give conditions for the weak limit

µ x,y t → δ γ (3)
where γ is a path of minimal energy in Ω x,y . We wish to understand, in particular, what can be said without symmetry or ellipticity of the operator L, and without compactness or even completeness of the underlying space M . The heat kernel and the bridge measures have a global dependence on L, while the limit objects have a more local character, so the limits depend on some localization of diffusion in small time. We will give two sufficient conditions for this localization, the first generalizing from the Riemannian case a criterion of Hsu [START_REF] Hsu | Heat kernel on noncomplete manifolds[END_REF] and the second requiring a 'sector condition' which ensures that the asymmetry in L is not too strong. We will thus give new conditions for the validity of ( 2) and (3), which do not require completeness, symmetry or ellipticity, nor do they require any condition on the measure ν. In a companion paper [START_REF] Bailleul | Small-time fluctuations for the bridge of a sub-Riemannian diffusion[END_REF], we have further investigated the small-time fluctuations of the diffusion bridge around the minimal path γ, which reveal a Gaussian limit process.

In this section, we state our three asymptotic results. In the next, we discuss related prior work. Later in the paper, we state three further results. The first of these, Proposition 4.1, shows that the dual characterization for complete sub-Riemannian metrics, proved by Jerison and Sanchez-Calle [START_REF] Jerison | Subelliptic, second order differential operators[END_REF], extends to the incomplete case. Then Propositions 5.1 and 5.2 give Gaussian-type upper bounds, for heat kernels and hitting probabilities respectively, from which the asymptotic results are deduced.

Let K be a closed set in M and set U = M \ K. Write p U for the Dirichlet heat kernel of L in U , extended by 0 outside U × U . Define p(t, x, K, y) = p(t, x, y) -p U (t, x, y).

Then

p(t, x, K, y) = p(t, x, y)µ x,y t ({ω ∈ Ω x,y : ω s ∈ K for some s ∈ [0, 1]}) . We call p(t, x, K, y) the heat kernel through K. In the case where U is relatively compact, we write p(t, x, K) for the hitting probability for K, given by 4

p(t, x, K) = P x (T t) = 1 - ¢ U p U (t, x, y)ν(dy) where T = inf{t ∈ [0, ζ) : B t ∈ K}.
Recall that the sub-Riemannian distance is given by

d(x, y) = inf{ I(γ) : γ ∈ Ω x,y }
where I(γ) denotes the energy 5 of γ associated to the bilinear form a. It is known that d defines a metric on M which is compatible with the topology of M . Set

d(x, K) = inf{d(x, z) : z ∈ K} d(x, K, y) = inf{d(x, z) + d(z, y) : z ∈ K}. 4 Note that p(t, x, K) ¢ M p(t, x, K, y)ν(dy) ¢ K p(t, x, y)ν(dy)
and the first inequality is strict if the process explodes, while the second inequality is always strict because the process returns to U with positive probability after hitting K. 5 For an absolutely continuous path γ : [0, 1] → M , the energy I(γ) is given by

I(γ) = inf ¢ 1 0 ξ t , a(γ t )ξ t dt
where the infimum is taken over all measurable paths ξ : [0, 1] → T * M such that ξ t ∈ T * γt M for all t and, for almost all t, γt = a(γ t )ξ t .

If γ is not absolutely continuous or there is no such path ξ, then we set

I(γ) = ∞.
Note that

d(x, K) + d(y, K) d(x, K, y). Define 6 d(x, ∞) = sup{d(x, K) : K closed and M \ K relatively compact}. Theorem 1.1. Suppose that there is a C ∞ 1-form β on M such that Lf = 1 2 div(a∇f ) + a(β, ∇f ) (4)
where the divergence is understood with respect to ν. Then, for all x, y ∈ M and any closed set K in M with M \ K relatively compact, we have

lim sup t→0 t log p(t, x, K) -d(x, K) 2 /2 (5)
and lim sup

t→0 t log p(t, x, K, y) -(d(x, K) + d(y, K)) 2 /2. (6) Moreover, if there is a constant λ ∈ [0, ∞) such that sup x∈M a(β, β)(x) λ 2 (7)
then, for any closed set K in M , lim sup

t→0 t log p(t, x, K, y) -d(x, K, y) 2 /2. ( 8
)
Moreover, all the above upper limits hold uniformly in x and y in compact subsets of M \ ∂K.

The sector condition (7) limits the strength of the asymmetry of L with respect to ν. We will deduce from Theorem 1.1 the small-time logarithmic asymptotics of the heat kernel. Theorem 1.2. Suppose that L has the form (4). Define

S = {(x, y) ∈ M × M : d(x, y) d(x, ∞) + d(y, ∞)}.
Then, as t → 0, uniformly on compacts in S,

t log p(t, x, y) → -d(x, y) 2 /2. ( 9
)
Moreover, if L satisfies (7), then (9) holds uniformly on compacts in M × M .

We will deduce from Theorem 1.1 also the following concentration estimate for the bridge measures µ x,y t on Ω x,y . A path γ ∈ Ω x,y is minimal if I(γ) < ∞ and I(γ) I(ω) for all ω ∈ Ω x,y .

We will say that γ is strongly minimal if, in addition, there exist δ > 0 and a relatively compact open set U in M such that7 I(γ) + δ I(ω) for all ω ∈ Ω x,y which leave U .

(

) 10
Theorem 1.3. Suppose that L has the form (4). Let x, y ∈ M and suppose that there is a unique minimal path γ ∈ Ω x,y . Suppose either that

d(x, y) < d(x, ∞) + d(y, ∞),
or that L satisfies (7) and γ is strongly minimal. Write δ γ for the unit mass at γ. Then µ x,y t → δ γ weakly on Ω x,y as t → 0.

The authors would like to thank Michel Ledoux and Laurent Saloff-Coste for helpful discussions. JN would like to acknowledge the hospitality of Université Paul Sabatier, Toulouse, where this work was completed.

Discussion and review of related work

The small-time logarithmic asymptotics for the heat kernel (2) were proved by Varadhan [START_REF] Varadhan | On the behavior of the fundamental solution of the heat equation with variable coefficients[END_REF] in the case when M = R d and a is uniformly bounded and uniformly positive-definite. Azencott [START_REF] Azencott | Géodésiques et diffusions en temps petit[END_REF] considered the case where a is positive-definite but M is possibly incomplete for the associated metric d. He showed [START_REF] Azencott | Géodésiques et diffusions en temps petit[END_REF]Chapter 8,Proposition 4.4], that the condition

d(x, y) < max{d(x, ∞), d(y, ∞)} (11)
is sufficient for a Gaussian-type upper bound which then implies [START_REF] Bailleul | Large deviation principle for bridges of sub-Riemannian diffusion processes[END_REF]. In particular, completeness is sufficient. He showed also [ fail without a suitable global condition on the operator L. Hsu [START_REF] Hsu | Heat kernel on noncomplete manifolds[END_REF] showed that Azencott's condition [START_REF] Jerison | Subelliptic, second order differential operators[END_REF] for (2) could be relaxed to

d(x, y) d(x, ∞) + d(y, ∞) (12)
and gave an example to show that (2) can fail without this condition. The methods in [START_REF] Azencott | Géodésiques et diffusions en temps petit[END_REF] and [START_REF] Hsu | Heat kernel on noncomplete manifolds[END_REF] work 'outwards' from relatively compact subdomains U in M and make essential use of the following identity, which allows to control p in terms of p U . See [1, Chapter 2, Theorem 4.2]. Let U, V be open sets in M with V compactly contained in U . Then, for x ∈ M and y ∈ V , we have the decomposition

p(t, x, y) = 1 U (x)p U (t, x, y) + ¢ [0,t)×∂V p U (t -s, z, y)µ x (ds, dz) (13)
where

µ x = ∞ n=1 µ n x , µ n x ([0, t] × A) = P x (B Tn ∈ A, T n t)
where we set S 0 = 0 and define recursively for n 1

T n = inf{t S n-1 : B t ∈ V }, S n = inf{t T n : B t ∈ U }.
This can be combined with the estimate

µ x ([0, t] × ∂V ) C(U, V )t, C(U, V ) < ∞
to obtain estimates on p(t, x, y) from estimates on p U (t, x, y). The same identity ( 13) is also used elsewhere to deduce estimates under local hypotheses from estimates requiring global hypotheses. See for example [START_REF] Kusuoka | Applications of the Malliavin calculus[END_REF] on hypoelliptic heat kernels, and [START_REF] Grigor | Localized upper bounds of heat kernels for diffusions via a multiple Dynkin-Hunt formula[END_REF] on Hunt processes. Varadhan's asymptotics (2) were extended to the sub-Riemannian case by Léandre [START_REF] Léandre | Majoration en temps petit de la densité d'une diffusion dégénérée[END_REF][START_REF] Léandre | Minoration en temps petit de la densité d'une diffusion dégénérée[END_REF] under the hypothesis M = R d and X 0 , X 1 , . . . , X m are bounded with bounded derivatives of all orders. [START_REF] Léandre | Minoration en temps petit de la densité d'une diffusion dégénérée[END_REF] Here, X 0 is the vector field on M which appears when we write L in Hörmander's form

L = 1 2 m =1 X 2 + X 0 .
Our Theorem 1.2 extends (2) to a general sub-Riemannian manifold, subject either to Hsu's condition [START_REF] Kusuoka | Applications of the Malliavin calculus[END_REF], understood for the sub-Riemannian metric, or to the sector condition [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF].

A powerful approach to analysis of the heat equation emerged in the work of Grigor'yan [START_REF] Grigor'yan | The heat equation on noncompact Riemannian manifolds[END_REF] and Saloff-Coste [START_REF] Saloff-Coste | A note on Poincaré, Sobolev, and Harnack inequalities[END_REF][START_REF] Saloff-Coste | Parabolic Harnack inequality for divergence-form second-order differential operators[END_REF]. They showed that a local volume-doubling inequality, combined with a local Poincaré inequality, implies a local Sobolev inequality, which then allows to prove regularity properties for solutions of the heat equation by Moser's procedure, and then heat kernel upper bounds by the Davies-Gaffney argument. This was taken up in the general context of Dirichlet forms by Sturm who proved a Gaussian upper bound [START_REF] Sturm | Analysis on local Dirichlet spaces. II. Upper Gaussian estimates for the fundamental solutions of parabolic equations[END_REF]Theorem 2.4] under such local conditions, without completeness and for non-symmetric operators. Moreover, in this bound, the intrinsic metric appears with the correct constant in the exponent, which allows to deduce the correct logarithmic asymptotic upper bound (2). This intrinsic metric corresponds in our context to the dual formulation of the sub-Riemannian metric. Our Gaussian upper bounds can be seen as applications of Sturm's result. For greater transparency, we will re-run part of the argument in our context, rather than embed in the general framework and check the necessary hypotheses. The approach thus adopted no longer relies on working outwards from well-behaved heat kernels using ( 13), but reduces the global aspect to a certain sort of L 2 -estimate for solutions of the heat equation, which requires no completeness in the underlying space. One finds that the sector condition ( 7) is enough to prevent pathologies in the L 2 -estimate, thus dispensing with the need for condition [START_REF] Kusuoka | Applications of the Malliavin calculus[END_REF]. This is a significant extension: for example, ( 7) is satisfied trivially by all symmetric operators Lf = 1 2 div(a∇f ), without any control on the diffusivity a or the symmetrizing measure ν near infinity.

The small-time convergence of bridge measures is known in the case of Brownian motion in a complete Riemannian manifold by a result of Hsu [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF]. For a compact sub-Riemannian manifold, it was shown by Bailleul [START_REF] Bailleul | Large deviation principle for bridges of sub-Riemannian diffusion processes[END_REF]. It is also known under the assumption [START_REF] Léandre | Minoration en temps petit de la densité d'une diffusion dégénérée[END_REF] and subject to the condition that a(x) is positive-definite by work of Inahama [START_REF] Inahama | Large deviation principle of Freidlin-Wentzell type for pinned diffusion processes[END_REF]. While the limit is the expected one, given the well-known small-time large deviations behaviour of diffusions, a statement such as Theorem 1.3 appears new, both for incomplete manifolds and in the non-compact sub-Riemannian case.

We have not attempted to minimize regularity assumptions for coefficients but note that their use for upper bounds is limited to certain basic tools. The analysis [START_REF] Nagel | Balls and metrics defined by vector fields. I. Basic properties[END_REF] of metric balls, in particular the volume-doubling inequality [START_REF] Nagel | Balls and metrics defined by vector fields. I. Basic properties[END_REF], is done for the case where X 1 , . . . , X m are C ∞ . Also the Poincaré inequality [START_REF] Sturm | Analysis on local Dirichlet spaces. II. Upper Gaussian estimates for the fundamental solutions of parabolic equations[END_REF] is proved in [START_REF] Jerison | The Poincaré inequality for vector fields satisfying Hörmander's condition[END_REF] in this framework. These points aside, for upper bounds, the C ∞ assumptions on a, ν and β are used only to imply local boundedness. While the dual characterization of the distance function is unaffected by modification of a on a Lebesgue null set, the definition as an infimum over paths is more fragile, and current proofs that these give the same quantity rely on the continuity of a. In contrast to the Riemannian case [START_REF] Norris | Heat kernel asymptotics and the distance function in Lipschitz Riemannian manifolds[END_REF], for lower bounds in the sub-Riemannian case, in particular for Léandre's argument using Malliavin calculus, current methods demand more regularity.

Review of some analytic prerequisites

The set-up of Section 1 is assumed. Nagel, Stein & Wainger's analysis [START_REF] Nagel | Balls and metrics defined by vector fields. I. Basic properties[END_REF] of the sub-Riemannian distance and of the volume of sub-Riemannian metric balls implies the following statements. There is a covering of M by charts φ : U → R d such that, for some constants α(U ) ∈ (0, 1] and

C(U ) ∈ [1, ∞), for all x, y ∈ U , C -1 |φ(x) -φ(y)| d(x, y) C|φ(x) -φ(y)| α . ( 15
)
Moreover, there is a covering of M by open sets U such that, for some constant C(U ) ∈ (1, ∞), for all x ∈ U and all r ∈ (0, ∞) such that B(x, 2r) ⊆ U , we have the volume-doubling inequality ν(B(x, 2r)) Cν(B(x, r)).

Moreover, in [16, Theorem 1], a uniform local equivalent for ν(B(x, r)) is obtained, which implies that, for all x ∈ M ,

lim r→0 log(ν(B(x, r))) log r = N (x). (17)
Here, N (x) is given by

N (x) = N 1 (x) + 2N 2 (x) + 3N 3 (x) + . . . (18)
where

N 1 (x) + • • • + N k (x)
is the dimension of the space spanned at x by brackets of the vector fields X 1 , . . . , X m of length at most k. While the limit ( 17) is in general not locally uniform, there is also the following uniform asymptotic lower bound on the volume of small balls, for any compact set F in M , lim sup

r→0 sup x∈F log(ν(B(x, r))) log r N (F ) (19)
where

N (F ) = sup x∈F N (x) < ∞.
We recall also the local Poincaré inequality proved by Jerison [START_REF] Jerison | The Poincaré inequality for vector fields satisfying Hörmander's condition[END_REF]. There is a covering of M by open sets U such that, for some constant C(U ) < ∞, for all x ∈ U and all r ∈ (0, ∞) such that B(x, 2r) ⊆ U , for all f ∈ C ∞ c (M ), we have

¢ B(x,r) |f -f B(x,r) | 2 dν Cr 2 ¢ B(x,2r) a(∇f, ∇f )dν ( 20
)
where

f B = ¡ B f dν/ν(B)
is the average value of f on B. As Saloff-Coste claimed [START_REF] Saloff-Coste | Parabolic Harnack inequality for divergence-form second-order differential operators[END_REF]Theorem 7.1], the validity of Moser's argument, given [START_REF] Nagel | Balls and metrics defined by vector fields. I. Basic properties[END_REF] and [START_REF] Sturm | Analysis on local Dirichlet spaces. II. Upper Gaussian estimates for the fundamental solutions of parabolic equations[END_REF], extends with minor modifications to suitable non-symmetric operators. This leads to the following parabolic mean-value inequality. Proposition 3.1. Let L be given as in equation (4) and let U be a relatively compact open set in M . Then there is a constant C(U ) < ∞ with the following property. For any non-negative weak solution u of the equation (∂/∂t)u t = Lu t on (0, ∞) × U , for all x ∈ U , all t ∈ (0, ∞) and all r ∈ (0, ∞) such that B(x, 2r) ⊆ U and r 2 t/2, we have

u t (x) 2 C 2 t t-r 2 2 B(x,r) u 2 s dνds. (21)
Moreover, the same estimate holds if L is replaced by its adjoint L under ν.

For a detailed proof, the reader may check the applicability of the more general results [4, Theorem 1.2] or [START_REF] Lierl | Parabolic Harnack inequality for time-dependent nonsymmetric Dirichlet forms[END_REF]Theorem 4.6].

Dual formulation of the sub-Riemannian distance

In Riemannian geometry, the distance function has a well known dual formulation in terms of functions of sub-unit gradient. Jerison & Sanchez-Calle [START_REF] Jerison | Subelliptic, second order differential operators[END_REF] showed that this dual formulation extends to complete sub-Riemannian manifolds. We now show that such a dual formulation holds without completeness, and for the distances to and through a given closed set. where F is the set of all locally Lipschitz functions w on M such that a(∇w, ∇w) 1 almost everywhere.

Proof. Denote the right hand sides of ( 22) and (23) by δ(x, K, y) and δ(x, K) for now. First we will show that δ(x, K, y) d(x, K, y). Let ω ∈ Ω x,y and suppose that ω is absolutely continuous with driving path ξ and that ω t ∈ K. Let w -, w + ∈ F, with w + = w -on K. It will suffice to consider the case where ω| [0,t] and ω| [t,1] are simple, and then to choose relatively compact charts U 0 and U 1 for M containing ω| [0,t] and ω| [t,1] respectively. Then, given ε > 0, since a is continuous, for i = 1, 2, we can find

C ∞ functions f - i , f + i on U i such that |f ± i (z) -w ± (z)| ε and a(∇f ± i , ∇f ± i )(z) 1 + ε for all z ∈ U i . Then w + (y) -w -(x) = w + (y) -w + (ω t ) + w -(ω t ) -w -(x) f + 1 (y) -f + 1 (ω t ) + f - 0 (ω t ) -f - 0 (x) + 4ε and f + 1 (y) -f + 1 (ω t ) + f - 0 (ω t ) -f - 0 (x) = ¢ t 0 ∇f - 0 (ω s ), ωs ds + ¢ 1 t ∇f + 1 (ω s ), ωs ds = ¢ t 0 ∇f - 0 (ω s ), a(ω s )ξ s ds + ¢ 1 t ∇f + 1 (ω s ), a(ω s )ξ s ds ¢ t 0 a(∇f - 0 , ∇f - 0 )(ω s )ds + ¢ 1 t a(∇f + 1 , ∇f + 1 )(ω s )ds 1/2 ¢ 1 0 a(ξ s , ξ s )ds 1/2
(1 + ε)I(ω).

Hence w + (y) -w -(x) I(ω). On taking the supremum over w ± and the infimum over ω, we deduce that δ(x, K, y) d(x, K, y).

For w ∈ F with w = 0 on K and for y ∈ K, we can take w -= -w and w + = 0 in (22) to see that δ(x, K) δ(x, K, y). Hence, on taking the infimum over y ∈ K in (24), we obtain δ(x, K) d(x, K).

Now we prove the reverse inequalities. Consider a C ∞ symmetric bilinear form ā on T * M such that ā a and ā is everywhere positive-definite. Write Ī for the associated energy function and write d and δ for the distance functions obtained by replacing a by ā in the definitions of d and δ. Set

w + (z) = d(x, K, z), w -(z) = d(x, z), w(x) = d(x, K).
Note that w + = w -and w = 0 on K. Since ā is positive-definite, the functions w -, w + and w are locally Lipschitz, and their weak gradients ∇w ± and ∇w satisfy, almost everywhere,

ā(∇w ± , ∇w ± ) 1, ā(∇w, ∇w) 1. Hence d(x, K, y) = w + (y) -w -(x) δ(x, K, y) δ(x, K, y), d(x, K) = w(x) δ(x, K) δ(x, K).
We will show that, for all ε > 0 and all d * ∈ [1, ∞), we can choose ā so that, for all x, y ∈ M with d(x, y) d * , d(x, y) d(x, y) + ε.

Then, for this choice of ā, we have also, for all closed sets K with d(x, K), d(y, K) d * -1,

d(x, K, y) d(x, K, y) + 2ε, d(x, K) d(x, K) + ε.
Since ε and d * are arbitrary, this completes the proof. The idea in choosing ā is as follows.

While we have no control over the behaviour of a near ∞, neither do we have any constraint on how small we can choose ā -a near ∞. Given ε > 0, this will allow us to choose ā so that, for any path γ ∈ Ω x,y with Ī(γ) < ∞, we can construct another path γ ∈ Ω x,y with I(γ) Ī(γ) + ε.

It will be convenient to fix C ∞ vector fields Y 1 , . . . , Y p on M which span the tangent space at every point, so that

a 0 (x) = p i=1 Y i (x) ⊗ Y i (x)
is a positive-definite symmetric bilinear form on T * M . There exists an exhaustion of M by open sets (U n : n ∈ N), such that U n is compactly contained in U n+1 for all n. Set U 0 = ∅. Let (δ n : n ∈ N) be a sequence of constants, such that δ n ∈ (0, 1] for all n, to be determined. There exists a positive C ∞ function f on M such that f δ n on M \ U n-2 for all n. We take ā = a + f 2 a 0 . Write d 0 and I 0 for the distance and energy functions associated with a + a 0 . Recall that we write d and Ī for the distance and energy functions associated with ā. Then

[0, 1] → R p such that, for almost all t, ωt = m =1 X (ω t ) ḣ t + p i=1 f (ω t )Y i (ω t ) ki t and ¢ 1 0 | ḣt | 2 dt + ¢ 1 0 | kt | 2 dt = Ī(ω).
By reparametrizing ω if necessary, we may assume that

| ḣt | 2 + | kt | 2 = Ī(ω)
for almost all t.

Consider for now the case where ω t ∈ U n+1 \ U n-1 for all t for some n and define a new path

γ by γt = m =1 X (γ t ) ḣ t , γ 0 = x.
Then I(γ) Ī(ω). By Gronwall's lemma, there is a constant A n ∈ [1, ∞), depending only on n and on the open sets (U n : n ∈ N) and the vector fields X 1 , . . . , X m and Y 1 , . . . , Y p , such that

d 0 (γ 1 , y) d * A n δ n provided that d * A n δ n ε n-2 ∧ ε n+1 . ( 25
)
We will ensure that (25) holds, and hence that γ 1 ∈ U n+2 . Then

d(x, y) d(x, γ 1 ) + d(γ 1 , y) Ī(ω) + C n d 0 (γ 1 , y) αn Ī(ω) + C n d * A n δ αn n .
We return to the general case. Then there is an integer k 1 and there is a sequence of times t 0 t 1 . . . t k and there is a sequence of positive integers n 1 , . . . , n k such that t 0 = 0, t k = 1, and |n j+1 -n j | = 1 and ω t j ∈ ∂U n j+1 for j = 1, . . . , k -1, and

ω t ∈ Ūn j +1 \ U n j -1 for all t ∈ [t j-1 , t j ] and all j = 1, . . . , k, and, if k 2, ω t ∈ ∂U n 1 for some t ∈ [t 0 , t 1 ]. Set S n = {t j : j ∈ {1, . . . , k -1} and n j+1 = n}, χ n = |S n |.
Since ω must hit either ∂U n+1 or ∂U n-1 immediately prior to any time in S n , we have

(ε n-1 ∧ ε n )χ n d * .
We have shown that

d(ω t j-1 , ω t j ) (t j -t j-1 ) Ī(ω) + C n j d * A n j δ αn j n j so d(x, y) k j=1 d(ω t j-1 , ω t j ) Ī(ω) + C n 1 d * A n 1 δ αn 1 n 1 + ∞ n=1 C n d * A n χ n δ αn n .
Now we can choose the sequence (δ n : n ∈ N) so that (25) holds and

2 ∞ n=1 C n d * 2 A n δ αn n ε n-1 ∧ ε n ε.
Then, on optimizing over ω, we see that d(x, y) d(x, y) + ε whenever d(x, y) d * , as required.

Gaussian-type upper bounds

Recall from Section 1 the notions of distance and heat kernel through a given closed set K. Proposition 5.1. Let L be given as in equation (4) and suppose that L satisfies [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF]. Then there is a continuous function C : M × M → (0, ∞) such that, for all x, y ∈ M and all t ∈ (0, ∞), for

r = min t d(x, y) , t 4 , d(x, ∞) 4 , d(y, ∞) 4
we have

p(t, x, y) C(x, y) ν(B(x, r)) ν(B(y, r)) exp - d(x, y) 2 2t + λ 2 t 2 . ( 26
)
Moreover, for any closed set K = M \ D in M , there is a continuous function C(., ., K) : D × D → (0, ∞) such that, for all x, y ∈ D and all t ∈ (0, ∞), for

r = min t d(x, K, y) , t 4 , r(x, K) 4 , r(y, K) 4 , r(x,
K) = min{d(x, ∞), d(x, K)} we have p(t, x, K, y) C(x, y, K) ν(B(x, r)) ν(B(y, r)) exp - d(x, K, y) 2 2t + λ 2 t 2 . ( 27
)
The statements above remain true with the constant 4 replaced by 2, by the local volumedoubling inequality. The value 4 will be convenient for the proof.

Proof. We omit the proof of (26), which is a simpler version of the proof of (27). For (27), we will show that the argument used in [17, Theorem 1.2], for the case where a is positivedefinite and β = 0, generalizes to the present context8 . Consider the set M = M -∪ M + , where M ± = K ∪ D ± and D -, D + are disjoint copies of D = M \ K. Write π for the obvious projection M → M . For functions f defined on M , we will write f also for the function f • π on M . Thus we will sometimes consider a as a symmetric bilinear form on T * D ± and β as a 1-form on D ± . Define a measure ν on M by

ν(A) = ν(A ∩ K) + 1 2 ν(π(A ∩ D -)) + 1 2 ν(π(A ∩ D + )). Note that ν = ν • π -1 . Now define p(t, x, y) =      p(t, x, y) + p D (t, x, y), if x, y ∈ D ± , p(t, x, y) -p D (t, x, y), if x ∈ D ± and y ∈ D ∓ , p(t, x, y), if x ∈ K or y ∈ K.
Given bounded measurable functions f -, f + on M with f -= f + on K, write f for the function on M such that f = f ± • π on M ± , and set f = (f -+ f + )/2 and

f D = (f + -f -)/2
¢ D φ D u D t dν = - 1 2 ¢ M a(∇ φ, ∇ū t )dν + ¢ M a( φβ, ∇ū t )dν - 1 2 ¢ D a(∇φ D , ∇u D t )dν + ¢ D a(φ D β, ∇u D t )dν = - 1 2 ¢ M a(∇φ, ∇u t )dν + ¢ M a(φβ, ∇u t )dν. (28)
Let (w -, w + ) be a pair of bounded locally Lipschitz functions on M such that w -= w + on K and a(∇w ± , ∇w ± ) 1 almost everywhere. Define a function w on M by setting w = w ± • π on M ± . Fix θ ∈ (0, ∞) and set ψ = θw. We deduce from (28) by a standard argument that

d dt ¢ M (e -ψ u t ) 2 dν = - ¢ M a(∇(e -2ψ u t ), ∇u t )dν + 2 ¢ M a(βe -2ψ u t , ∇u t )dν = - ¢ M a(∇u t , ∇u t )e -2ψ dν + 2 ¢ M a((β + ∇ψ)u t , ∇u t )e -2ψ dν ¢ M a(β + ∇ψ, β + ∇ψ)(e -ψ u t ) 2 dν ρ ¢ M (e -ψ u t ) 2 dν
Recall that r 2 t/4. For each s ∈ [t -r 2 , t], we can take f + = cp(s, x, K, .)1 B(y,r) , where c is chosen so that ¡ M f 2 dν = 1. For this choice of f + , we have

u s (x -) 2 = ¢ B + p(s, x -, z) 2 ν(dz).
Hence

p(t, x, K, y) 2 C(V ) ν(B(y, r)) 2 t t-r 2 u s (x -) 2 ds C(U )C(V ) ν(B(x, r))ν(B(y, r))
e -2θ(w + (y)-w -(x)-2r)+ρt .

Here, we applied (31) with τ = s, noting that s 3t/4, so r 2 t/4 s/2. We optimize over (w -, w + ) and take θ = d(x, K, y)/t to obtain

p(t, x, K, y) C(U, V, x, y) ν(B(x, r)) ν(B(y, r)) exp - d(x, K, y) 2 2t + λ 2 t 2 where C(U, V, x, y) = e 2+λd(x,K,y)/2 C(U )C(V ).
Finally, since U is locally finite, there is a continuous function C(., ., K) : D × D → (0, ∞) such that C(U, V, x, y) C(x, y, K) for all U, V ∈ U and all x ∈ U and y ∈ V .

Proposition 5.2. Let L be given as in equation (4). Let U be a relatively compact open set in M and set K = M \ U . There is a constant C(U ) < ∞ with the following property. For all x ∈ U and all t ∈ (0, ∞), and for r = t/d(x, K),

p(t, x, K) C ν(B(x, r)) exp - d(x, K) 2 2t . (32)
Proof. We adapt the argument of the proof of Proposition 5.1. Since ν(U ) < ∞ and p(t, x, K) 1, it will suffice to consider the case where d(x, K) 2 2t. We modify the measure ν and the 1-form β on K, if necessary, by multiplication by suitable C ∞ functions, so that ν(K) 1 and a(β, β)(x) λ 2 for all x ∈ M , for some λ < ∞. This does not affect the value of p(t, x, K) for

x ∈ U . Set f = 1 + 1 U + -1 U -and define, for x ∈ M , u t (x) = ¢ M p(t, x, y)f (y)ν(dy).
Then p(t, x, K) = u t (x -) for x ∈ U . Fix a locally Lipschitz function w on M such that w = 0 on K ∪ U + and a(∇w, ∇w)

1 almost everywhere. Then, as we showed at (30), for all For t > 0 sufficiently small, we have r(t, z) = t/d(z, y) for all z ∈ ∂U , and then ν(B(z, r(t, z))) ν(B(y, r(t, z))) t C(U, y)

θ ∈ [0, ∞), ¢ M (e θw u t ) 2 dν e ρt ¢ M f 2 dν = e ρt (2ν(U ) + ν(K)) where µ x ([0, t] × ∂U ) C(U, V )t, C(U, V ) < ∞. For all z ∈ ∂U , p V (t, z, y) C V (z, y) ν(B(z, r(t, z))) ν(B(y, r(t, z))) exp - d(z, y) 2 2t + λ 2 V t 2 where r(t, z) = min t d(z, y) , t 4 , d(z, ∂V ) 4 , d(y, ∂V ) 4 , λ 2 V = sup z∈V a(β, β)(z) < ∞.
N ( Ū )+1
.

Hence, for t > 0 sufficiently small, and all z ∈ ∂U ,

p V (t, z, y) C U,V (y) t N ( Ū )+1 exp - d(y, K) 2 2t + λ 2 V t 2 where C U,V (y) = sup z∈∂U C V (z, y) × C(U, y) N ( Ū )+1 .
This estimate, along with (32), allows us to short-cut some steps in Hsu's argument. On substituting the estimates into (33) and using the elementary [8, Lemma In the following proof, we introduce an auxiliary real Brownian bridge, from 0 to 1 of speed ε. This is known to converge weakly to a uniform drift as ε → 0. So this auxiliary process provides a new coordinate which acts as a surrogate for time, thereby allowing us to lift the small-time estimate for the heat kernel to a weak convergence result for the associated bridge.

Proof of Theorem 1.3. Consider first the case where L satisfies [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF] and γ is strongly minimal. We will show, for all δ > 0, for Γ t (δ) = {ω t : ω ∈ Ω x,y , I(ω) < d(x, y) 2 + δ} and for r = δ 1/4 (d(x, y) 2 + δ) 1/2 that we have lim sup ε→0 ε log µ x,y ε ({ω ∈ Ω x,y : d(ω t , Γ t (δ)) r for some t ∈ [0, 1]}) -δ/2. (38)

Then, since γ is the unique minimal path in Ω x,y and γ is strongly minimal, for all ρ > 0, there exists δ > 0 such that, for all ω ∈ Ω x,y , we have I(ω) d(x, y) 2 + δ whenever d(ω t , γ t ) ρ for some t ∈ [0, 1]. Hence d(z, γ t ) < ρ for all z ∈ Γ t (δ) and all t ∈ [0, 1]. Then it follows from (38) that, as ε → 0, µ x,y ε ({ω ∈ Ω x,y : d(ω t , γ t ) < r + ρ for all t ∈ [0, 1]}) → 1

showing that µ x,y ε → δ γ weakly on Ω x,y . Consider the operator L and measure ν on M = M × R given by L = L + 1 2

∂ ∂τ Then K is closed in M . Write β 0,1 ε for the law on Ω 0,1 (R) of a Brownian bridge from 0 to 1 of speed ε. Then, with obvious notation, p(t, x, ỹ) = p(t, x, y) 1 √ 2π e -1/(2t) , μx,ỹ ε (dω, dτ ) = µ x,y ε (dω)β 0,1 ε (dτ ).

Proposition 4 . 1 .

 41 For all x, y ∈ M and any closed subset K of M , we have d(x, K, y) = sup{w + (y) -w -(x) : w -, w + ∈ F with w + = w -on K} (22) and d(x, K) = sup{w(x) : w ∈ F with w = 0 on K} (23)

d 0 d d .

 d Set ε n = d 0 (∂U n , ∂U n+1 ). By the sub-Riemannian distance estimate[START_REF] Lierl | Parabolic Harnack inequality for time-dependent nonsymmetric Dirichlet forms[END_REF], there are constants α n ∈ (0, 1] and C n < ∞, depending only on n and on the open sets (U n : n ∈ N) and the vector fields X 1 , . . . , X m and Y 1 , . . . , Y p , such that, for all x, y ∈ U n+2 , d(x, y) C n d 0 (x, y) αn .Fix a constant d * ∈ [1, ∞). Fix x, y ∈ M with d(x, y) d * and suppose that ω ∈ Ω x,y satisfies Ī(ω) d * 2 . There exist absolutely continuous paths h : [0, 1] → R m and k :

  Now inf z∈∂U d(z, y) = d(y, K), sup z∈∂U d(z, y) = C(U, y) < ∞ and, for r > 0 sufficiently small inf z∈ Ū ν(B(z, r)) r N ( Ū )+1 .

2 ,|

 2 ν(dx, dτ ) = ν(dx)dτ where τ denotes the coordinate in R.Then Lf = 1 2 div(ã∇f ) + ã( β, ∇f )where div is the divergence associated to ν and whereã(x, τ ) = a(x) + ∂ ∂τ ⊗ ∂ ∂τ , β(x, τ ), v ± ∂ ∂τ = β(x), v , v ∈ T x M.Moreover, ã has a sub-Riemannian structure andã( β, β)(x, τ ) = a(β, β)(x) λ 2Write Ω 0,1 (R) for the set of continuous paths σ : [0, 1] → R such that σ 0 = 0 and σ 1 = 1. Forσ ∈ Ω 0,1 (R), define σt | 2 dt, if σ is absolutely continuous, ∞, otherwise.Set x = (x, 0) and ỹ = (y, 1), and defineK = M \ Ũ , Ũ = {(γ t , σ t ) : (γ, σ) ∈ Γ(δ), t ∈ [0, 1]}where Γ(δ) = (γ, σ) ∈ Ω x,y × Ω 0,1 (R) : I(γ) + I(σ) < d(x, y) 2 + 1 + δ .

  . Let φ -and φ + be C ∞ functions on M , of compact support, with φ -= φ + on K and define φ on M and φ and φ D on M similarly. For t ∈ (0, ∞), define functions u t on M , ūt on M and Then ūt and u D t solve the heat equation with Dirichlet boundary conditions in M and D respectively. It is straightforward to check that u t = u ± t • π on M ± , where u ±

	u D t on D by					¢
						u t (x) =	M p(t, x, y)f (y)ν(dy)
	and			¢			¢
		ūt (x) =		p(t, x, y) f (y)ν(dy), u D t (x) =	p D (t, x, y)f D (y)ν(dy).
					M			M
								t = ūt ± u D t and
	we extend u D t by 0 on K. Hence ¢	¢	¢
						M φu t dν =	M	φū t dν +	D	φ D u D t dν
	and so d dt ¢	M φu t dν =	d dt	¢	M	φū t dν +	d dt

  This implies (34), since p U (t, x, y) p(t, x, y) and ε and η may be chosen arbitrarily small. In the case where L satisfies[START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF], this follows from Theorem 1.1 by taking K = M . On the other hand, given ε > 0 and a compact set F in S, there is a relatively compact open set U in M such that, for K = M \ U and all (x, y) ∈ F ,

	It remains to show the upper bound
		lim sup	t log p(t, x, y) -d(x, y) 2 /2	(35)
		t→0
	locally uniformly in x and y. d(x, y) -ε d(x, K) + d(y, K).
	Now the restriction of L to U satisfies (7), so
	lim sup
				2.1], we conclude as
	claimed that
	lim sup	t log p(t, x, K, y) -(d(x, K) + d(y, K)) 2 /2.
	t→0

t→0 t log p U (t, x, y) -d U (x, y) 2 /2 -d(x, y) 2 /2 (36) uniformly in (x, y) ∈ F , while, by Theorem 1.1, lim sup t→0 t log p(t, x, K, y) -(d(x, K) + d(y, K)) 2 /2

(37)

also uniformly in (x, y) ∈ F . Since p(t, x, y) = p U (t, x, y) + p(t, x, K, y) and ε is arbitrary, (35) follows from (36) and (37).

It is clear that d(., ∞) is either finite or identically infinite. By the sub-Riemannian version of the Hopf-Rinow theorem, the second case occurs if and only if M is complete for the sub-Riemannian metric. Note that the triangle inequality does not apply 'at K' or 'at ∞', and d(x, K) may exceed d(x, ∞) if M \ K is not relatively compact.

When M is complete for the sub-Riemannian distance, all metric balls are relatively compact, so every minimal path is strongly minimal. Also, if there is a unique minimal path γ ∈ Ω x,y , which is strongly minimal, then, by a weak compactness argument, for all relatively compact domains U containing γ, there is a δ > 0 such that (10) holds.

The idea is to combine a standard argument for heat kernel upper bounds with a reflection trick. In terms of Markov processes, we give a random sign to each excursion of the diffusion process into D, viewing it as taking values in D -or D + . Then a generalization of the classical reflection principle for Brownian motion allows to express the density for paths from x to y via K in terms of this enhanced process. In fact the heat kernel p for this process may be written in terms of p and p D , and we find it technically simpler to define p in those terms, rather than set up the enhanced process.

 [START_REF] Bailleul | Small-time fluctuations for the bridge of a sub-Riemannian diffusion[END_REF]Research supported by EPSRC grant EP/103372X/1 1 for some C ∞ functions b i .

where ρ = a(β + ∇ψ, β + ∇ψ) ∞ (λ + θ) 2 .

Then, by Gronwall's inequality, ¢ M (e -ψ u t ) 2 dν e ρt ¢ M (e -ψ f ) 2 dν.

(29) There exists a locally finite cover U of D by sets of the form B(x, r(x, K)/4), where we recall that r(x, K) = min{d(x, ∞), d(x, K)}. For U = B(x, r(x, K)/4) ∈ U, set Ũ = B(x, 7r(x, K)/8). Then Ũ is a relatively compact open subset of D. By the triangle inequality, for all U ∈ U and all x ∈ U , we have B(x, r(x, K)/2) ⊆ Ũ . Fix U, V ∈ U and write C(U ) and C(V ) for the constants appearing in the parabolic mean-value inequality for L on Ũ and for L on Ṽ . Fix x ∈ U , y ∈ V and t ∈ (0, ∞), and recall that we set

Write x -and y + for the unique points in D -and D + respectively such that π(x -) = x and π(y + ) = y. Set

Take f -= 0 and choose f + 0 supported on B(y, r) and such that

Note that w w -(x) + r on B -and w w + (y) -r on B + . Hence we obtain from (29), for all s 0,

s dν e ρs e -2θ(w + (y)-r) .

Since u - t 0 and (∂/∂t)u - t = Lu - t on (0, ∞) × D, by the parabolic mean-value inequality, for all τ ∈ (0, ∞) such that r 2 τ /2,

Set v s (z) = p(s, x, K, z), then v s 0 and (∂/∂s)v s = Lv s on (0, ∞) × D. By the parabolic mean-value inequality again,

where

By the same argument as that leading to (31), there is a constant C(U ) < ∞ with the following property. For all x ∈ U and all t ∈ (0, ∞), for all r ∈ (0, ∞) such that B(x, 2r) ⊆ U and r 2 t/2, we have

where

and, by optimizing over ε, θ and w, using Proposition 4.1, we obtain

6 Proofs of Theorems 1.1, 1.2 and 1.3

Proof of Theorem 1.1. The asymptotic upper bound [START_REF] Hsu | Heat kernel on noncomplete manifolds[END_REF] for the heat kernel through K, under condition [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF], follows directly from the Gaussian upper bound (27) and the asymptotic lower bound [START_REF] Saloff-Coste | Parabolic Harnack inequality for divergence-form second-order differential operators[END_REF] for the volume of small balls, on letting t → ∞. Similarly, the asymptotic upper bound (5) for the hitting probability for K, when M \ K is relatively compact, follows from (32) and [START_REF] Saloff-Coste | Parabolic Harnack inequality for divergence-form second-order differential operators[END_REF]. It remains to show [START_REF] Grigor | Localized upper bounds of heat kernels for diffusions via a multiple Dynkin-Hunt formula[END_REF]. For this, we adapt an argument of Hsu [START_REF] Hsu | Heat kernel on noncomplete manifolds[END_REF] for the Riemannian case. Consider the L-diffusion process (B t :

We use the identity

Note that P x (T t) = p(t, x, K) and the estimate (32) applies. We estimate p(t, z, y) for z ∈ K using [START_REF] Léandre | Majoration en temps petit de la densité d'une diffusion dégénérée[END_REF]. Choose V relatively compact containing the closure of U . Then, for z ∈ ∂U ,

Proof of Theorem 1.2. First we will show the lower bound lim inf t→0 t log p(t, x, y) -d(x, y) 2 /2 (34) locally uniformly in x and y. Given ε > 0, there exists a simple path γ ∈ Ω x,y , with driving path ξ say, such that I(γ) d(x, y)+ε. We can and do parametrize γ so that a(ξ t , ξ t ) = I(γ) for almost all t ∈ [0, 1]. Fix δ > 0 and consider the open set

We can and do choose δ so that U is compactly contained in the domain of a chart. Choose n 1 such that d(x, y) + ε n δ and fix η ∈ (0, δ/4). For k = 0, 1, . . . , n, set t k = k/n and x k = γ t k and suppose that

We can identify the chart with a subset of R d and choose extensions X0 , X1 , . . . , Xm to R d of the restrictions of X 0 , X 1 , . . . , X m to U such that the extended vector fields are all bounded with bounded derivatives of all orders, such that X1 , . . . , Xm is a sub-Riemannian structure on R d , and such that X0 = χX 0 for some C ∞ function χ vanishing outside the chart. Then, by Léandre's lower bound [START_REF] Léandre | Minoration en temps petit de la densité d'une diffusion dégénérée[END_REF]Theorem II.3] in R d , for k = 1, . . . , n, uniformly in y k-1 and y k , lim inf

On the other hand, by Theorem 1.1, for k = 1, . . . , n, uniformly in y k-1 and y k , lim sup

Hence, by our choice of n and η, uniformly in y k-1 and y k , lim inf

Now, by a standard chaining procedure, we obtain, uniformly in y 0 and y n , lim inf

By Theorem 1.1, we have lim sup

where we have used the lower bound from Theorem 1.2. By standard estimates, we also have

Suppose then that ω ∈ Ω x,y and τ ∈ Ω 0,1 (R) satisfy (ω t , τ t ) ∈ Ũ and |τ t -t| < √ δ/2 for all t ∈ [0, 1]. Then, for each t ∈ [0, 1], there exist s ∈ [0, 1] and γ ∈ Ω x,y and σ ∈ Ω 0,1 (R) such that ω t = γ s , τ t = σ s , I(γ) < d(x, y) 2 + δ, I(σ) < 1 + δ.

Then |σ s -s| √ δ/2 so |t -s| √ δ and so

The estimates (39) and ( 40

Set Ω x,y U = {ω ∈ Ω x,y : ω t ∈ U for all t ∈ [0, 1]}. Then γ is the unique minimal path in Ω x,y U , γ is strongly minimal in Ω x,y U , and p(ε, x, y)1 Ω x,y U (ω)µ x,y ε (dω) = p U (ε, x, y)µ x,y,U ε (dω). (43)

Consider the limit ε → 0. Since the restriction of L to U satisfies [START_REF] Hsu | Brownian bridges on Riemannian manifolds[END_REF], by the first part of the proof, we have µ x,y,U ε → δ γ weakly on Ω x,y U . Since p(ε, x, y) = p U (ε, x, y) + p(ε, x, K, y) it follows from (41) and (42) that p U (ε, x, y)/p(ε, x, y) → 1. Hence, on letting ε → 0 in (43), we see that also µ x,y ε → δ γ weakly on Ω x,y .


domain: spi phys

text:
Moussa Sofiane Karoui

Fatima Zohra Benhalouche

Yannick Deville
email: yannick.deville@irap.omp.eu

Khelifa Djerriri
email: kdjerriri@cts.asal.dz

Xavier Briottet
email: xavier.briottet@onera.fr

Arnaud Le Bris
email: arnaud.le-bris@ign.fr

Detection and area estimation for photovoltaic panels in urban hyperspectral remote sensing data by an original nmf-based unmixing method

Keywords: Hyperspectral imaging, hyperspectral unmixing, partial/informed nonnegative matrix factorization, detection and area estimation, photovoltaic panels

INTRODUCTION

The number of operating renewable energy systems has been increasing rapidly in recent years, particularly in urban areas of developed countries [START_REF] Malof | Automatic detection of solar photovoltaic arrays in high resolution aerial imagery[END_REF]. Photovoltaic panels constitute a significant portion of such systems. Indeed, these solar installations can generate greener energy from non-polluting resources, which has led government agencies, electricity grid operators and decision makers to encourage them by funding and tax reduction. Therefore, and to prevent malicious frauds with these alternative energies, many organizations are increasingly interested in detailed information about these solar installations including their localization and energy production [START_REF] Malof | Automatic detection of solar photovoltaic arrays in high resolution aerial imagery[END_REF]. One of the methods for obtaining this detailed information, at least for the localization, on photovoltaic panels consists in making field surveys. However, this approach is costly and time consuming, which leads to the use of other faster and less expensive methods such as the integration of remote sensing data. Recently, some approaches, using remote sensing data, were proposed for the automatic detection of photovoltaic installations and their localization [START_REF] Malof | Automatic detection of solar photovoltaic arrays in high resolution aerial imagery[END_REF], [START_REF] Puttemans | Detection of photovoltaic installations in RGB aerial imaging: A comparative study[END_REF]. These methods use high spatial resolution airborne/spaceborne RGB images, therefore with a limited number of spectral bands. Nevertheless, such data do not allow an effective detection of photovoltaic panels mainly due to their material properties. Indeed, depending on the irradiation configurations, specular reflections can modify the visual properties of these panels, which makes their detection difficult. To overcome the above limitations, high spectral resolution hyperspectral data may be considered. In fact, hyperspectral imaging sensors simultaneously gather images from hundreds of narrow and contiguous bands of the electromagnetic spectrum, which allows precise material recognition [START_REF] Chang | Hyperspectral data processing: Algorithm design and analysis[END_REF]. However, and principally due to the low spatial resolution of these sensors, mixed pixels, which are characterized by mixtures of spectra of more than one pure material (also called endmember), may occur in the acquired images and prevent direct detection of endmembers without advanced processing to unmix these spectra. Linear Spectral Unmixing (LSU) is one of the most used techniques for processing hyperspectral data. It aims at linearly unmixing all spectra into a set of endmember spectra, and a collection of associated abundance fractions [START_REF] Bioucas-Dias | Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches[END_REF].

In this paper, a hyperspectral-unmixing method is proposed to automatically detect photovoltaic panels and to estimate their areas. This approach is based on an original algorithm with known photovoltaic panel spectra, which has relationships with multiplicative Nonnegative Matrix Factorization (NMF) [START_REF] Cichocki | Nonnegative matrix and tensor factorizations. Applications to exploratory multi-way data analysis and blind source separation[END_REF]. Therefore, the proposed method can be considered as a partial/informed NMF technique [START_REF] Limem | Non-negative matrix factorization using weighted beta divergence and equality constraints for industrial source apportionment[END_REF], [START_REF] Tong | Nonnegative-Matrix-Factorization-Based Hyperspectral Unmixing with Partially Known Endmembers[END_REF]. Experiments based on realistic synthetic and real data are conducted to evaluate the performance of the proposed approach and of the standard multiplicative NMF one. This paper is structured in the following manner. The mathematical data model used in the proposed approach is defined in Section 2. The designed algorithm is presented in Section 3. In Section 4, the conducted experiments are described. Section 5 concludes this paper.

MATHEMATICAL DATA MODEL

As mentioned above, in this section the mathematical data model, used in the LSU techniques for unmixing hyperspectral remote sensing images, is described. Each observed nonnegative reflectance pixel spectrum in considered data is supposed to be a linear mixture of the nonnegative reflectance endmember spectra contained in that pixel. As a result, the entire observed nonnegative reflectance hyperspectral image ∈ can be modeled, in matrix form, as [START_REF] Karoui | A New Spatial Sparsity-Based Method for Extracting Endmember Spectra from Hyperspectral Data with Some Pure Pixels[END_REF] ,

(

where each row of the matrix X represents one spectral band of the considered hyperspectral image. The K pixels of this image are restructured as a one-dimensional array. N represents the number of spectral bands of that image. Each column of the matrix ∈ represents one endmember spectrum and each row of the matrix ∈ represents all nonnegative abundance fractions, in all pixels, of one endmember. These coefficients must satisfy the well-known abundance sum-to-one constraint [START_REF] Bioucas-Dias | Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches[END_REF]. L represents the number of endmembers. As suggested above, the proposed approach assumes that some photovoltaic panel spectra are known and that they are similar to those involved in the observed data. Therefore, the matrix A can be expressed as

, ( 2
)
where the matrix ∈ contains, in its first L 1 columns, the known photovoltaic panel spectra. The rest of this first matrix is composed of zeros. The matrix ∈ contains, in its last L 2 columns, the other, unknown, endmember spectra and the rest of this second matrix is composed of zeros. Obviously, L = L 1 + L 2 . Consequently, Equation (1) becomes .

(3)

In the latter equation, the unknown (to be estimated by the proposed algorithm, described hereafter) matrices are A 2 and, in particular, S.

PROPOSED APPROACH

The method proposed to automatically detect photovoltaic panels in urban hyperspectral remote sensing data is based on the NMF-unmixing technique [START_REF] Cichocki | Nonnegative matrix and tensor factorizations. Applications to exploratory multi-way data analysis and blind source separation[END_REF]. The variables involved in this method consist of two nonnegative matrices, and that respectively aim at estimating A 2 and S, such that .

(4)

Thus, the proposed method can be considered as a partial/informed NMF one [START_REF] Limem | Non-negative matrix factorization using weighted beta divergence and equality constraints for industrial source apportionment[END_REF], [START_REF] Tong | Nonnegative-Matrix-Factorization-Based Hyperspectral Unmixing with Partially Known Endmembers[END_REF]. The designed algorithm minimizes the following cost function

1 2 (5)
by using iterative and multiplicative gradient-based update rules. ‖.‖ F represents the Frobenius norm.

To easily obtain the gradient expressions, the cost function J is rewritten as follows

where Tr(.) and (.) T , respectively, denote the matrix trace and the matrix transpose. Considering the properties provided in [START_REF] Petersen | The matrix Cookbook[END_REF], the gradient expressions of J, with respect to the two matrix variables, are in matrix form ,

.

Since the proposed minimization algorithm is a gradientbased one, the following iterative update rule is considered

← ⊙ , (9)
where ⊙ represents the element-wise multiplication. denotes one of the two matrix variables ( or ), and is a learning rate in matrix form. To guarantee the nonnegativity constraint, an iterative and multiplicative rule is derived, from the above gradient-based one, in the following manner: as it appears in [START_REF] Tong | Nonnegative-Matrix-Factorization-Based Hyperspectral Unmixing with Partially Known Endmembers[END_REF] and [START_REF] Karoui | A New Spatial Sparsity-Based Method for Extracting Endmember Spectra from Hyperspectral Data with Some Pure Pixels[END_REF], the gradient expressions of J with respect to can be written as the difference of two nonnegative functions such that . The nonnegative function corresponds to the terms of ( 7) or ( 8) preceded by a plus sign, whereas corresponds to the terms preceded by a minus sign. The nonnegativity constraint can be satisfied by initializing with a nonnegative value and choosing the value of the learning rate matrix according to

⊘ , (10)
where ⊘ stands for element-wise division. Then, the update rule (9) becomes

← ⊙ ⊘ . ( 11
)
Accordingly, the final proposed iterative and multiplicative update rules for the and matrix variables read

← ⊙ ⊘ , (12)
← ⊙ ⊘ . ( 13
)
In the latter expressions, (a very small and positive value) is added to the denominator of each update rule to avoid division by zero. It should here be noted that the above update rules guarantee the nonnegativity constraint (when the considered matrix variables are nonnegatively initialized) and experimentally yield a descent algorithm, which is such that J decreases throughout its optimization.

As standard NMF methods, the designed algorithm is not assured to provide a unique result and its convergence point depends on its initialization. Thus, and to avoid random initialization (from the viewpoint of the proposed algorithm) of the considered variables, initial estimates of unknown endmember spectra (stacked in the last L 2 columns of matrix ) are calculated by the Vertex Component Analysis (VCA) method [START_REF] Bioucas-Dias | Hyperspectral Unmixing Overview: Geometrical, Statistical, and Sparse Regression-Based Approaches[END_REF]. Also, each initial value of the matrix is set to 1/L. An additional constraint to be considered is the abundance sum-to-one property. To this end, the technique described in [START_REF] Heinz | Fully Constrained Least Squares Linear Spectral Mixture Analysis Method for Material Quantification in Hyperspectral Imagery[END_REF] is used. The adaptation of the considered matrix variables is stopped when the number of iterations reaches a predefined maximum value.

EXPERIMENTAL RESULTS

Tests using realistic synthetic and real data are conducted to evaluate the performance of the designed algorithm. Also, the obtained results are compared to those of the standard Lee and Seung multiplicative NMF method [START_REF] Cichocki | Nonnegative matrix and tensor factorizations. Applications to exploratory multi-way data analysis and blind source separation[END_REF].

Tested data

Four spectra of urban environment components (photovoltaic panel, tiles, grass and tree) are used to create a realistic synthetic hyperspectral image according to the considered linear mixing model [START_REF] Malof | Automatic detection of solar photovoltaic arrays in high resolution aerial imagery[END_REF], except that some "random noise" is first added to the photovoltaic panel spectrum in each image pixel, in order to validate the robustness of the tested methods to spectral variability. The above four pure spectra contain 214 samples, ranging from 0.4 to 2.5 µm, and are obtained by ground measurements using a spectrometer. The abundance fractions are created from a land cover classification of a real urban image, with the above-mentioned four materials, by averaging pixel classification values on a non-overlapping sliding window. Two real sub-images from a full hyperspectral scene are also used in the conducted tests. The airborne hyperspectral image of 1.6m spatial resolution and 214 spectral bands covers a part of Toulouse (France) urban area.

In the conducted tests, the above ground-measured photovoltaic panel spectrum is considered as a known spectrum (i.e. L 1 =1) for both realistic synthetic and real data.

Performance evaluation criteria

For the realistic synthetic data, the Normalized Mean Square Error (NMSE), between the original and estimated photovoltaic panel abundance fractions, is used to evaluate the performance of the tested methods. A smaller value of this criterion indicates a better detection /area estimation for photovoltaic panels.

For the real data, and in addition to a visual inspection of the obtained results, areas of photovoltaic panels, computed from the obtained thresholded abundance fractions (to avoid false detections), are compared with those manually digitized on very high spatial resolution (0.1m) orthoimages of the same zones.

Results and discussion

The proposed approach, called "Multi-Part-NMF" for "Multiplicative Partial NMF", and the standard multiplicative NMF (Multi-NMF) technique [START_REF] Cichocki | Nonnegative matrix and tensor factorizations. Applications to exploratory multi-way data analysis and blind source separation[END_REF] are applied to the above-defined data. These two tested methods are initialized with the same matrix values. The next table shows the values of the performance criterion considered for the realistic synthetic data. The table below shows the photovoltaic panel areas obtained from the two real hyperspectral sub-images. The above table shows that the proposed algorithm is able to automatically estimate photovoltaic panel areas much more accurately than the tested literature algorithm.

CONCLUSION

In this investigation, an original unmixing-based method was proposed to automatically detect photovoltaic panels in urban hyperspectral remote sensing data and to estimate their areas. The designed algorithm, which uses iterative and multiplicative update rules, is based on nonnegative matrix factorization (NMF) with known photovoltaic panel spectra. Compared with the standard multiplicative NMF technique, and according to the obtained results, the proposed algorithm proves to be very attractive in the framework of the detection of photovoltaic panels in urban hyperspectral remote sensing data, and can be of great help to electricity grid operators and decision makers. Future extensions of this investigation will particularly consist in improving the initialization stage of the proposed algorithm by using advanced LSU techniques to set initial unknown endmember spectra. Other extensions of this work will consist in performing further tests, by applying the proposed approach to other synthetic/real data.

Figure 1 .

 1 Figure 1. (a) Ortho-image. Estimated photovoltaic panel abundance fraction maps: (b) Multi-Part-NMF, (c) Multi-NMF. Red circles surround photovoltaic panels.

Table I .

 I Photovoltaic panel NMSE values (%). This figure confirms, by visual inspection, that the proposed algorithm is capable to better detect photovoltaic panels than the tested literature technique.

	images.
	Multi-Part-NMF Multi-NMF
	16.62	81.76
	The above table shows that the proposed approach achieves
	much better detection of photovoltaic panels than the
	standard multiplicative NMF algorithm. Indeed, this table
	reveals that the proposed Multi-Part-NMF algorithm is able
	to reach 16.62% NMSE for photovoltaic panels, whereas the
	literature Multi-NMF technique attains 81.76%.
	The next figure shows the ortho-image and estimated (by
	the two tested algorithms) photovoltaic panel abundance
	fraction maps for the two considered real hyperspectral sub-

Table II .

 II Photovoltaic panel areas (m²).

		Manual digitization Multi-Part-NMF Multi-NMF
	Sub-image 1	15	15.27	21.26
	Sub-image 2	12	13.06	1.10

project on "HYperspectral imagery for Environmental urban Planning" (HYEP, no. ANR 14-CE22-0016-01).