Alice Le Brigant
email: alice.lebrigant@math.u-bordeaux.fr

COMPUTING DISTANCES AND GEODESICS BETWEEN MANIFOLD-VALUED CURVES IN THE SRV FRAMEWORK

Computing distances and geodesics between manifold-valued curves in the SRV framework

Alice Le Brigant

Introduction

Computing distances between shapes of open or closed curves is of interest in many applications, from medical imaging to radar signal processing, as soon as one wants to compare, classify or statistically analyze trajectories or contours of objects. While the shape of an organ or the trajectory of an object on a short distance can be modeled by a curve in the plane R 2 or in the ambient space R 3 , some applications provide curves in an intrinsically non-flat space. Simple examples in positive curvature include trajectories on the sphere where points represent positions on the earth, and a negatively-curved space of interest in signal processing is the hyperbolic half plane, which as we will explain later coincides with the statistical manifold of Gaussian densities. We are motivated by the study of curves in the latter for signal processing purposes, however our framework is more general.

Here we consider open oriented curves in a Riemannian manifold M , more precisely the space of smooth immersions c : [0, 1] → M ,

M = Imm([0, 1], M ) = {c ∈ C ∞ ([0, 1], M ), c (t) = 0 ∀t ∈ [0, 1]}.
To compare or average elements of this space, one way to proceed is to equip M with a Riemannian structure, that is to locally define a scalar product G on its tangent space T M. A property that is usually required of this metric is reparametrization invariance, that is that the metric be the same at all points of M representing curves that are identical modulo reparametrization. Two curves are identical modulo reparametrization when they pass through the same points of M but at different speeds. Reparametrizations are represented by increasing diffeomorphisms φ : [0, 1] → [0, 1] (so that they preserve the end points of the curves), and their set is denoted by Diff + ([0, 1]). Elements h, k ∈ T c M of the tangent space in c ∈ M are infinitesimal deformations of c and can be seen as vector fields along the curve c in M (this results from the so called "Exponential law" for smooth functions, see e.g. [START_REF] Kriegl | Aspects of the theory of infinite dimensional manifolds[END_REF], Theorem 5.6.). The Riemannian metric G is reparametrization invariant if the action of Diff + ([0, 1]) is isometric for G (1)

G c•φ (h • φ, k • φ) = G c (h, k),
for any c ∈ M, h, k ∈ T c M and φ ∈ Diff + ([0, 1]). This is often called the equivariance property, and it guarantees that the induced distance between two curves c 0 and c 1 does not change if we reparametrize them by the same diffeomorphism φ

d(c o • φ, c 1 • φ) = d(c 0 , c 1 ).
What's more, a reparametrization invariant metric on the space of curves induces a Riemannian structure on the "shape space", where the space of reparametrizations is quotiented out. A shape can be seen as the equivalence class of all the curves that are identical modulo a change of parameterization, and the shape space as the associated quotient space

S = Imm([0, 1], M )/Diff + ([0, 1]).
While the space of immersions is an open submanifold of the Fréchet manifold C ∞ ([0, 1], M ) (see [START_REF] Michor | Manifolds of differentiable mappings[END_REF], Theorem 10.4.), the shape space is not a manifold and therefore the fiber bundle structure we discuss next is to be understood formally. We get a principal bundle structure π : M → S, which induces a decomposition of the tangent bundle T M = V M ⊕ HM into a vertical subspace V M = ker(T π) consisting of all vectors tangent to the fibers of M over S, and a horizontal subspace HM = (V M) ⊥ G defined as the orthogonal complement of V M according to the metric G that we put on M. If G verifies the equivariance property, then it induces a Riemannian metric Ĝ on the shape space, for which the geodesics are the projected horizontal geodesics of M for G. The geodesic distance d on S between the shapes [c 0 ] and [c 1 ] of two given curves c 0 and c 1 is then given by d Riemannian metrics on the space of curves lying in the Euclidean space R n , and especially closed curves c : S 1 → R n (S 1 is the circle), have been widely studied ( [START_REF] Younes | Computable elastic distances between shapes[END_REF], [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF], [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF], [START_REF] Bauer | Constructing reparametrization invariant metrics on spaces of plane curves[END_REF]). The most natural candidate for a reparametrization invariant metric is the L 2 -metric with integration over arc length d = c (t) dt

([c 0 ], [c 1 ]) = inf d (c 0 , c 1 • φ) | φ ∈ Diff + ([0, 1]) ,
G L 2 c (h, k) = h, k d ,
but Michor and Mumford have shown in [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF] that the induced metric on the shape space always vanishes. This has motivated the study of Sobolev-type metrics ( [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF], [START_REF] Mennucci | Properties of Sobolev-type metrics in the space of curves[END_REF], [START_REF] Younes | A Metric on shape space with explicit geodesics[END_REF], [START_REF] Bauer | Why use Sobolev metrics on the space of curves[END_REF]), where higher order derivatives are introduced. Local existence and uniqueness of geodesics for first and second-order Sobolev metrics on the space of closed plane curves were shown in [START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF] and completion results were given in [START_REF] Mennucci | Properties of Sobolev-type metrics in the space of curves[END_REF].

One first-order Sobolev metric where different weights are given to the tangential and normal parts of the derivative has proved particularly interesting for the applications ( [START_REF] Laga | Landmark-free statistical analysis of the shape of plant leaves[END_REF], [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF])

(2) G c (h, k) = D h N , D k N + 1 4 D h T , D k T d .
In that case c is a curve in R n , •, • denotes the Euclidean metric on R n , D h = h / c is the derivation of h according to arc length, D h T = D h, v v is the projection of D h on the unit speed vector field v = c / c , and

D s h N = D s h - D s h T .
This metric belongs to the class of so-called elastic metrics, defined by

G a,b c (h, k) = a 2 D h N , D k N + b 2 D h T , D k T d ,
for any weights a, b ∈ R + . The parameters a and b respectively control the degree of bending and stretching of the curve. Srivastava et al. introduced in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF] a convenient framework to study the case where a = 1 and b = 1/2 by showing that metric (2) could be obtained by pullback of the L 2 -metric via a simple transformation R called the Square Root Velocity Function (SRVF), which associates to each curve its velocity renormalized by the square root of its norm. A similar idea had been previously introduced in [START_REF] Younes | Computable elastic distances between shapes[END_REF] and then used in [START_REF] Younes | A Metric on shape space with explicit geodesics[END_REF], where a Sobolev metric is also mapped to an L 2 -metric. The general elastic metric G a,b with weights a and b satisfying 4b 2 ≥ a 2 can also be studied using a generalization of the SRVF [START_REF] Bauer | Constructing reparametrization invariant metrics on spaces of plane curves[END_REF].

The SRV framework can be extended to curves in a Lie groupe using translations [START_REF] Celledoni | Shape analysis on Lie groups with applications in computer animation[END_REF], and to curves in a general manifold using parallel transport. For manifoldvalued curves, this can be done in a way that enables to move the computations to the tangent space to the origin of one of the two curves under comparison [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF], [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF], [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] the authors consider the general elastic metric G a,b , but no Riemannian framework is given. In [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], a Riemannian framework is given for the case a = 1, b = 1/2, and the geodesic equations are derived. In this paper we also restrict to this particular choice of coefficients a and b for simplicity, but we propose another generalization of the SRV framework to manifold-valued curves. Instead of encoding the information of each curve within a tangent space at a single point as in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF] using parallel transport, the distance is computed in the manifold itself which enables us to be more directly dependent on its geometry. Intuitively, the data of each curve is no longer concentrated at any one point, and so the energy of the deformation between two curves takes into account the curvature of the manifold along the entire "deformation surface", not just along the path traversed by the starting point of the curve.

In the following section, we introduce our metric as the pullback of a quite natural metric on the tangent bundle T M, and show that it induces a fiber bundle structure over the manifold M seen as the set of starting points of the curves. In section 3, we give the induced geodesic distance and highlight the difference with respect to the distance introduced in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In section 4, we give the geodesic equations associated to our metric and exploit them to build the exponential map. Geodesics of the space of curves can then be computed using geodesic shooting. To this end, we describe the Jacobi fields on M. We test these algorithms on curves lying in the hyperbolic half-plane H, a choice that we motivate in section 5. Finally, in the setting of radar spectral analysis, we model locally stationary radar signals by curves in H and compute their mean.

Extension of the SRV framework to manifold-valued curves

2.1. Our metric on the space of curves. Let c : [0, 1] → M be a curve in M and h, k ∈ T c M two infinitesimal deformations. We consider the following first-order Sobolev metric on M

G c (h, k) = h(0), k(0) + ∇ h N , ∇ k N + 1 4 ∇ h T , ∇ k T d ,
where we integrate according to arc length d = c (t) dt, •, • and ∇ respectively denote the Riemannian metric and the associated Levi-Civita connection of the manifold M , ∇ h = 1 c ∇ c h is the covariant derivative of h according to arc length, and ∇ h T = ∇ h, v v and ∇ h N = ∇ h -∇ h T are its tangential and normal components respectively, with the notation v = c / c . If M is a flat Euclidean space, we obtain the metric (2) studied in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF], with an added term involving the origins. Without this extra term, the bilinear form G is not definite since it vanishes if h or k is covariantly constant along c. Here we show that G can be obtained as the pullback of a very natural metric G on the tangent bundle T M. We consider the square root velocity function (SRVF, introduced in [START_REF] Srivastava | Shape analysis of elastic curves in Euclidean spaces[END_REF]) on the space of curves in M ,

R : M → T M, c → c c ,
where • is the norm associated to the Riemannian metric on M . In order to define G, we introduce the following projections from T T M to T M . Let ξ ∈ T (p,u) T M and t → (x(t), U (t)) be a curve in T M that passes through (p, u) at time 0 at speed ξ. Then we define the vertical and horizontal projections vp (p,u) :

T (p,u) T M → T p M, ξ → ξ V := ∇ x (0) U, hp (p,u) : T (p,u) T M → T p M, ξ → ξ H := x (0).
The horizontal and vertical projections live in the tangent bundle T M and are not to be confused with the horizontal and vertical parts which live in the double tangent bundle T T M and will be denoted by ξ H , ξ V . Furthermore, let us point out that the horizontal projection is simply the differential of the natural projection T M → M , and that according to these definitions, a very natural metric on the tangent bundle T M , the Sasaki metric ( [START_REF] Sasaki | On the differential geometry of tangent bundles of Riemannian manifolds[END_REF], [START_REF] Sasaki | On the differential geometry of tangent bundles of Riemannian manifolds II[END_REF]), can be written where •, • denotes the Riemannian metric on M . Now we can define the metric that we put on T M. Let us consider h ∈ T M and ξ, η ∈ T h T M . We define

g S (p,u) (ξ, η) = ξ H , η H + ξ V , η V ,
Gh (ξ, η) = ξ(0) H , η(0) H + 1 0 ξ(t) V , η(t) V dt,
where ξ(t) H = hp(ξ(t)) ∈ T M and ξ(t) V = vp(ξ(t)) ∈ T M are the horizontal and vertical projections of ξ(t) ∈ T T M for all t. Then we have the following result.

Proposition 1. The metric G on the space of curves M can be obtained as pullback of the metric G by the square root velocity function R, that is

G c (h, k) = GR(c) (T c R(h), T c R(k)) ,
for any curve c ∈ M, and vector fields h, k ∈ T c M along c.

Notation.

Here and in all the paper we will denote by s the parameter of paths in the space of curves M and by t the parameter of a curve in M . For any path of curves s → c(s, •) the corresponding derivatives will be denoted by c s = ∂c/∂s and c t = ∂c/∂t , and we will also use the notations ∇ s = ∇ ∂c/∂s and ∇ t = ∇ ∂c/∂t .

Proof of Proposition 1. For any t ∈ [0, 1], we have

T c R(h)(t) H = h(t) and T c R(h) V = ∇ h R(c)(t).
To prove this proposition, we just need to compute the latter. Let s → c(s, •) be a curve in M such that c(0, •) = c and c s (0, •) = h . Then

∇ h R(c)(t) = 1 c 1/2 ∇ h c + h c -1/2 c = 1 c t 1/2 ∇ s c t + ∂ s c t , c t -1/4 c t = 1 c t 1/2 ∇ t c s - 1 2 c t , c t -5/4 ∇ s c t , c t c t = c 1/2 (∇ h) N + 1 2 ∇ h , c c c c ,
where we used twice the inversion ∇ s c t = ∇ t c s .

2.2. Fiber bundle structures. This choice of metric induces two fiber bundle structures. While the second one is an actual fiber bundle structure between two manifolds, the first structure is over the shape space which, as discussed in the introduction, is not a manifold, and so it should be understood formally. Note that we could obtain a manifold structure by restricting ourselves to so-called "free" immersions, that is elements of M on which the diffeomorphism group acts freely, see [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF].

Principal bundle over the shape space. Just as in the planar case, the fact that the square root velocity function R satisfies

R(c • φ) = φ (R(c) • φ) ,
for all c ∈ M, h, k ∈ T c M and φ ∈ Diff + ([0, 1]), guarantees that the integral part of G is reparametrization invariant. Remembering that the reparametrizations φ ∈ Diff + ([0, 1]) preserve the origins of the curves, we notice that G is constant along the fibers and verifies the equivariance property (1). We then have a formal principal bundle structure over the shape space

π : M = Imm([0, 1], M ) → S = M/Diff + ([0, 1]).
which induces a decomposition T M = V M ⊥ ⊕ HM . There exists a Riemannian metric Ĝ on the shape space S such that π is (formally) a Riemannian submersion from (M, G) to (S, Ĝ)

G c (h H , k H ) = Ĝπ(c) (T c π(h), T c π(k)) ,
where h H and k H are the horizontal parts of h and k, as well as the horizontal lifts of T c π(h) and T c π(k), respectively. This expression does in fact define Ĝ in the sense that it does not depend on the choice of the representatives c, h and k. For more details, see the theory of Riemannian submersions and G-manifolds in [START_REF] Michor | Topics in Differential geometry[END_REF].

Fiber bundle over the starting points. The special role played by the starting point in the metric G induces another fiber bundle structure, where the base space is the manifold M , seen as the set of starting points of the curves, and the fibers are composed of the curves which have the same origin. The projection is then

π ( * ) : M → M, c → c(0).
It induces another decomposition of the tangent bundle in vertical and horizontal bundles

V ( * ) c M = ker T π ( * ) = { h ∈ T c M | h(0) = 0 } , H ( * ) c M = V ( * ) c M ⊥ G .
Proposition 2. We have the usual decomposition

T M = V ( * ) M ⊥ ⊕ H ( * ) M, the horizontal bundle H ( * )
c M consists of parallel vector fields along c, and π ( * ) is a Riemannian submersion for (M, G) and (M, •, • ).

Proof. Let h be a tangent vector. Consider h 0 the parallel vector field along c with initial value h 0 (0) = h(0). It is a horizontal vector, since its vanishing covariant derivative along c assures that for any vertical vector l we have G c (h 0 , l) = 0. The difference h = h -h 0 between those two vectors has initial value 0 and so it is a vertical vector, which gives a decomposition of h into a horizontal vector and a vertical vector. The definition of H ( * ) M as the orthogonal complement of V ( * ) M guaranties that their sum is direct. Now if k is another tangent vector, then the scalar product between their horizontal parts is

G c (h H , k H ) = h H (0), k H (0) c(0) = h(0), k(0) c(0) = T c π ( * ) (h H ), T c π ( * ) (k H ) π ( * ) ,
which proves that π ( * ) is a Riemannian submersion and completes the proof.

Induced distance on the space of curves

Here we give an expression of the geodesic distance induced by the metric G. We show that it can be written similarly to the product distance given in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], with an added curvature term. Let us consider two curves c 0 , c 1 ∈ M, and a path of curves s → c(s, c(s,•)) the image of this path of curves by the SRVF R. Note that q is a vector field along the surface c in M . Let now q be the "raising" of q in the tangent space T c(0,0) M defined by [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF] q(s, t) = P s,0 c(•,0) • P t,0 c(s,•) (q(s, t)) , where we denote by P t1,t2 γ : T γ(t1) M → T γ(t2) M the parallel transport along a curve γ from γ(t 1 ) to γ(t 2 ). Notice that q is a surface in a vector space, as illustrated in Figure 1. Lastly, we introduce a vector field (a, τ ) → ω s,t (a, τ ) in M , which parallel translates q(s, t) along c(s, •) to its origin, then along c(•, 0) and back down again, as shown in Figure 1. More precisely ( 4)

•) linking them in M c(0, t) = c 0 (t), c(1, t) = c 1 (t), for all t ∈ [0, 1]. We denote by q(s, •) = R (
ω s,t (a, τ ) = P 0,τ c(a,•) • P s,a c(•,0) • P t,0 c(s,•) (q(s, t)) ,
for all b, s. That way the quantity ∇ a ω s,t (s, t) measures the holonomy along the rectangle of infinitesimal width shown in Figure 1.

Proposition 3. With the above notations, the geodesic distance induced by the Riemannian metric G between two curves c 0 and c 1 on the space M = Imm([0, 1], M ) of parameterized curves is given by

(5) d(c 0 , c 1 ) = inf c(0,•)=c0,c(1,•)=c1 1 0 c s (s, 0) 2 + 1 0 ∇ s q(s, t) 2 dt ds,
where q = R(c) is the Square Root Velocity representation of the curve c and the norm is the one associated to the Riemannian metric on M . It can also be written as a function of the "raising" q of q in the tangent space T c0(0) M defined by (3), ( 6)

d(c 0 , c 1 ) = inf c(0,•)=c0,c(1,•)=c1 1 0 c s (s, 0) 2 + 1 0 qs (s, t) + Ω(s, t) 2 dt ds,
where Ω is a curvature term measuring the holonomy along a rectangle of infinitesimal width

Ω(s, t) = P s,0 c(•,0) • P t,0 c(s,•) ∇ a ω s,t (s, t) = P s,0 c(•,0) t 0 P τ,0 c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ ,
if R denotes the curvature tensor of the manifold M and ω s,t is defined by (4). Remark 1. The second expression (6) highlights the difference with respect to the distance given in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. In the first term under the square root we can see the velocity vector of the curve c(•, 0) linking the two origins, and in the second the velocity vector of the curve q linking the TSRVF-images of the curves -Transported Square Root Velocity Function, as introduced by Su et al. in [START_REF] Su | Statistical analysis of trajectories on Riemannian manifolds: bird migration, hurricane tracking and video surveillance[END_REF]. If instead we equip the tangent bundle TM with the metric

G h (ξ, ξ) = ξ(0) H 2 + 1 0 ξ(t) V - t 0 P τ,t c R(c , ξ H )P t,τ c q(t) dτ 2 dt,
for h ∈ T M and ξ, η ∈ T h T M, then the curvature term Ω vanishes and the geodesic distance on M becomes ( 7)

d (c 0 , c 1 ) = inf c(0,•)=c0,c(1,•)=c1 1 0 c s (s, 0) 2 + qs (s, •) 2 L 2 ds,
which corresponds exactly to the geodesic distance introduced by Zhang et al. in [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF] on the space

C = ∪ p∈M L 2 ([0, 1], T p M ).
The difference between the two distances ( 5) and ( 7) resides in the curvature term Ω, which measures the holonomy along the rectangle of infinitesimal width shown in Figure 1, and arises from the fact that in the first one, we compute the distance in the manifold, whereas in the second, it is computed in the tangent space to one of the origins of the curves. Therefore, the first one takes more directly into account the "relief" of the manifold between the two curves under comparison. For example, if there is a "bump" between two curves in an otherwise relatively flat space, the second distance ( 7) might not see it, whereas the first one ( 5) should thanks to the curvature term.

Remark 2. Let us briefly consider the flat case : if the manifold M is flat, e.g. M = R n , then the two distances ( 5) and ( 7) coincide. If two curves c 0 and c 1 in R n have the same starting point p, the first summand under the square root vanishes and the distance becomes the L 2 -distance between the two SRV representations q 0 = R(c 0 ) and q 1 = R(c 1 ). If two R n -valued curves differ only by a translation, then the distance is simply the distance between their origins.

Remark 3. Note that this distance is only local in general, that is, only works for curves that are "close enough". Indeed, if we consider two curves c 1 , c 2 in M = R 2 with the same origin, the distance between them is the length of the L 2 geodesic between their SRV representations q 1 and q 2 in C ∞ ([0, 1], R 2 \{0}). If the minimizing geodesic between those two (in C ∞ ([0, 1], R 2 )) passes through 0, then there is no minimizing geodesic between q 1 and q 2 in C ∞ ([0, 1], R 2 \{0}).

Proof of Proposition 3. Since G is defined by pullback of G by the SRVF R, we know that the lengths of c in M and of q = R(c) in TM are equal and so that

d(c 0 , c 1 ) = inf c(0,•)=c0,c(1,•)=c1 1 0 G (q s (s, •), q s (s, •)) ds, with G (q s (s, •), q s (s, •)) = c s (s, 0) 2 + 1 0 ∇ s q(s, t) 2 dt.
To obtain the second expression of this distance we need to express ∇ s q as a function of the derivative qs . Let us fix t ∈ [0, 1], and consider the vector field ν along the surface (s, τ ) → c(s, τ ) that is parallel along all curves c(s, •) and takes value ν(s, t) = q(s, t) in τ = t for any s ∈ [0, 1], that is

ν(s, τ ) = P t,τ c(s,•) (q(s, t)) ,
for all s ∈ [0, 1] and τ ∈ [0, 1]. With this definition we have ∇ s ν(s, t) = ∇ s q(s, t). Since ν(•, 0) : s → P t,0 c(s,•) (q(s, t)) is a vector field along c(•, 0), we can write

∇ s ν(s, 0) = ∇ s P t,0 c(s,•) q(s, t) = P 0,s c(•,0) ∂ ∂s P s,0 c(•,0) •P t,0 c(s,•) (q(s, t)) = P 0,s c(•,0) qs (s, t).
Noticing that we additionally have ∇ τ ν(s, τ ) = 0 for all s, τ ∈ [0, 1], and using

∇ τ ∇ s ν = ∇ s ∇ τ ν + R(c τ , c s )ν, the covariant derivative in τ = t can be written ∇ s ν(s, t) = P 0,t c(s,•) (∇ s ν(s, 0)) + t 0 P τ,t c(s,•) (∇ τ ∇ s ν(s, τ )) dτ = P 0,t c(s,•) • P 0,s c(•,0) (q s (s, t)) + t 0 P τ,t c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ. ( 8
)
Now let us fix s ∈ [0, 1] as well. Notice that the vector field ω s,t defined above as ω s,t (a, τ

) = P 0,τ c(a,•) • P s,a c(•,0) • P t,0 c(s,•) (q(s, t)) verifies ∇ τ ω s,t (s, τ ) = 0 ∀τ ∈ [0, 1], (9)
∇ a ω s,t (a, 0) = 0 ∀a ∈ [0, 1]. ( 10
)
Note that unlike ν, we do not have ∇ a ω s,t (s, t) = ∇ s q(s, t) because ω s,t (a, t) = q(a, t) is only true for a = s. Using Equations ( 9) and ( 10) we get

∇ a ω s,t (s, t) = P 0,t c(s,•) ∇ a ω s,t (s, 0) + t 0 P τ,t c(s,•) ∇ τ ∇ a ω s,t (s, τ ) dτ = t 0 P τ,t c(s,•) ∇ a ∇ τ ω s,t (s, τ ) + R(c τ , c s )ω s,t (s, τ ) dτ, = t 0 P τ,t c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ,
which is the same integral as the one in [START_REF] Bauer | Constructing reparametrization invariant metrics on spaces of plane curves[END_REF]. Finally, since ∇ s q(s, t) = ∇ s ν(s, t) = P s,0 c(•,0) • P t,0 c(s,•) (∇ s ν(s, t)) we obtain

∇ s q(s, t) = qs (s, t) + P s,0 c(•,0) • P t,0 c(s,•) ∇ a ω s,t (s, t) = qs (s, t) + P s,0 c(•,0) t 0 P τ,0 c(s,•) R(c τ , c s )P t,τ c(s,•) q(s, t) dτ
which gives Equation ( 6) and completes the proof.

Computing geodesics

4.1. Geodesic equations on M. To be able to compute the distance given by ( 5) between two curves, we first need to compute the optimal deformation s → c(s, •) from one to the other. That is, we need to characterize the geodesics of M for our metric. In order to do so, taking inspiration from [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF], we use the variational principle. In what follows, we use the lighter notation u(t 1 ) t1,t2 = P t1,t2 c (u(t 1 )) to denote the parallel transport of a vector u(t 1 ) ∈ T c(t1) M along a curve c from c(t 1 ) to c(t 2 ), when there is no ambiguity on the choice of c. We also denote by w T = w, v v the tangential component of any vector field w along a curve c, that is its projection on the unit speed vector field v = c / c . Proposition 4. Let [0, 1] s → c(s, •) ∈ M be a path of curves. It is a geodesic of M if and only if it verifies the following equations ∇ s c s (s, 0) + r(s, 0) = 0, ∀s (11a) ∇ s ∇ s q(s, t) + q(s, t) r(s, t) + r(s, t) T = 0, ∀t, s (11b) where q = c t / c t is the SRV representation of c, the vector field r is given by

r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, and r T = r, v v with v = c t / c t , is the tangential component of r.
Proof. The path c is a geodesic if and only if it is a critical point of the energy functional

E : C ∞ ([0, 1], M) → R + , E(c) = 1 2 1 0 G ∂c ∂s , ∂c ∂s ds.
Let a → ĉ(a, •, •), a ∈ (-, ), be a proper variation of the path s → c(s, •), meaning that it coincides with c in a = 0, and it preserves its end points

ĉ(0, s, t) = c(s, t) ∀s, t, ĉa (a, 0, t) = 0 ∀a, t, ĉa (a, 1, t) = 0 ∀a, t.
Then c is a geodesic of M if and only if d da a=0 E(ĉ(a, •, •)) = 0 for any proper variation ĉ. If we denote by E(a) = E(ĉ(a, •, •)), for a ∈ (-, ), the energy of a proper variation ĉ, then we have

E(a) = 1 2
ĉs (a, s, 0), ĉs (a, s, 0) ds + ∇ s q(a, s, t), ∇ s q(a, s, t) dt ds, where q = ĉt / ĉt is the SRV representation of ĉ. Its derivative is given by E (a) = ∇ a ĉs (a, s, 0), ĉs (a, s, 0) ds + ∇ a ∇ s q(a, s, t), ∇ s q(a, s, t) dt ds.

Considering that the variation preserves the end points, integration by parts gives ∇ a ĉs , ĉs ds = -∇ s ĉs , ĉa ds ∇ s ∇ a q, ∇ s q ds = -∇ s ∇ s q, ∇ a q ds, and so we obtain

E (a) = - ∇ s ĉs , ĉa | t=0 ds + R(c a , c s )q + ∇ s ∇ a q, ∇ s q dt ds = - ∇ s ĉs , ĉa | t=0 ds - R(q, ∇ s q)ĉ s , ĉa + ∇ s ∇ s q, ∇ a q dt ds.
This quantity has to vanish in a = 0 for all proper variations ĉ

∇ s c s | t=0 , ĉa | a=0,t=0 ds + R(q, ∇ s q)c s , ĉa | a=0 + ∇ s ∇ s q, ∇ a q| a=0 dt ds = 0.
We cannot yield any conclusions at this point, because ĉa (0, s, t) and ∇ a q(0, s, t) cannot be chosen independently, since q is not any vector field along ĉ but its image via the Square Root Velocity Function. Computing the covariant derivative of q = ĉt / ĉt 1/2 according to a gives ∇ a q = ĉt -1/2 (∇ a ĉt -1 2 ∇ a ĉt T ), and projecting both sides on v = c t / c t results in ∇ a qT = 1 2 q -1 ∇ a ĉt T . We deduce

∇ a ĉt = q ∇ a q + ∇ a qT ,
and since ∇ t ĉa = ∇ a ĉt , we can express the variation ĉa as follows ĉa (0, s, t) = ĉa (0, s, 0) 0,t + t 0 q(0, s, τ ) ∇ a q(0, s, τ ) + ∇ a qT (0, s, τ ) τ,t dτ.

Inserting this expression in the derivative of the energy we obtain the following, where we omit to write that the variations ĉ and q are always taken in a = 0 for the sake of readability,

1 0 ∇ s c s (s, 0), ĉa (s, 0) ds + 1 0 1 0 R(q, ∇ s q)c s (s, t), ĉa (s, 0) 0,t dt ds + 1 0 1 0 R(q, ∇ s q)c s (s, t), t 0 q(s, τ ) ∇ a q(s, τ ) + ∇ a qT (s, τ ) τ,t dτ dt ds + 1 0 1 0 ∇ s ∇ s q(s, t), ∇ a q(s, t) dt ds = 1 0 ∇ s c s (s, 0) + 1 0 R(q, ∇ s q)c s (s, τ ) τ,0 dτ , ĉa (s, 0) ds + 1 0 1 0 1 t
R(q, ∇ s q)c s (s, τ ) τ,t , q(s, t) (∇ a q(s, t) + ∇ a qT (s, t)) dτ dt ds

+ 1 0 1 0 ∇ s ∇ s q(s, t), ∇ a q(s, t) dt ds = 1 0 ∇ s c s (s, 0) + r(s, 0) , ĉa (s, 0) ds + 1 0 1 0 ∇ s ∇ s q(s, t) + q(s, t) (r(s, t) + r(s, t) T ), ∇ a q(s, t) dt ds =0,
with the previously given definition of r. Since the variations ĉa (0, s, 0) and ∇ a q(0, s, t) can be chosen independently and take any value for all s and all t, we obtain the desired equations.

+ , t) = exp M c(s,t) ( c s (s, t)) ,
for all t ∈ [0, 1] and the second requires the computation of the variation ∇ s c s (s, •)

c s (s + , t) = [c s (s, t) + ∇ s c s (s, t)] s,s+ ,
for all t ∈ [0, 1] where once again, we use the notation w(s) s,s+ = P s,s+ c (w(s)) for the parallel transport of a vector field s → w(s) along a curve s → c(s) in M . If we assume that at time s we have c(s, •) and c s (s, •) at our disposal, then we can estimate c t (s, •) and ∇ t c s (s, •), and deduce q(s, •) = c t (s, •)/ c t as well as [START_REF] Costa | Fisher information distance : a geometrical reading[END_REF] ∇ s q(s,

•) = ∇ s c t |c t | (s, •) - 1 2
∇ s c t , c t |c t | 5/2 c t (s, •),
using the fact that ∇ s c t = ∇ t c s . The variation ∇ s c s (s, •) can then be computed in the following way

(13) ∇ s c s (s, t) = ∇ s c s (s, 0) 0,t + t 0 [∇ s ∇ s c t (s, τ ) + R(c t , c s )c s (s, τ )] τ,t dτ for all t ∈ [0, 1],
where ∇ s c s (s, 0) is given by equation (11a), the second order variation ∇ s ∇ s c t (s, •) is given by

∇ s ∇ s c t = |c t | 1/2 ∇ s ∇ s q + ∇ t c s , c t |c t | 2 ∇ t c s + ∇ s ∇ s q, c t |c t | 3/2 - 3 2
∇ t c s , c t 2 |c t | 4 + |∇ t c s | 2 |c t | 2 c t , (14)
and ∇ s ∇ s q can be computed via equation (11b).

Algorithm 1 (Exponential Map).

Input : c 0 , u. Initialization : Set c(0, t) = c 0 (t) and c s (0, t) = u(t) for all t ∈ [0, 1]. Heredity : For j = 0, . . . , m -1, set s = j with = 1/m and (1) compute for all t c t (s, t) = lim

δ→0 1 δ log M c(s,t) c(s, t + δ), ∇ t c s (s, t) = lim δ→0 1 δ c s (s, t + δ) t+δ,t -c s (s, t) ,
where log M denotes the inverse of the exponential map on M , and compute q(s, t) = 1 √ ct c t (s, t) and ∇ s q(s, t) using equation ( 12).

(2) Compute r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ, and ∇ s ∇ s q(s, t) = -q(s, t) r(s, t) + r(s, t) T , and deduce ∇ s ∇ s c t (s, t) using equation ( 14) for all t ∈ [0, 1].

(3) Initialize ∇ s c s (s, 0) = -r(s, 0) and compute ∇ s c s (s, •) using equation [START_REF] Carmo | Riemannian geometry[END_REF].

(4) Finally, for all t ∈ [0, 1], set

c(s + , t) = exp M c(s,t) ( c s (s, t)) , c s (s + , t) = [ c s (s, t) + ∇ s c s (s, t) ] s,s+ ,
where exp M is the exponential map on the manifold M . Output : c = exp M c0 u. The last step needed to compute the optimal deformation between two curves c 0 and c 1 is to find the appropriate initial speed u, that is the one that will connect c 0 to c 1 . Since we do not have an explicit expression for this appropriate initial speed, we compute it iteratively using geodesic shooting. 4.3. Geodesic shooting and Jacobi fields. The aim of geodesic shooting is to compute the geodesic linking two points p 0 and p 1 of a manifold N , knowing the exponential map exp N . More precisely, the goal is to iteratively find the initial speed u p0,p1 such that exp N p0 (u p0,p1 ) = p 1 . An initial speed vector u ∈ T p0 N is chosen, and is iteratively updated after evaluating the gap between the point p = exp N p0 u obtained by taking the exponential map at point p 0 in u -that is, by "shooting" from p 0 in the direction u -and the target point p 1 . Assuming that the current point p is "not too far" from the target point p 1 , and that there exists a geodesic linking p 0 to p 1 , we can consider that the gap between p and p 1 is the extremity of a Jacobi field J : [0, 1] → N in the sense that it measures the variation between the geodesics s → exp N p0 (su) and s → exp N p0 (su p0,p1 ). Since both geodesics start at p 0 , this Jacobi field has value J(0) = 0 in 0. Then, the current speed vector can be corrected by u ← u + J(0), as shown in Figure 2. Let us briefly explain why. If c(a, s), a ∈ (-, ), s ∈ [0, 1], is a family of geodesics starting from the same point p 0 at different speeds u(a) ∈ T p0 N , i.e. c(a, s) = exp N p0 (su(a)), and J(s) = c a (0, s), s ∈ [0, 1] measures the way that these geodesics spread out, then we have

J(0) = ∂ ∂s s=0 ∂ ∂a a=0 exp N p0 (su(a)) = ∂ ∂a a=0 ∂ ∂s s=0 exp N p0 (su(a)) = u(0).
In the context of geodesic shooting between two curves c 0 and c 1 in M, the speed vector u can be initialized using the L 2 logarithm map, the inverse of the exponential map for the L 2 -metric (these maps are simply obtained by post-composition of mappings with the finite-dimensional maps exp M and log M ). That is, we set u = log L 2 c0 (c 1 ). The L 2 logarithm map also allows us to approximate the gap between the current point and the target point. This amounts to minimizing the functional F (u) = dist L 2 (exp M c0 (u), c 1 ). We summarize as follows. The function φ associates the last value J(1) of a Jacobi field with initial value J(0) = 0 to the initial speed J(0), and can be deduced from Algorithm 3, which describes the function associating J(1) to the initial conditions J(0) and J(0). To find the inverse of this function, we consider the image of a basis of the tangent vector space in which J(0) lives. Now, let us characterize the Jacobi fields of M to obtain the function φ. A Jacobi field is a vector field that describes the way geodesics spread out on a manifold. Consider a → c(a, •, •), a ∈ (-, ), a family of geodesics in M, that is for each a ∈ (-, ), [0, 1] s → c(a, s, •) is a geodesic of M. Then for all a, c(a, •, •) verifies the geodesic equations ∇ s c s (a, s, 0) + r(a, s, 0) = 0, ∀s (15a) ∇ s ∇ s q(a, s, t) + q(a, s, t) r(a, s, t) + r(a, s, t) T = 0, ∀t, s, (15b) where q = c t / c t is the SRV representation of c and r is given by r(a, s, t) = 1 t R(q, ∇ s q)c s (a, s, τ ) τ,t dτ.

Recall that we use the notation w T = w, v(a, s, t) v(a, s, t) with v = c t / c t for the tangential component of a tangent vector w ∈ T c(a,s,t) M . To characterize the way these geodesics spread out, we consider the Jacobi field J : [0, 1] → T M,

J(s, •) = ∂ ∂a a=0 c(a, s, •).
By decomposing ∇ s ∇ s J(s, 0) and ∇ t ∇ s ∇ s J(s, τ ) we can write the second order variation of J as

∇ s ∇ s J(s, t) = ∇ a ∇ s c s + R(c s , J)c s a=0,t=0 0,t + t 0 ∇ s ∇ s ∇ t J + R(c t , c s )∇ s J + ∇ s (R(c t , c s )J) a=0,t=τ τ,t dτ. (16)
The term ∇ s ∇ s ∇ t J can be expressed as a function of ∇ s ∇ s ∇ a q by twice differentiating the equation ∇ a q = c t -1/2 (∇ a c t -1 2 ∇ a c t T ) according to s. This gives

∇ s ∇ s ∇ a q = ∇ s ∇ s ( c t -1/2 )(∇ a c t -1 2 ∇ a c t T ) + 2∇ s ( c t -1/2 ) ∇ s ∇ a c t -1 2 ∇ s ∇ a c t T -1 2 ∇ a c t , ∇ s v v -1 2 ∇ a c t , v ∇ s v + c t -1/2 ∇ s ∇ s ∇ a c t -1 2 ∇ s ∇ s ∇ a c t T -∇ s ∇ a c t , ∇ s v v -∇ s ∇ a c t , v ∇ s v -∇ a c t , ∇ s v ∇ s v -1 2 ∇ a c t , ∇ s ∇ s v -1 2 ∇ a c t , v ∇ s ∇ s v . Since ∇ a c t = ∇ t c a =
∇ a J for a = 0, we know that the term we are looking for is

∇ s ∇ s ∇ t J = ∇ s ∇ s ∇ a c t . Noticing that W = Z -1 2 Z T is equivalent to Z = W +W T we get (17) ∇ s ∇ s ∇ t J = W + W T , W = ∇ s ∇ t J, ∇ s v v + ∇ s ∇ t J, v ∇ s v + ∇ t J, ∇ s v ∇ s v + 1 2 ∇ t J, ∇ s ∇ s v v + 1 2 ∇ t J, v ∇ s ∇ s v + c t ∇ a ∇ s ∇ s q -2∇ s c t -1/2 ∇ s ∇ t J -1 2 ∇ s ∇ t J T -1 2 ∇ t J, ∇ s v v -1 2 ∇ t J, v ∇ s v -∇ s ∇ s c t -1/2 ∇ t J -1 2 ∇ t J T + R(c s , J)∇ s q + ∇ s (R(c s , J)q) . ( 18
)
The terms ∇ a ∇ s c s (0, s, 0) and ∇ a ∇ s ∇ s q(0, s, τ ) for all τ ∈ [0, 1] can be obtained by differentiating the geodesic equations (15a) and (15b) ∇ a ∇ s c s (0, s, 0) + ∇ a r(0, s, 0) = 0, ∀s ∇ a ∇ s ∇ s q(0, s, t) + ∇ a q(0, s, t) r(0, s, t) + r(0, s, t) T + q(0, s, t) ∇ a r(0, s, t) + ∇ a r(0, s, t) T = 0, ∀t, s. The first one gives

(19) ∇ a ∇ s c s (s, 0) = -∇ a r(s, 0),
and for all s and t we get

∇ a ∇ s ∇ s q = - c t ∇ a r + ∇ a r T - 1 c t r,∇ t J v + r, v ∇ t J + 1 2 ∇ t J, v r -3r T . (20)
The only term left to compute is the variation ∇ a r, which is by definition

∇ a r(0, s, t) = 1 t ∇ a V t (0, s, τ ) dτ,
if we define V t for any fixed t by

V t (a, s, τ ) = [R(q, ∇ s q)c s (a, s, τ )] τ,t , τ ∈ [t, 1].
Since the covariant derivative of V t in τ vanishes, we can write for any t ≤ τ ≤ 1

∇ a V t (0, s, τ ) = ∇ a V t (0, s, t) + τ t R(c t , J)| t (V t (0, s, u)) du.
Integrating this equation according to τ from t to 1 we obtain

(21) ∇ a r(0, s, t) = (1 -t)∇ a V t (0, s, t) + R(c t , J)| t 1 t (1 -τ )V t (0, s, τ ) dτ ,
where, since V t (0, s, t) = R(q, ∇ s q)c s (0, s, t), we get for τ = t

(22) ∇ a V t | t = ∇ J R(q, ∇ s q)c s + R(∇ a q, ∇ s q)c s + R(q, ∇ a ∇ s q)c s + R(q, ∇ s q)∇ s J,
with finally

∇ a q = 1 c t ∇ t J -1 2 ∇ t J T , (23a) ∇ a ∇ s q = 1 c t ∇ s ∇ t J -1 2 ∇ s ∇ t J T -1 2 ∇ t J, ∇ s v v -1 2 ∇ t J, v ∇ s v (23b) + ∇ s c t -1/2 ∇ t J -1 2 ∇ t J T + R(J, c s )q.
We can notice that, however complicated, the numbered equations ( 16) to (23b) when put together define a partial differential equation verified by the Jacobi field J. They allow us to iteratively compute J(s + , •) and ∇ s J(s + , •), for a fixed step > 0, knowing J(s, •) and ∇ s J(s, •). Indeed, we can estimate ∇ t J(s, •) since J(s, t) is known for all t, as well as ∇ t ∇ s J(s, •) since ∇ s J(s, t) is known for all t, and finally ∇ s ∇ t J = ∇ t ∇ s J + R(c s , c t )J. Assuming that we are able to compute the covariant derivative ∇ J R of the curvature tensor, for example if we are in a symmetric space (then it is zero), we obtain an algorithm to compute the Jacobi fields in the space of curves. To summarize : Algorithm 3 (Jacobi fields in the space of curves in a symmetric space). Input : c, J 0 and W . Initialization : Set J(0, t) = J 0 (t) and ∇ s J(0, t) = W (t) for t ∈ [0, 1]. Heredity : For j = 0, . . . , m -1, set s = j with = 1/m and

(1) for all t, set

∇ t J(s, t) = lim δ→0 1 δ J(s, t + δ) t+δ,t -J(s, t) , ∇ t ∇ s J(s, t) = lim δ→0 1 δ ∇ s J(s, t + δ) t+δ,t -∇ s J(s, t) , ∇ s ∇ t J(s, t) = ∇ t ∇ s J(s, t) + R(c s , c t )J(s, t).
(

) Compute r(s, t) = 1 t R(q, ∇ s q)c s (s, τ ) τ,t dτ for all t ∈ [0, 1]. ( 2
) Compute ∇ a q(0, s, t), ∇ s ∇ a q(0, s, t) and ∇ a V t (0, s, t) for all t using (23a), (23b) and [START_REF] Michor | Riemannian geometries on spaces of plane curves[END_REF], and deduce ∇ a r(0, s, t) using Equation [START_REF] Michor | Vanishing geodesic distance on spaces of submanifolds and diffeomorphisms[END_REF].

(4) Compute ∇ a ∇ s c s (0, s, 0) and ∇ a ∇ s ∇ s q(0, s, t) for all t using equations [START_REF] Mennucci | Properties of Sobolev-type metrics in the space of curves[END_REF] and (20).

(5) Compute W (s, t) using (18) and ∇ s ∇ s ∇ t J(s, t) = W (s, t) + W T (s, t), and deduce ∇ s ∇ s J(s, t) for all t using Equation ( 16).

(6) Finally, for all t ∈ [0, 1], set

J(s + , t) = [ J(s, t) + ∇ s J(s, t) ] s,s+ , ∇ s J(s + , t) = [ ∇ s J(s, t) + ∇ s ∇ s J(s, t) ] s,s+ .
Output : J(1).

Using a discretization of Algorithms 1, 2 and 3, we are able to compute an approximation of the optimal deformation between two curves, as shown in a toy example in Figure 3. In this simple case we perform geodesic shooting between two geodesics c 0 and c 1 (in red) of the hyperbolic half-plane H. The reasons for our interest in this particular space, as well as the tools needed to work in it, are given in the next section. The first two lines of Figure 3 show the different steps of the first iteration of geodesic shooting, and the last line gives only the last step of the following two iterations. In this simple case, we can see that we converge in only three iterations. Further toy examples are given in Figure 4, where we show the optimal deformations between pairs of geodesics of H (in blue), compared to the L 2 -geodesics (in green). We can see in the first image that our metric has a tendency to "shrink" the curves in the center of the deformation compared to the L 2 -metric. We also show the influence of the orientation of the curves, on which the deformations depend. We do not give the details of the discretization used for these examples or the ones in the following section here. A detailed description of this discrete model will be given in a forthcoming paper.

Example : curves in the hyperbolic half-plane H

In this section we consider the case where the base manifold M is a symmetric manifold of negative curvature, the two-dimensional hyperbolic space H. We first explain why this space can be interesting for applications, namely as it coincides with the statistical manifold of Gaussian densities equipped with the Fisher Information metric. Then we give some basic tools -exponential map, logarithm map, curvature tensor -needed to implement the previous algorithms in H. Finally, we consider a specific application of curve analysis in that space, for the statistical study of locally stationary radar signals. We explain how this framework gives us curves lying in the hyperbolic plane, and we present some simulation results.

5.1.

The hyperbolic half-plane as a statistical manifold. It is possible to adopt a geometrical point of view to solve problems in various fields such as statistical inference, information theory or signal processing [START_REF] Fréchet | Sur l'extension de certaines évaluations statistiques au cas de petits chantillons[END_REF], [START_REF] Burbea | Entropy differential metric, distance and divergence measures in probability spaces: A unified approach[END_REF], [START_REF] Amari | Differential Geometry in Statistical Inference[END_REF]. This framework is given by information geometry. Each element of a parametric family of probability densities {f (•, θ), θ ∈ Θ} can be seen as a point in the manifold of parameters Θ. Intuitively, it is easy to see that the Euclidean metric is not always appropriate to compare two probability distributions in that space. For example, each univariate Gaussian distribution can be identified with its mean and standard deviation (m, σ) in the upper half-plane R × R * + . As explained in [START_REF] Costa | Fisher information distance : a geometrical reading[END_REF], two univariate Gaussian densities N (m 1 , σ 1 ) and N (m 2 , σ 1 ) with different means but the same standard deviation get "closer" to each other as their common standard deviation increases, meaning that intuitively the distance between the points of coordinates (m 1 , σ 1 ) and (m 2 , σ 1 ) in the upper half-plane should be greater than the distance between the points (m 1 , σ 2 ) and (m 2 , σ 2 ) for σ 2 > σ 1 . A more pertinent Riemannian structure on the space of parameters Θ is the one induced by the Fisher information metric, defined in its matrix form as the Fisher information. If the parameter θ ∈ R d is d-dimensional and E denotes the expected value,

g ij (θ) = I(θ) ij = E ∂ ∂θ i ln f (X; θ) ∂ ∂θ j ln f (X; θ) ,
for any 1 ≤ i, j ≤ d and θ = (θ 1 , . . . , θ d ). This metric is chosen, among other reasons, because it has statistical meaning : in parameter estimation, the Fisher information measures the "amount of information" on the parameter θ contained in data sampled from the density f (•, θ); it also gives a fundamental limit to the precision at which one can estimate this θ, in the form of the Cramer-Rao bound.

In the case of univariate Gaussian densities N (m, σ 2 ), Fisher geometry amounts to hyperbolic geometry. More precisely, the space of parameters (m, σ) equipped with the Fisher Information metric is in bijection with the hyperbolic half-plane via the change of variables (m, σ) → ( m √ 2 , σ). Indeed, with this rescaling of the mean, the Fisher Information matrix becomes which defines the Riemannian metric of the well-known hyperbolic half-plane. This is coherent with the example given above, since in the hyperbolic half-plane the distance between the points of coordinates (m 1 , σ) and (m 2 , σ) decreases as σ increases for fixed values of m 1 , m 2 . The differential geometry of Gaussians has proved useful for applications, e.g. in image processing where in the image model, each pixel is represented by a univariate Gaussian distribution [START_REF] Angulo | Morphological processing of univariate Gaussian distribution-valued images based on Poincaré upper-half plane representation[END_REF], and in radar signal processing [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF], [START_REF] Barbaresco | Koszul Information Geometry and Souriau Geometric Temperature/Capacity of Lie Group Thermodynamics[END_REF], [START_REF] Pilté | Tracking quality monitoring based on information geometry and geodesic shooting[END_REF], [START_REF] Brigant | Geometric barycenters of time/Doppler spectra for the recognition of non-stationary targets[END_REF], as we will see in Section 5.3.

g ij (m, σ) = 1 σ 2 0 0 1 σ 2 , - 10 -8 -6 -4

5.2.

Geometry of the hyperbolic half-plane. First, let us give a few tools which are necessary to work in the hyperbolic half-plane representation. Along with the Poincaré disk, the Klein model and others, the hyperbolic half-plane H = {z = x + iy ∈ C, y > 0} is one of the representations of two-dimensional hyperbolic geometry. The Riemannian metric is given by

ds 2 = dx 2 + dy 2 y 2 .
This means that the scalar product between two tangent vectors u

= u 1 + iu 2 and v = v 1 + iv 2 at a point z = x + iy is u, v = u 1 v 1 + u 2 v 2 y 2 .
Using the usual formula (see e.g. [START_REF] Carmo | Riemannian geometry[END_REF]) to compute the Christoffel symbols, we can easily compute the covariant derivative of a vector field

v(t) = v 1 (t) + iv 2 (t) along a curve c(t) = x(t) + iy(t) in H. It is given by ∇ ċ(t) v = X(t) + iY (t) where (24) X = v1 - ẋv 2 + ẏv 1 y , Y = v2 + ẋv 1 -ẏv 2 y .
Let us now remind a well-known expression [START_REF] Carmo | Riemannian geometry[END_REF] for the Riemann curvature tensor in a manifold of constant sectional curvature. Recall that H has constant sectional curvature K = -1.

Proposition 5 (Curvature tensor). Let X, Y, Z be three vector fields on a manifold of constant sectional curvature K. The Riemann curvature tensor can be written

R(X, Y )Z = K ( Y, Z X -X, Z Y ) .
For the algorithms described above, we need to be able to compute the geodesic starting from a point p ∈ H at speed u 0 ∈ T p H -in other words, the exponential map u → exp H p (u) -as well as the geodesic linking two points p and q, with the associated initial vector speed -the inverse q → log H p (q) of the exponential map. The geodesics of the hyperbolic half-plane are vertical segments and half-circles whose origins are on the x-axis, as shown in Figure 5, and they can be obtained as images of the vertical geodesic following the y-axis by a Moebius transformation z → az+b cz+d , with ad -bc = 1. To be complete, we give the proofs of the three following propositions in the appendix.

Proposition 6 (Geodesics of H and logarithm map). Let z 0 = x 0 + iy 0 and z 1 = x 1 + iy 1 be two elements of H.

• If x 0 = x 1 , then the geodesic going from z 0 to z 1 is the segment γ(t) = iy(t) with y(t) = y 0 e t ln y 1 y 0 , and the logarithm map is given by log H z0 (z 1 ) = iy 0 ln y 1 y 0 .

• If x 0 = x 1 , the geodesic is given by γ(t) = x(t) + iy(t) with

x(t) = bd + acȳ(t) 2 d 2 + c 2 ȳ(t) 2 , y(t) = ȳ(t) d 2 + c 2 ȳ(t) 2 , t ∈ [0, 1],
where the coefficients of the Moebius transformation can be deduced from the center x Ω and the radius R of the semi-circle going through z 0 and z 1 : The logarithm map is in turn given by

a = 1 2 xΩ R + 1 , b = x Ω -R, c = 1 2R , d = 1,
log H z0 (z 1 ) = 2cdK ȳ0 2 (d 2 + c 2 ȳ0 2 ) 2 + i K ȳ0 (d 2 -c 2 ȳ0 2 ) (d 2 + c 2 ȳ0 2 ) 2 .
We now give the exponential map in H.

Proposition 7 (Exponential map in H). Let z 0 = x 0 + iy 0 be an element of H and u 0 = ẋ0 + i ẏ0 a tangent vector. Then the exponential map is given by exp H z0 (u 0 ) = γ(1), where

• if ẋ0 = 0, γ(t) = iy 0 e t ẏ0 y 0 ,

• if ẋ0 = 0, γ(t) = x(t) + iy(t) with

x(t) = bd + acȳ(t) 2 d 2 + c 2 ȳ(t) 2 , y(t) = ȳ(t) d 2 + c 2 ȳ(t) 2 , t ∈ [0, 1].
The coefficients a, b, c, d of the Moebius transformation can be computed as previously from the center x Ω = x 0 +y 0 ẏ0 ẋ0 and the radius R = (x 0 -x Ω ) 2 + y 2 0 of the semi-circle of the geodesic, and for all t ∈ [0, 1], ȳ(t) = ȳ0 e t ẏ0 ȳ0 , with ȳ0 = -i az 0 + b cz 0 + d and ẏ0 = ẋ0 (d 2 + c 2 ȳ0 2 ) 2 2cd ȳ0 .

Finally, we give the expression of parallel transport along a geodesic in the hyperbolic plane.

Proposition 8 (Parallel transport in H). Let t → γ(t) be a curve in H with coordinates x(t), y(t), and u 0 ∈ T γ(t0) H a tangent vector. The parallel transport of u 0 along γ from t 0 to t is given by u(t) = y(t) y(t 0 ) cos θ(t 0 , t) sin θ(t 0 , t) -sin θ(t 0 , t) cos θ(t 0 , t) u 0 , where θ(t i , t f ) =

t f ti ẋ(τ )
y(τ ) dτ . If γ is a vertical segment then θ(t i , t f ) = 0, and if it is a portion of a circle, we get

θ(t i , t f ) = 2 (arg(d + icȳ(t f )) -arg(d + icȳ(t i ))) ,
where the coefficients c and d of the Moebius transformation can be computed as explained previously, and γ = iȳ is the pre-image of γ by that transformation. Now that we have these explicit formulas at our disposal, we are able to test the algorithms described above in the simple case where the base manifold M has constant sectional curvature K = -1. Note that computations are further simplified by the existence of a global chart. 5.3. Spectral estimation of locally stationary radar signals. In radar signal processing, given an observation of a signal, it is useful to estimate the spectrum of the underlying process, as it is indicative of its structure. If we are interested in the temporal modulations of that signal, we can estimate several spectra for that same signal and study their evolution in time. The study of these time-frequency spectra, or spectrograms, is at the heart of micro-Doppler analysis. Here we explain how a time series of spectra can be represented as a curve in the Poincaré polydisk.

The data we use for this example is synthetic data generated by a simulator of helicopter signatures. Using this simulator, we obtain a series z = (z 1 , . . . , z N ) ∈ C N of N complex numbers that simulates the reflected signal received by a fixed radar antenna after sending a burst of N pulses in the direction of a fixed helicopter. Given this vector of N observations, the goal is to study the temporal evolutions of the underlying process. To do so, we consider that this process is locally stationary and Gaussian, and we estimate a spectrum for each stationary portion. More precisely, using a gliding window of size n < N to be adjusted, we estimate a high resolution spectrum for each position of the window of size n on the vector of size N . This gives us a time series of N -n + 1 spectra, which we index by 1 ≤ i ≤ N -n + 1.

Let us now explain how we estimate and represent these spectra, each of which corresponds to the observation z (i) = (z i , . . . , z i+n-1 ) of a centered stationary Gaussian time series Z (i) . To overcome the low resolution issues of the classical FFT-based spectral estimation methods, Burg suggested in the 1970s an alternative method based on autoregressive processes [START_REF] Burg | Maximum entropy spectral analysis[END_REF]. Given the partial knowledge of the autocorrelation function of a stationary and Gaussian process Z, Burg showed that the process which maximizes the entropy -that is, which adds the fewest assumptions on the data -is an autoregressive process of the appropriate order. Following this maximum entropy approach, we estimate an autoregressive spectrum for each stationary portion z (i) using the so-called Burg algorithm, see e.g. [START_REF] Arnaudon | Riemannian medians and means with applications to radar signal processing[END_REF]. Burg also showed that the second order statistics of such a process Z can be equivalently represented by its covariance matrix Σ n = E(ZZ * ), a Toeplitz (because of the stationarity) Hermitian Positive Definite (THPD) matrix of size n, or the so-called reflection coefficients (P 0 , µ 1 , . . . , µ n ) of the autoregressive model, as there exists a bijection [START_REF] Trench | An algorithm for the inversion of finite Toeplitz matrices[END_REF], [START_REF] Verblunsky | On positive harmonic functions: a contribution to the algebra of Fourier series[END_REF]

Ψ : T n → R * + × D n-1 , Ψ : Σ n → (P 0 , µ 1 , . . . , µ n-1 ),
between the space T n of THPD matrices of size n and the product space R * + ×D n-1 where these reflection coefficients live. Here D = {z ∈ C | |z| < 1} is the unit disk of the complex plane. This means that each stationary portion z (i) of size n (the size of the gliding window) can be equivalently parameterized by its covariance matrix in T n or by an element of the product space R * + × D n-1 . We choose to work with the latter representation, because we can select a convenient metric on that space. Indeed, the Legendre dual of the Fisher Information metric [START_REF] Barbaresco | Koszul Information Geometry and Souriau Geometric Temperature/Capacity of Lie Group Thermodynamics[END_REF], defined on T n as the hessian of minus the entropy,

Φ(Σ n ) = -ln(det Σ n ) -n ln(2πe), Σ n ∈ T n , has a nice expression in R * + ×D n-1
, in which the Riemannian metric of the Poincaré disk appears [START_REF] Barbaresco | Interactions between symmetric cone and information geometries : Bruhat-Tits and Siegel Spaces Models for High Resolution Autoregressive Doppler Imagery[END_REF]

ds 2 = n dP 0 P 0 2 + n-1 k=1 (n -k) |dµ k | 2 (1 -|µ k | 2 ) 2 ,
with (P 0 , µ 1 , . . . , µ n-1 ) = Ψ(Σ n ). In other words, equipped with the Legendre dual of the Fisher information metric, the space of reflection coefficients becomes the product manifold R * + × D n-1 , where D is the Poincaré disk. This means that each stationary portion Z (i) of our radar signal can be parameterized in the product manifold R * + × D n-1 by a set of coefficients (P 0 (i), µ 1 (i), . . . , µ n-1 (i)), and so the entire locally stationary radar signal is represented by a time series in that space, which corresponds to a set of observations of the "real" evolution (P 0 (t), µ 1 (t), . . . , µ n-1 (t)) of the locally stationary signal.

With this choice of representation, comparing two vectors of radar observations can be carried out by computing the distance between the two corresponding curves in the product manifold R * + × D, which, thanks to the product metric, is the same as comparing their components separately -that is, pairwise comparing the evolutions of each reflection coefficient µ k (t) in the Poincaré disk. More generally, the representation of a vector of radar observation in a Riemannian manifold enables to do basic statistics on these objects, such as defining the mean, median and variance of a set, or performing classification. This can naturally be useful in target detection as well as target recognition. Here we use the algorithms presented in the previous sections to compute the Fréchet mean μ(t) of p curves µ (1) (t), . . . , µ (p) (t). The Fréchet mean, also called intrinsic mean, is defined by

μ = argmin µ∈M p j=1 d(µ, µ (j) ) 2 ,
if M is the space of curves in D and d the distance on M. Since it is defined as a minimizer of a functional, this intrinsic mean can be found by a gradient descent type procedure, summarized as follows.

Algorithm 4 (Mean of a set of curves). Input : (µ (j) (t), 1 ≤ j ≤ p). Initialize μ. Repeat until convergence : (1) For j = 1, . . . , p, compute the geodesic γ (j) (s) linking μ to µ (j) using geodesic shooting (Algorithm 2) and its initial tangent vector γ Using our radar data, we compute the Fréchet mean for 4 sets of 11 curves tracing the evolutions of one reflection coefficient of 11 signals -we choose to represent only one of the coefficients µ k , 1 ≤ k ≤ n -1. These signals are generated using the helicopter signature simulator, and correspond to the observation at 4 different times of 11 helicopters which differ only in their rotor rotation speeds. We consider small variations (less than 1%) around the mean value of 390 RPM (rotations per minute), and show the obtained curves in Figure 6. Theses curves are shown in the hyperbolic half-plane representation, which is equivalent to the Poincaré disk in terms of geometry. In each case, the red extremity of the colormap corresponds to the helicopter with the highest rotation speed, the blue extremity to the lowest rotation speed, and the mean curve is shown in black. This can be used to construct a "reference signature" for a given type of helicopter, for target recognition purposes.

Conclusion

We have studied a first-order Sobolev metric G on the space of manifold-valued curves and its induced geometry. The metric G can be obtained as the pullback of a natural metric on the tangent bundle TM by the square root velocity function, and as such it is reparametrization invariant. The special role that G gives to the starting points of the curves induces a fiber bundle structure over the manifold M seen as the set of starting points of the curves, for which the projection is a Riemannian submersion. The geodesic distance induced by G takes into account the distance between the origins of the curve in M and the L 2 -distance between the SRV representations, without parallel transporting the computations to a unique tangent plane as in [START_REF] Brigant | Reparameterization invariant distance on the space of curves in the hyperbolic plane[END_REF] and [START_REF] Zhang | Video-based action recognition using rate-invariant analysis of covariance trajectories[END_REF]. This should allow us to take into account a greater amount of information on the geometry of the manifold M . Using the pullback form of G, explicit equations can be obtained for the geodesics, as well as for Jacobi fields, which allow us to construct the optimal deformation between two curves by geodesic shooting. Once we can compute geodesics in the space of curves, we can also compute the mean of a set of curves and conceivably more. We considered the case where the base manifold M is the hyperbolic half-plane, whose geometry coincides with the Fisher geometry of gaussian densities, and tested the algorithms on simulated radar data for the spectral analysis of locally stationary gaussian radar signals. Future work will include applications on the sphere for the statistical analysis of large trajectories. and z 1 on iȳ 1 , then the segment corresponding to the portion of γ linking z 0 to z 1 is γ(t) = ȳ0 e t ln ȳ1 ȳ0 . Taking the derivative of (25) in t = 0 gives the logarithm map.

Proof of Proposition 7. The exponential map uses the same equations as the logarithm map with the difference that u 0 is known instead of z 1 . The proof is very similar to the the proof of Proposition 6 and is not detailed here.

Proof of Proposition 8. Parallel transporting a vector u 0 ∈ T γ(t0) H along a curve γ from t 0 to t is gives a vector field u satisfying ∇ γ(t) u = 0 and u(t 0 ) = u 0 . Using Equation [START_REF] Michor | Topics in Differential geometry[END_REF], this can be rewritten u = Au where

A = 1 y ẏ ẋ -ẋ ẏ .
A is of the form aI + bK where I is the identity matrix and K = 0 1 -1 0 and so the solution is u(t) = u(t 0 ) exp The matrix B(t) is diagonalizable and therefore its exponential can be written exp B(t) = e a(t0,t) cos θ(t 0 , t) sin θ(t 0 , t)

-sin θ(t 0 , t) cos θ(t 0 , t) ,

where a(t 0 , t) = ln y(t) y(t0) and θ(t 0 , t) = t t0 ẋ(τ ) y(τ ) dτ . This gives us the desired formula.

  and d verifies the stronger propertyd(c 0 • φ, c 1 • ψ) = d(c 0 , c 1 ),for any reparametrizations φ, ψ ∈ Diff + ([0, 1]). This motivates the choice of a reparametrization invariant metric on M.

Figure 1 .

 1 Figure 1. Illustration of the distance between two curves c 0 and c 1 in the space of curves M.

4. 2 .

 2 Exponential map. Now that we have the geodesic equations, we are able to describe an algorithm which allows us to compute the geodesic s → c(s, •) starting from a point c ∈ M at speed u ∈ T c M. This amounts to finding the optimal deformation of the curve c in the direction of the vector field u according to our metric. We initialize this path s → c(s, •) by setting c(0, •) = c and c s (0, •) = u, and we propagate it using iterations of fixed step > 0. The aim is, given c(s, •) and c s (s, •), to deduce c(s + , •) and c s (s + , •). The first is obtained by following the exponential map on the manifold M c(s

Figure 2 .

 2 Figure 2. Geodesic shooting in the space of curves M

Figure 3 .

 3 Figure 3. Geodesic shooting between two geodesics of the upper half-plane H. The first six subfigures constitute the steps of the first iteration of the algorithm, and the last two show only the result of the two following iterations.

Figure 4 .

 4 Figure 4. Optimal deformations between pairs of geodesics (in black) of the upper half-plane H, for our metric (in blue) and for the L 2 -metric (in green). The orientation of the right-hand curve is inverted in the second image compared to the first, and in the fourth compared to the third.

Figure 5 .

 5 Figure 5. Geodesics of the hyperbolic half-plane.

  and for all t ∈ [0, 1], ȳ(t) = ȳ0 e Kt , with K = ln ȳ1 ȳ0 , ȳ0 = -i az 0 + b cz 0 + d and ȳ1 = -i az 1 + b cz 1 + d .

Figure 6 .

 6 Figure 6. Computation of the mean curve (in black) for 4 sets of 11 curves in the hyperbolic half-plane, constructed from simulated helicopter radar data.

2 )

 2 Update the mean μ in the direction of the sum of the initial speed vectors μ ← exp M μ

  t t0 A(τ )dτ , that is, u(t) = u(t 0 ) exp B(t) with B(t) = ln y

.

Acknowledgments

This research was supported by Thales Air Systems and the french MoD DGA.

Appendix

Proof of Proposition 6. The geodesic γ(t) = x(t)+iy(t) linking two points vertically aligned z 0 = x 0 + iy 0 and z 1 = x 0 + iy 1 is a vertical segment γ(t) = iy(t). It verifies the geodesic equation ∇ γ(t) γ(t) = 0. Using the expression [START_REF] Michor | Topics in Differential geometry[END_REF] of the covariant derivative of a vector field in H, this gives the equation ÿy = ẏ2 , which can be rewritten as ÿ ẏ = ẏ y . Integrating twice, we find that y(t) = y 0 e t ln y 1 y 0 . Now if z 1 = x 1 + iy 1 with x 1 = x 0 , the geodesic γ is the image by a Moebius transformation a → az+b cz+d (with ad -bc = 1) of a vertical line γ(t) = iȳ(t), which gives [START_REF] Pilté | Tracking quality monitoring based on information geometry and geodesic shooting[END_REF] x

We know that γ describes a semi-circle Ω whose origin x Ω is on the x-axis, and that one end of the vertical line γ is sent to the point a/c and the other to the point b/d. This implies that the center of the semi-circle is half-way between the two x Ω = ad+bc 2cd , and that the radius is R = 1 2cd . These two equations as well as the condition ad -bc = 1 gives a system of equations for the coefficients a, b, c and d, which, if we choose to set d = 1, yields the desired expressions. If the extremity z 0 is sent by the inverse of the obtained Moebius transformation on iȳ 0 ,


Ricardo Vilaça
email: rmvilaca@di.uminho.pt

Francisco Cruz
email: fmcruz@di.uminho.pt

José Pereira

Rui Oliveira

An Effective Scalable SQL Engine for NoSQL Databases

Keywords: SQL, NoSQL, Cloud Computing, Middleware Relational CF 999999999 666666666 Relational Column Family

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

With cloud-based databases as part of platform-as-a-service offerings, such as Google's BigTable [START_REF] Chang | Bigtable: a distributed storage system for structured data[END_REF], Amazon's DynamoDB [START_REF] Decandia | Dynamo: Amazon's highly available key-value store[END_REF], and Yahoo!'s PNUTS [START_REF] Cooper | PNUTS: Yahoo!'s hosted data serving platform[END_REF], HBase [START_REF] George | HBase: The Definitive Guide[END_REF] and Cassandra [START_REF] Lakshman | Cassandra -A Decentralized Structured Storage System[END_REF], NoSQL databases become attractive for a larger and more diverse set of applications.

However, their lack of SQL support represents a major hurdle for a wider adoption. This arises at different levels due to the prevalence of SQL as the standard, widely mastered and efficient query language for databases. Most web scale applications are purposely kept SQL-based for their core data management [START_REF] Rys | Scalable SQL[END_REF]. Any application at least two, three years old is directly or indirectly (e.g. on top of an object-relational mapping) based on an SQL interface and its migration is usually not straightforward. A large number of tools and middleware coupled to SQL have been developed and matured over the years and are currently at the basis of most application development frameworks.

It is therefore not surprising that SQL compliance has been one of the most requested additions to the Google App Engine platform [15]. This state of affairs has sparked a number of proposals for middleware that exposes higher-level query interfaces on top of the barebones key-value primitives of NoSQL databases. Many of these aim at approximating the traditional SQL abstraction, ranging from shallow SQL-like syntax for simple key-value queries (e.g. CQL for Cassandra, Phoenix [START_REF]Phoenix: A SQL layer over HBase[END_REF], PIQL [START_REF] Armbrust | PIQL: success-tolerant query processing in the cloud[END_REF]) to the translation of analytical queries into map-reduce jobs [START_REF]Google[END_REF][2][START_REF] Abouzeid | HadoopDB: an architectural hybrid of MapReduce and DBMS technologies for analytical workloads[END_REF][START_REF] Lin | Tenzing: A SQL implementation on the MapReduce framework[END_REF]. However, due to the complexity of SQL existing solutions are limited to a subset of the language, thus not allowing to leverage exiting SQL applications and tools.

In this paper, we present a distributed query engine (DQE) for running SQL queries on top of a NoSQL database, while preserving its scalability and schema flexibility. The DQE allows to combine the expressiveness and performance of a Relational Database Management System (RDBMS) with the scalability and schema flexibility of a NoSQL database. DQE is a transaction-less database [START_REF] Harizopoulos | OLTP through the looking glass, and what we found there[END_REF], preserving the isolation semantics provided by the underlying NoSQL database.

Enabling scalable SQL processing atop of a NoSQL database poses several architectural challenges to the query engine [START_REF] Stonebraker | The end of an architectural era: (it's time for a complete rewrite)[END_REF][START_REF] Stonebraker | 10 rules for scalable performance in 'simple operation' datastores[END_REF]. On one hand, traditional RDBMS architectures include several legacy components such as on-disk data structures, log-based recovery, and buffer management, that were developed years ago but are not suited to modern hardware. Those components impose an huge overhead to transaction processing [START_REF] Harizopoulos | OLTP through the looking glass, and what we found there[END_REF] limiting its scalability. On the other, large scale NoSQL databases have a simple data model, using a simple key-value store or at most variants of the entity-attribute-value (EAV) model [START_REF] Nadkarni | Data Extraction and Ad Hoc Query of an Entity-Attribute-Value Database[END_REF] which strongly limit the expressiveness of data representation. Therefore, a first challenge consists in addressing the impedance mismatches [START_REF] Meijer | A co-relational model of data for large shared data banks[END_REF] while supporting SQL queries. Moreover, it implies mapping relational tables and indexes to the NoSQL database tables in such way that the processing capabilities of the database are fully exploited. The other major challenge regards the basic key-value store interface of NoSQL databases which only allows applications to insert, query, and remove individual tuples or, at most, issue range queries based on the primary key of the tuple. These range queries allow for fast iteration over ranges of rows and also allow to limit the number and what columns are returned. However, NoSQL databases do not support partial key scans, but SQL index scans must perform equality and range queries on all or part of the fields of the index.

Contributions. This paper makes the following three contributions. First, we propose an architecture that allows to combine the expressiveness and performance of RDBMS with the scalability of a NoSQL database. The resulting query processing component is stateless regarding application data and can thus be seamlessly replicated and embedded in the client application. Then, we show how the basic key-value operations and data models can be mapped to scan operators within a traditional (SQL enabled) RDBMS. And finally, we describe a complete implementation of DQE with full SQL support, using Apache Derby's [START_REF] Foundation | Apache derby[END_REF] query engine and HBase as the NoSQL database.

Roadmap. The remainder of the paper is structured as follows. Section 2 introduces the proposed architecture of the DQE. Section 3 describes how it is implemented using Apache Derby components and HBase. Section 4 presents the experimental evaluation. Section 5 compares our approach to other proposals for query processing on a NoSQL database, and Section 6 concludes the paper. The proposed architecture is shown in Figure 1(c), in the context of a scalable NoSQL and a traditional RDBMS. A major motivation for a NoSQL database is scalability. As depicted in Figure 1(a), a typical NoSQL database builds on a distributed setting with multiple nodes of commodity hardware. By adding more nodes to the system we can increase both the overall performance and capacity of the system and also its resilience through data replication. By allowing clients to directly contact multiple fragments and replicas, the system can scale also in terms of clients connected. To make this possible, NoSQL databases provide a simple data model as well as basic querying and searching capabilities, that allow applications to insert, query, and remove individual items or, at most, issue range queries based on the primary key of the item [START_REF] Vilaça | On the expressiveness and trade-offs of large scale tuple stores[END_REF].

Architecture

In sharp contrast, a RDBMS is organized as tables (also called relations). We seldom are concerned with the storage structure but instead express queries in a high-level language, invariably SQL. SQL allows applications to realize complex operations and processing capabilities, such as filtering, joining, grouping, ordering and counting.

Our current proposal builds on rewriting the internal architecture of a typical RDBMS, by reusing some of its components and adding new components atop of a NoSQL database. To understand how components can be reused we examine the internals of a RDBMS roughly splitting it in a query processor block and a storage manager block (Figure 1(b)).

The query processor is responsible for offering an SQL based API to applications, and for translating the application queries, through compilation and execution, to the underlying storage manager.

The architecture proposed in Figure 1(c) reuses a number of components from the SQL query processor (shown in light gray). In detail, these are: the JDBC driver and client connection handler; the compiler and the optimizer, and a set of generic relational operator implementations. These components can be shielded from changes, as they depend only on components that are reimplemented (shown in medium gray) providing the same interfaces as those that, in the RDBMS, embody the centralized storage functionality (shown in dark gray) and that is removed from our architecture. The components to be reimplemented are the following:

-A mapping from the relational model to the data model of a NoSQL database.

This includes: atomic data types and their representation, representation of rows and tables, and representation of indexes. -A mapping from the relational schema to that of the NoSQL database, which allows data to be interpreted as relational tables (see Section 3); -Implementation of sequential and index scan operators. This includes: matching the interface and data representation of the database and leveraging the indexing and filtering capabilities in the NoSQL database to minimize data network traffic;

The proposed architecture has the key advantage of being stateless regarding application data. Data manipulation language (DML) statements (SELECT, INSERT, UPDATE and DELETE) can be executed without any coordination among different DQE instances. As a result, the system should retain the scaleout capabilities of the supporting NoSQL database.

In addition, this architecture also offers the possibility to take advantage of the flexible schema exposed by the underlying NoSQL database. That is, each application applies its own view of the schema over the NoSQL database.

Implementation

The current DQE prototype was built by reusing Apache Derby components and uses HBase as the NoSQL database. In the following, we start with an overview of both systems.

HBase is a key-value based distributed data storage system based on Bigtable [START_REF] Chang | Bigtable: a distributed storage system for structured data[END_REF]. In HBase, data is stored in the form of HBase tables (HTable) that are multidimensional sorted maps. The index of the map is the row's key, column's name, and a timestamp. Columns are grouped into column families. Column families must be created before data can be stored under any column key in that family.

Data is maintained in lexicographic order by row key. Finally, each column can have multiple versions of the same data indexed by their timestamp.

A read or write operation is performed on a row using the row-key and one or more column-keys. Update operations on a single row are atomic, i.e. concurrent writes on a single row are serialized. Any update performed is immediately visible to any subsequent reads. HBase exports a non-blocking key-value interface on the data: put, get, delete, and scan operations.

HBase closely matches the scale-out properties assumed for NoSQL databases. HTables are horizontally partitioned in regions that are assigned to Region-Servers nodes. In turn, each region is stored as an appendable file in the distributed file system, Hadoop File System (HDFS) [START_REF] Shvachko | The Hadoop Distributed File System[END_REF] based on GFS [START_REF] Ghemawat | The Google file system[END_REF]. By default, HBase uniformly distributes data among all available nodes in the systems.

Apache Derby is an open source relational database implemented entirely in Java. Derby has a small footprint, about 2.6 megabytes for the base engine and an embedded JDBC driver. In addition, it is easy to install, deploy and use.

Besides providing a complete implementation of SQL and JDBC, Derby has the advantage of already providing an embedded mode, which eases its use as a middleware layer.

The store layer of Derby is split into two main areas, access and raw. The access layer presents a conglomerate (table or index)/row based interface to the SQL layer. It handles table scans, index scans, index lookups, indexing, sorting, locking policies, transactions, isolation levels. The access layer sits on top of the raw store, which provides the raw storage of rows in pages in files, transaction logging, transaction management. Following the architecture proposed in the previous chapter, the raw store layer was removed in our prototype and some components of the access layer were replaced.

Prototype components

The system is composed of the following layers: (i) query engine, (ii) storage and (iii) file system. Applications issue SQL requests to any query engine node. The query engine node communicates with storage nodes, executes queries and returns the results to applications. We use Derby components to implement the query engine. We reuse its query processing sub-system, both the compiler and the optimizer components. Two new operators for index and sequential data scans have been added to the set of Derby's generic relational operators. These operators leverage HBase's indexing and filtering capabilities to minimize the amount of data that needs to be fetched. The SQL advanced operators such as JOIN and aggregations are not supported by HBase and are implemented at the query engine. The query engine translates the user queries into some appropriate put, get, delete, and scan operations to be invoked on HBase. Each HBase region is stored as an appendable file in the distributed file system. The distributed file system could be implemented using any scalable file system that supports append-only files such as HDFS.

Relational-tuple store mapping

Relational tables and secondary indexes are mapped to the HBase's data model. We have adopted a simple mapping from a relational table to a HTable. There is a one-to-one mapping where the HBase row's key is the relational primary key (simple or compound) and all relational columns are mapped into a single column family. Since relational columns are not multi-valued, each relational column is mapped to a HTable column. The schema of relational tables is rigid, i.e., every row in the same table must have the same set of columns. However, the value for some relational columns can be NULL and thus an HTable column for a given row exists only if its original relational column for that row is not NULL.

A secondary index of the relational model is mapped into an additional HTable. The additional table is necessary so that data is ordered by the indexed attributes. For each indexed attribute a HTable row is added and its row's key is the indexed attribute. For unique indexes the row has a single column with its value being the key of the matching indexed row in the primary key table. For non-unique indexes there is one column per matching indexed row with the name of the column being the matching row's key. Figure 2(a) depicts a relational table example. The column Number is the primary key and the table has two additional indexes: one unique index on attribute Telephone and a non-unique index on column Address. Therefore, the mapping will have three HTables: base data -Figure 2(b), unique index on column Telephone -Figure 2(c), and non-unique index on column Address -Figure 2(d).

Due to the simple mapping from the relational data model to HBase, the user can take advantage of the flexible schema exposed by the underlying NoSQL database and directly use its data in HBase for simple queries or complex mapreduce jobs.

Reducing data transfer

In order to reduce network traffic between the query engine and HBase, the implementation of sequential and index scan operators takes advantage of the indexing and filtering capabilities of HBase.

For index scans data is maintained ordered by one or more columns. This allows to restrict the desired rows for a given scan by optionally specifying the start and the stop keys. In a relational table each column is typed (e.g., char, date, integer, decimal, varchar) and data is ordered according to the natural order of the indexed column data type. However, row keys in HBase are plain byte arrays and neither Derby or HBase byte encoding preserve the data type's natural order. In order to build and store indexes in HBase maintaining the data type's order we need to map row keys into plain bytes in such a way that when HBase compares them the order of the data type is preserved. This mapping has been implemented for integer, decimal, char, varchar and date types. As indexes may be composite, besides each specific data type encoding, we also needed to define a way to encode multiple indexed columns in the same byte array. We do so by simply concatenating them from left to right, according to the order they are defined in the index using a pre-defined separator.

In HBase the start and stop keys of a scan must always refer to all the columns defined in the index. However, when using compound indexes the DQE may generate scans using subsets of the index columns. Indeed, an index scan can use equality conditions on any prefix of the indexed columns (from left to right) and at most one range condition on the rightmost queried column. In order to map these partial scans, the default start and stop keys in HBase are not used but instead the scan expression is run through HBase's BinaryPrefixComparator filter.

The aforementioned mechanisms reduce the traffic between the query engine and HBase by only retrieving the rows that match the range of the index scan. However, a scan can also select non-indexed columns. A naive implementation of this selection would fetch all rows from the index scan and test the relevant columns row by row. In detail, doing so on top of HBase would require a full table scan, which means fetching all the table rows from the different regions and possible different RegionServers. The full table would therefore be brought to the query engine instance and only then discard those rows not selected by the filter. To mitigate this performance overhead, particularly for low selective queries, that this approach may incur, the whole selection is pushed down into HBase. This is done by using the SingleColumnValueFilter filter to test a single column and, to combine them respecting the conjunctive normal form, using the FilterList filter. The latter represents an ordered list of filters that are evaluated with a specified boolean operator FilterList.Operator.MUST PASS ALL (AND) or FilterList.Operator.MUST PASS ONE (OR).

Metadata

Derby must also store information about the in-memory representation of tables and index, conglomerates, in a persistent manner.

For this we use a special HBase table, ConglomerateInfo, with information for all available conglomerates. ConglomerateInfo has a single column family, MetaInfo, and data for each conglomerate is stored in a row whose key is the conglomerate's identifier. Three columns are used to save information: value, a byte array with the encoding of all information; name, the name of the index or table to which this conglomerate matches; and size, to store the estimate of the current size (number of rows) of the conglomerate.

The information stored in ConglomerateInfo is mainly modified by Data Definition Language (DDL) statements. However, the size attribute is modified by any update DML statement (INSERT, UPDATE or DELETE) and therefore this attribute may change frequently. Thus, we used a simple mechanism to update the value in a distributed and efficient manner. The size of the table/index is used only for the query optimizer to decide the best query plan and therefore don't need to be the most recent value. In most RDBMS, this value is maintained probabilistically manner and thus we update its value in each query engine instance asynchronously.

Each query engine instance maintains the last estimate it has for the global size (shared by all instances), updates it accordingly to the local changes (when some update DML statement occurs), and periodically updates the value stored in ConglomerateInfo (shared by all instances) and refresh its local estimate. For, this it maintains the delta of the size after the last update to Conglomer-ateInfo, and in the next update it increments or decrements the size stored in ConglomerateInfo with the delta value, using a special HBase method for this purpose (incrementColumnValue ). Then, it resets the delta value for the next time window.

Evaluation

The evaluation of our prototype addresses two performance aspects: the overhead imposed by the DQE in terms of added latency, and the system scale-out in terms of the achieved throughput by increasing the number of nodes.

Overhead

We measure the increased latency of the system resulting from using the DQE instead of a standard HBase client.

Test workload We used a workload typical of NoSQL databases, Yahoo! Cloud Serving Benchmark [START_REF] Cooper | Benchmarking cloud serving systems with YCSB[END_REF]. YCSB was designed to benchmark the performance of NoSQL databases under different workloads. It has client implementations for several NoSQL databases and a JDBC client for RDBMS. We have used the HBase and JDBC clients without any modifications.

Experimental setting

The machines used for the experiments had a 2.4 GHz Dual-Core AMD Opteron(tm) Processor, with 4GB of RAM and a local SATA hard-disk. The machines were interconnected by a switched Gigabit local area network.

For the experiments 2 machines were used: one to run the workload generator, using an embedded connection to the modified Derby; and another running HBase. HBase was run in standalone mode, meaning that the HBase master and HBase RegionServer were collocated in the same machine using the local filesystem.

The YCSB database was populated with 100,000 rows (185MB) and the workload consisted of 1,000,000 operations. The proportion for the different types of operations was 60% reads, 20% updates and 20% scans. The operations were distributed uniformly over the database rows. The size of the scan operator was also a uniform random number between 1 and 10. We measured two runs where each client had 1 and 50 threads, respectively. In both, the target throughput was 100 operations per second.

Results

The average latency (in milliseconds) for the YCSB workload is shown in Table 1, for the standard HBase client and DQE.

The results for the insert operations correspond to the loading of the database while other operations are due to the mix of operations generated by the workload itself.

The results show that for all types of operations the query engine can be embedded in the application with an overhead totally in line with the base figures. The additional latency is due to the SQL processing and required marshalling/unmarshaling. Moreover, the overhead of DQE decreases with the increasing number of concurrent clients (threads).

Scale-out

We measured the increased throughput of the system when varying the number RegionServer nodes from 1 to 30.

Test Workload For the evaluation of the scale-out we used the load of an industry standard on-line transaction processing SQL benchmark, TPC-C. 1 It mimics a wholesale supplier with a number of geographically distributed sales districts and associated warehouses. The warehouses are hotspots of the system and the benchmark defines 10 client per warehouse. TPC-C specifies five transactions: NewOrder with 44% of the occurrences; Payment with 44%; OrderStatus with 4%; Delivery with 4%; and StockLevel with 4%. The NewOrder, Payment and Delivery are update transactions while the others are read-only. The traffic is a mixture of 8% read-only and 92% update transactions and therefore is a write intensive workload.

We have used both an existing SQL implementation,2 without modifications, to drive the DQE, and an existing NoSQL implementation optimized for HBase. 3Briefly, in the NoSQL implementation, TPC-C columns are grouped into column families, named differently for optimization, and the data storage layout has been optimized.

Experimental Setting

We ran all the experiments on a cluster of 42 machines with 3.10GHz GHz Dual-Core i3-2100 CPU, with 4GB of RAM and a local SATA hard-disk. The machines were interconnected by a switched Gigabit local area network.

The TPC-C workload has run from a varying number of machines. For our proposal, we varied the number of client machines from 1 to 10, each running 150 client threads. Each client machine as also running a DQE instance as a middleware layer.

One machine was used to run the HDFS namenode, HBase Master and Zookeper [START_REF] Hunt | Zookeeper: wait-free coordination for internet-scale systems[END_REF]. The remaining machines were RegionServers, each configured with a heap of 3GB, and also running a HDFS DataNode instance.

The TPC-C database was populated according to the number of Region-Servers, ranging from 5 warehouses, for a single RegionServer, to 150 warehouses, for 30 RegionServers. All TPC-C tables, were partitioned and distributed so there were 5 warehouses per RegionServer each handling a total of 50 clients. With 150 warehouses, the size of the database was about 75GB.

Results

The system throughput for a varying configuration of 1, 6, 12, 18, 24, and 30 RegionServers is depicted in Figure 3(a). The results show that DQE presents linear scalability. This is mainly due to the scale independence of the query processing layer and the scalability of the NoSQL database layer.

Furthermore, while with the DQE the system has a slightly lower throughput up tp 6 RegionServers its scalability is better than that of HBase under the NoSQL TPC-C load generator. This can be mainly attributed to the optimizations achieved by Derby's query engine that take advantage of relational operators, filtering and secondary indexes, while manual optimizations and denormalization still incur on greater complexity resulting in a greater overhead and worse scalability.

In this specific case, the query engine can greatly restrict the amount of data retrieved from HBase by, as previously explained, also taking advantage of HBase filters to select non-indexed columns. As a matter of fact, the network traffic when using the DQE is much lower than using the NoSQL TPC-C implementation. In fact, this prevented us from getting results for PyTPCC, the implementation optimized for HBase, with more than 18 RegionServers due to the saturation of the network. The latency figures in Figure 3(b) reflect the problem very clearly.

Related Work

Existing NoSQL databases rely on simplified and heterogenous query interfaces, that constitute a barrier on their adoption. Projects like BigQuery [START_REF]Google[END_REF], Hive [2] or Tenzing [START_REF] Lin | Tenzing: A SQL implementation on the MapReduce framework[END_REF] try to mitigate this constraint by providing an interface based on SQL over a MapReduce framework and a key/value store, but are mainly intended for data warehousing and analytical purposes.

BigQuery is a Google web service, built on top of BigTable. As query language, it uses a SQL dialect, that is a variation of the standard. It only offers a subset of the standard operators like selection, projection, aggregation and ordering. Joins and other more complex operators are not supported. In addition, data is immutable once uploaded to BigTable. Hive is built on top of Hadoop, a project that encompasses HDFS and the MapReduce framework. Hive also defines a simple SQL-like query language to query data but it offers more complex operators such as equi-joins, which are converted into MapReduce jobs, and unions. Likewise, Tenzing relies on a MapReduce framework to provide a SQL query execution engine, offering a mostly complete SQL implementation. The Hadapt commercial system4 (previously HadoopDB [START_REF] Abouzeid | HadoopDB: an architectural hybrid of MapReduce and DBMS technologies for analytical workloads[END_REF]) is also an analytical driven database, but it takes a slightly different approach by providing a hybrid system. Like Hive, it uses an SQL interface over the MapReduce framework from Hadoop, but replaces the HDFS layer with a cluster of single-node relational databases.

Like Hadapt, the CloudDB [START_REF] Hacigümüs | CloudDB: One Size Fits All Revived[END_REF] project is a hybrid system. However, it supports both OLAP and OLTP by providing three types of data storage systems: a relational database, a NoSQL database and a database oriented for OLAP. Data is stored in the database, according to the guarantees of data consistency required by the user.

Similarly to our approach, PIQL [START_REF] Armbrust | PIQL: success-tolerant query processing in the cloud[END_REF] and Megastore [START_REF] Baker | Megastore: Providing Scalable, Highly Available Storage for Interactive Services[END_REF] propose an architecture with higher-level processing functionality via a database library.

In PIQL, the application issues queries in a new declarative language that is based on a subset of SQL, but extended with new statements and some new operators to always achieve a predictable performance independently of the database size (i.e. scale-independence). However, there are several restrictions on the supported operations. For instance, a table scan is not scale-independent and has to be appended with a limit statement to bound the results. While this is done to achieve predictable performance it is a major impediment to run legacy applications.

MegaStore is built on top of BigTable and implements some of the features of RDBMS, such as secondary indexing. Nonetheless, join operations must be implemented at the application side. Therefore, applications must be written specifically for MegaStore using its data model and queries are restricted to scans and lookups.

The use of a library-centric component to offer higher-level processing functionality in our prototype is similar to the architecture of PIQL, Megastore, as well as that by Brantner et al. [START_REF] Brantner | Building a database on S3[END_REF]. However, our architecture allows to combine the expressiveness and performance of RDBMS, taking advantage of its query optimizer and full SQL support. This allows our proposal to run existing SQL applications and tools while retaining the scale-out capabilities of the NoSQL database.

On a different perspective, other approaches make use of object mapping tools that allow to bypass the database lower level interfaces. By using [START_REF] Gomes | An object mapping for the Cassandra distributed database[END_REF], the user has at her disposal the generic object interfaces like JPA and JDO that allow her to use NoSQL databases in an almost transparent way, leveraging the knowledge already existent in the area. These solutions have also the advantage of aiding the migration of existent solutions based on object to relation mappers allowing the mix of different types of databases under the same code base.

Conclusions and Future Work

We presented a distributed query engine allowing to execute ANSI compliant SQL queries on top of a NoSQL database. The query engine allows to combine the expressiveness and performance of a Relational Database Management System with the scalability and schema flexibility of a NoSQL database.

The developed prototype offers a standard JDBC client interface and can be embedded in the client application as a middleware layer. It is stateless with respect to application data making its replication straightforward. For the execution of data management language statements it does not require any kind of coordination which prevents any undesirable impact on the scalability of the underlying NoSQL database. The feasibility of the approach is demonstrated by the performance results obtained with the YCSB and TPC-C benchmarks.

Moreover, the comparison with a NoSQL TPC-C implementation optimized for HBase shows that the presented prototype, through the use of the full-fledged query engine from Apache Derby achieves impressive performance.

Fig. 1 .

 1 Fig. 1. Data management architectures.

Fig. 2 .

 2 Fig. 2. Data Model Mapping

Fig. 3 .

 3 Fig. 3. TPC-C results

Table 1 .

 1 Overhead results (ms)

	Workload 1 thread, 100 tps 50 threads, 100 tps
	Operation HBase DQE HBase	DQE
	Insert	0.58	0.93	1.04	1.98
	Update 0.51	1.3	2.66	3.1
	Read	0.53	0.79	1.63	1.7
	Scan	1.43	2.9	4.64	6.1

Since the current system does not include a transaction manager, all transactional contexts of the benchmark are simply discarded.

BenchmarkSQLhttp://sourceforge.net/projects/benchmarksql/

https://github.com/apavlo/py-tpcc/wiki/HBase-Driver

http://www.hadapt.com/

Acknowledgment

This work is part-funded by: ERDF -European Regional Development Fund through the COMPETE Programme (operational programme for competitiveness) and by National Funds through the FCT -Fundação para a Ciência e a Tecnologia (Portuguese Foundation for Science and Technology) within project Stratus/FCOMP-01-0124-FEDER-015020; and European Union Seventh Framework Programme (FP7) under grant agreement n • 257993 (CumuloNimbo).