Sandro Franceschi

Irina Kurkova

Kilian Raschel

Analytic approach for reflected Brownian motion in the quadrant

Keywords: Reflected Brownian motion in the quadrant, (Random) Walks in the quarter plane, Stationary distribution, Laplace transform, Generating function, Boundary value problem, Asymptotic analysis

Random walks in the quarter plane are an important object both of combinatorics and probability theory.

Of particular interest for their study, there is an analytic approach initiated by Fayolle, Iasnogorodski and Malyšev, and further developed by the last two authors of this note. The outcomes of this method are explicit expressions for the generating functions of interest, asymptotic analysis of their coefficients, etc. Although there is an important literature on reflected Brownian motion in the quarter plane (the continuous counterpart of quadrant random walks), an analogue of the analytic approach has not been fully developed to that context. The aim of this note is twofold: it is first an extended abstract of two recent articles of the authors of this paper, which propose such an approach; we further compare various aspects of the discrete and continuous analytic approaches.

Introduction

Random walks in the quarter plane

Since the seventies and the pioneered papers [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF]; [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF], random walks in the quarter plane (cf. Figure 1) are extensively studied. They are indeed an important object of probability theory and have been studied for their recurrence/transience, for their links with queueing systems [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]), representation theory [START_REF] Biane | Minuscule weights and random walks on lattices[END_REF]), potential theory. Moreover, the state space N 2 = {0, 1, 2, . . .} 2 offers a natural framework for studying any two-dimensional population; accordingly, quadrant walks appear as models in biology and in finance [START_REF] Cont | Price dynamics in a Markovian limit order market[END_REF]). Another interest of random walks in the quarter plane is that in the large class of random processes in cones, they form a family for which remarkable exact formulas exist. Moreover, quadrant walks are popular in combinatorics, see Bousquet-Mélou and Mishna (2010); [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF]; [START_REF] Kurkova | On the functions counting walks with small steps in the quarter plane[END_REF]. Indeed, many models of walks are in bijection with other combinatorial objects: maps, permutations, trees, Young tableaux, etc. In combinatorics again, famous models have emerged from quadrant walks, as Kreweras' or Gessel's ones, see Bousquet-Mélou and Mishna (2010); [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF]. Finally, walks in the quarter plane are interesting for the numerous tools used for their analysis: combinatorial (Bousquet-Mélou and Mishna (2010)), from complex analysis [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF]; [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]; [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF]; Kurkova andRaschel (2011, 2012); [START_REF] Bernardi | Counting quadrant walks via Tutte's invariant method[END_REF]), computer algebra [START_REF] Bostan | The complete generating function for Gessel walks is algebraic[END_REF]), for instance.

Issues and technicalities of the analytic approach

In the literature (see, e.g., [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF]; [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]; [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF]; Kurkova andRaschel (2011, 2012)), the analytic approach relies on six key steps:

(i) Finding a functional equation between the generating functions of interest;

(ii) Reducing the functional equation to boundary value problems (BVP);

(iii) Solving the BVP;

(iv) Introducing the group of the walk;

(v) Defining the Riemann surface naturally associated with the model, continuing meromorphically the generating functions and finding the conformal gluing function;

(vi) Deriving the asymptotics of the (multivariate) coefficients.

Before commenting these different steps, let us note that altogether, they allow for studying the following three main problems:

(P1) Explicit expression for the generating functions of interest (needs (i), (ii), (iii) and (v));

(P2) Algebraic nature of these functions (needs (iv) and (v));

(P3) Asymptotics of their coefficients in various regimes (needs (iii), (v) and (vi)).

The point (i) reflects the inherent properties of the model and is easily obtained. Point (ii), first shown in [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF], is now standard (see [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF]) and follows from algebraic manipulations of the functional equation of (i). Item (iii) uses specific literature devoted to BVP (our main reference for BVP is the book of [START_REF] Litvinchuk | Solvability Theory of Boundary Value Problems and Singular Integral Equations with Shift[END_REF]). The idea of introducing the group of the model (iv) was proposed in [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF], and brought up to light in the combinatorial context in Bousquet-Mélou and Mishna (2010). Point (v) is the most technical a priori; it is however absolutely crucial, as it allows to access key quantities (as a certain conformal gluing function which appears in the exact formulation of (iii)). Finally, (vi) uses a double refinement of the classical saddle point method: the uniform steepest descent method.

Reflected Brownian motion in the quarter plane

There is a large literature on reflected Brownian motion in quadrants (and also in orthants, generalization to higher dimension of the quadrant), to be rigorously introduced in Section 3. First, it serves as an approximation of large queuing networks (see [START_REF] Foddy | Analysis of Brownian motion with drift, confined to a quadrant by oblique reflection (diffusions, Riemann-Hilbert problem)[END_REF]; [START_REF] Baccelli | Analysis of models reducible to a class of diffusion processes in the positive quarter plane[END_REF]); this was the initial motivation for its study. In the same vein, it is the continuous counterpart of (random) walks in the quarter plane. In other directions, it is studied for its Lyapunov functions in [START_REF] Dupuis | Lyapunov functions for semimartingale reflecting Brownian motions[END_REF], cone points of Brownian motion in Le [START_REF] Gall | Mouvement brownien, cônes et processus stables[END_REF], intertwining relations and crossing probabilities in [START_REF] Dubédat | Reflected planar Brownian motions, intertwining relations and crossing probabilities[END_REF], and of particular interest for us, for its recurrence/transience in [START_REF] Hobson | Recurrence and transience of reflecting Brownian motion in the quadrant[END_REF]. The asymptotic behavior of the stationary distribution (when it exists) is now well known, see [START_REF] Harrison | Reflected Brownian motion in the quadrant: tail behavior of the stationary distribution[END_REF]; [START_REF] Dai | Reflecting Brownian motion in two dimensions: Exact asymptotics for the stationary distribution[END_REF]; [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF]. There exist, however, very few results giving exact expressions for the stationary distribution. Let us mention [START_REF] Foddy | Analysis of Brownian motion with drift, confined to a quadrant by oblique reflection (diffusions, Riemann-Hilbert problem)[END_REF] (dealing with the particular case of a Brownian motion with identity covariance matrix), [START_REF] Baccelli | Analysis of models reducible to a class of diffusion processes in the positive quarter plane[END_REF] (on a diffusion having a quite special behavior on the boundary), [START_REF] Harrison | Multidimensional reflected Brownian motions having exponential stationary distributions[END_REF]; [START_REF] Dieker | Reflected Brownian motion in a wedge: sum-of-exponential stationary densities[END_REF] (on the special case when stationary densities are exponential) and [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF] (on the particular case of an orthogonal reflection). We also refer to [START_REF] Burdzy | Obliquely reflected Brownian motion in non-smooth planar domains[END_REF] for the analysis of reflected Brownian motion in bounded planar domains by complex analysis techniques.

Main results and plan

This note is an extended abstract of the papers [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF]; [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF], whose main contributions are precisely to export the analytic method for reflected Brownian motion in the quarter plane. Our study constitutes one of the first attempts to apply these techniques to the continuous setting, after [START_REF] Foddy | Analysis of Brownian motion with drift, confined to a quadrant by oblique reflection (diffusions, Riemann-Hilbert problem)[END_REF]; [START_REF] Baccelli | Analysis of models reducible to a class of diffusion processes in the positive quarter plane[END_REF]. In addition of reporting about the works [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF]; [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF], we also propose a comparative study of the discrete/continuous cases.

Our paper is organized as follows: Section 2 concerns random walks and Section 3 Brownian motion. For clarity of exposition we have given the same structure to Sections 2 and 3: in Section 2.1/3.1 we first state the key functional equation (a kernel equation), which is the starting point of our entire analysis. We study the kernel (a second degree polynomial in two variables). In Section 2.2/3.2 we state and solve the BVP satisfied by the generating functions. We then move to asymptotic results (Section 2.3/3.3). In Section 2.4/3.4 we introduce the Riemann surface of the model and some important related facts.

Random walks in the quarter plane

This section is devoted to the discrete case and is based mainly on [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF].

Functional equation

One considers a piecewise homogeneous random walk with sample paths in N 2 . There are four domains of spatial homogeneity (the interior of N 2 , the horizontal and vertical axes, the origin), inside of which the transition probabilities (of unit size) are denoted by p i,j , p i,j , p i,j and p 0 i,j , respectively. See Figure 1. The inventory polynomial of the inner domain is called the kernel and equals K(x, y) = xy{ -1 i,j 1 p i,j x i y j -1}.

(1)

The inventory polynomials associated to the other homogeneity domains are

k(x, y) = x{ p i,j x i y j -1}, k(x, y) = y{ p i;j x i y j -1}, k 0 (x, y) = { p 0 i,j x i y j -1}.
Assuming the random walk ergodic (we refer to [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF], Theorem 1.2.1) for necessary and sufficient conditions), we denote the invariant measure by {π i,j } i,j 0 and introduce the generating functions π(x, y) = i,j 1 π i,j x i-1 y j-1 , π(x) = i 1 π i,0 x i-1 , π(y) = j 1 π 0,j y j-1 .

S. Franceschi, I. Kurkova and K. Raschel Writing the balance equations at the generating function level, we have (see (Fayolle et al., 1999, Equation (1.3.6)) for the original statement):

Lemma 1 The fundamental functional equation holds

-K(x, y)π(x, y) = k(x, y)π(x) + k(x, y) π(y) + k 0 (x, y)π 0,0 . (2)
Equation ( 2) holds a priori in the region {(x, y) ∈ C 2 : |x| 1, |y| 1}. Indeed, the π i,j sum up to 1, so that the generating functions π(x, y), π(x) and π(y) are well defined on the (bi)disc. The identity (2) is a kernel equation, and a crucial role will be played by the kernel (1). This polynomial K is of second order in both x and y; its roots X(y) and Y (x) defined by

K(X(y), y) = K(x, Y (x)) = 0 (3)
are thus algebraic of degree 2. Writing the kernel as K(x, y) = a(y)x 2 + b(y)x + c(y) and defining its discriminant d(y) = b(y) 2 -4a(y)c(y), one has obviously

X(y) = -b(y) ± d(y) 2a(y) .
The polynomial d has three or four roots, and exactly two of them are located in the unit disc, see (Fayolle et al., 1999, Lemma 2.3.8). They are named y 1 , y 2 , cf. Figure 2. On (y 1 , y 2 ) one has d(y) < 0, so that the two values (or branches) of X(y) (that we shall call X 0 (y) and X 1 (y)) are complex conjugate of one another. In particular, the set

M = X([y 1 , y 2 ]) = {x ∈ C : K(x, y) = 0 and y ∈ [y 1 , y 2 ]}
is symmetrical w.r.t. the real axis (Figure 2). This curve will be used to set a boundary condition for the unknown function π (Lemma 2).

Statement and resolution of the BVP

The analytic approach of [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF]; [START_REF] Fayolle | Two coupled processors: the reduction to a Riemann-Hilbert problem[END_REF]; [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF] proposes a way for solving the functional equation (2), by reduction to a BVP. Generally speaking, a BVP consists of a regularity condition and a boundary condition.

Lemma 2 The function π satisfies the following BVP:

• π is meromorphic in the bounded domain delimitated by M and has there identified poles;

• for any x ∈ M, k(x, Y 0 (x)) k(x, Y 0 (x)) π(x) - k(x, Y 0 (x)) k(x, Y 0 (x)) π(x) = k 0 (x, Y 0 (x)) k(x, Y 0 (x)) π 0,0 - k 0 (x, Y 0 (x)) k(x, Y 0 (x)) π 0,0 .
Proof: The regularity condition follows from Theorem 5, which provides a (maximal) meromorphic continuation of the function π. We now turn to the boundary condition. For i ∈ {0, 1}, we evaluate the functional equation ( 2) at (X i (y), y) and divide by k(X i (y), y). We then make the difference of the identities corresponding to i = 0 and i = 1. Finally, we substitute X 0 (y) = x and X 1 (y) = x, noting that when y ∈ [y 1 , y 2 ], x ∈ M by construction. Notice that we have chosen the segment [y 1 , y 2 ] connecting the points inside of the unit disc (Figure 2), in which we know that the generating function π is well defined. Although Lemma 2 could be written more precisely (by giving the number and the location of the poles of π), we shall prefer the above version, since we focus in this note on the main ideas of the analytic approach.

2 E T 1 y 1 y 2 y 3 y 4 t t t t t E T M t t t t
Lemma 2 happens to characterize the generating functions, as it eventually leads to an explicit expression for π, see Theorem 3. Before stating this central result (borrowed from (Fayolle et al., 1999, Theorem 5.4.3)), we need to introduce a function w called a conformal gluing function. By definition it satisfies w(x) = w(x) for x ∈ M and is one-to-one inside of M. This function will be constructed in Theorem 6 of Section 2.4.

Theorem 3 There exist two functions f and g, constructed from k, k, k 0 and w, such that the following integral formulation for π holds:

π(x) = f (x) M g(u) w (u) w(u) -w(x) du.
A similar contour integral representation exists for π, and eventually the functional equation ( 2) provides us with an explicit expression for the bivariate function π(x, y).

Asymptotics of the stationary probabilities

The asymptotics of coefficients π i,j of unknown generating functions satisfying the functional equation ( 2) has been obtained by [START_REF] Malyšev | Asymptotic behavior of the stationary probabilities for two-dimensional positive random walks[END_REF] via analytic arguments. He computed the asymptotics of the stationary probabilities π i,j as i, j → ∞ and j/i = tan α, for any given α ∈ (0, π/2). Let us briefly present these results. It is assumed in [START_REF] Malyšev | Asymptotic behavior of the stationary probabilities for two-dimensional positive random walks[END_REF] that the random walk is simple, meaning that

p -1,1 = p 1,1 = p -1,-1 = p 1,-1 = 0. (4)
It is also assumed that both coordinates of the interior drift vector are negative (as in Figure 3). For α ∈ (0, π/2), we define the point (x(α), y(α)) as follows. Introducing as in [START_REF] Kurkova | Random walks in (Z + ) 2 with non-zero drift absorbed at the axes[END_REF] the function P (u, v) = i,j p i,j e iu e jv on R 2 , the mapping

(u, v) → ∇P (u, v) |P (u, v)| is a homeomorphism between {(u, v) ∈ R 2 : P (u, v) = 1}
and the unit circle. The point (u(α), v(α)) is the unique solution to ∇P (u,v) |P (u,v)| = (cos α, sin α). Finally, (x(α), y(α)) = (e u(α) , e v(α) ). Following [START_REF] Malyšev | Asymptotic behavior of the stationary probabilities for two-dimensional positive random walks[END_REF], we introduce the sets of parameters P --= ({p i,j }, α) : k(ψ(x(α), y(α))) 0 and k(φ(x(α), y(α))) 0 ,

P +-= ({p i,j }, α) : k(ψ(x(α), y(α)
)) > 0 and k(φ(x(α), y(α))) 0 , and P -+ and P ++ accordingly. The automorphisms ψ and φ are defined in Section 2.4 by ( 8). The following theorem is proven in [START_REF] Malyšev | Asymptotic behavior of the stationary probabilities for two-dimensional positive random walks[END_REF].

Theorem 4 Let (i, j) = (r cos α, r sin α) with α ∈ (0, π/2). Then as r → ∞ we have

π i,j = (1 + o(1)) •          C0(α) √ r x -i (α)y -j (α) in P --, C 1 p -i 1 q -j 1 in P -+ , C 2 p -i 2 q -j 2 in P +-, C 1 p -i 1 q -j 1 + C 2 p -i 2 q -j 2 in P ++ , (5)
where C 0 , C 1 and C 2 are constants that can be expressed in terms of the functions π and π. The point (p 1 , q 1 ) is a solution of the system {K(x, y) = 0, k(ψ(x, y)) = 0} and similarly, (p 2 , q 2 ) is a solution of {K(x, y) = 0, k(φ(x, y)) = 0}.

Proof: The stationary probabilities π i,j are first written as two-dimensional Cauchy integrals, then reduced via the residue theorem to one-dimensional integrals along some contours. The asymptotics of these integrals is characterized either by the saddle point (x(α), y(α)) in the case of the set of parameters P --or by a pole (p 1 , q 1 ) or (p 2 , q 2 ) that is encountered when moving the integration contour to the saddle point; this happens for the sets of parameters P -+ , P +-and P ++ . 2

This approach has been applied for the analysis of the join the shortest queue problem in [START_REF] Kurkova | Malyshev's theory and JS-queues. Asymptotics of stationary probabilities[END_REF], and for computing the asymptotics of Green functions for transient random walks in the quarter plane reflected at the axes (see [START_REF] Kurkova | Martin boundary and elliptic curves. Markov Process[END_REF]) or killed at the axes (cf. [START_REF] Kurkova | Random walks in (Z + ) 2 with non-zero drift absorbed at the axes[END_REF]). Moreover, as illustrated in [START_REF] Kurkova | Random walks in (Z + ) 2 with non-zero drift absorbed at the axes[END_REF]; [START_REF] Kurkova | Malyshev's theory and JS-queues. Asymptotics of stationary probabilities[END_REF], the assumption (4) is not crucial for the applicability of the method. The limiting cases α = 0 and α = π/2 can also be treated via this approach, with some additional technical details (the saddle point then coincides with a branch point of the integrand), it is done in [START_REF] Kurkova | Random walks in (Z + ) 2 with non-zero drift absorbed at the axes[END_REF].

Riemann surface and related facts

In Section 2.1 the set

K = {(x, y) ∈ C 2 : K(x, y) = 0} = {(x, y) ∈ C 2 : p i,j x i y j = 1}
has appeared very naturally, since in order to state the BVP (our Lemma 2), we introduced the functions X(y) and Y (x), which by construction cancel the kernel, see (3).

In this section the central idea is to consider the (global) complex structure of K. The set K turns out to be a Riemann surface of genus 1, i.e., a torus. This simply comes from the reformulation of the identity

K(x, y) = 0 as {2a(y)x + b(y)} 2 = d(y).
Moreover, the Riemann surface of the square root of a polynomial of degree 3 or 4 is classically a torus (with this terminology, the roots of the discriminant are branch points). This new point of view on K brings powerful tools. Of particular interest is a parametrization of K in terms of Weierstrass elliptic functions:

K = {(x(ω), y(ω)) : ω ∈ C/(ω 1 Z + ω 2 Z)}. ( 6
)
This parametrization is totally explicit: x(ω) and y(ω) are rational functions in the ℘-Weierstrass function and its derivative ℘ (see (Fayolle et al., 1999, Lemma 3.3.1)); the periods ω 1 and ω 2 admit expressions as elliptic integrals in terms of {p i,j } (cf. (Fayolle et al., 1999, Lemma 3.3.2)), etc. Moreover, as any functions of x and/or y, the functions π(x) and π(y) can be lifted on K by setting

Π(ω) = π(x(ω)), Π(ω) = π(y(ω)). (7)

Group of the walk

Introduced in [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF] in a probabilistic context and further used in [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF]; Bousquet-Mélou and Mishna (2010), the group of the walk is a dihedral group generated by

ζ(x, y) = x, i p i,-1 x i i p i,+1 x i 1 y , η(x, y) = j p -1,j y j j p +1,j y j 1 x , y . (8)
(One easily verifies that these generators are idempotent:

ζ 2 = η 2 = 1.)
The group ζ, η can be finite or infinite, according to the order of the element ζ • η. The generator ζ (resp. η) exchanges the roots in y (resp. in x) of K(x, y) = 0. Viewed as a group of birational transformations in Bousquet-Mélou and Mishna (2010), we shall rather see it as a group of automorphisms of the Riemann surface K. This group has many applications. First, it allows for a continuation of the functions π(x) and π(y) (Theorem 5 below). It further connects the algebraic nature of the generating functions to the (in)finiteness of the group (Theorem 7). Finally, in the finite group case, elementary algebraic manipulations of the functional equations can be performed (typically, via the computation of certain orbit-sums) so as to eventually obtain D-finite expressions for the unknowns, see (Fayolle et al., 1999, Chapter 4) and Bousquet-Mélou and Mishna (2010).

Using the structure of the automorphisms of a torus, the lifted versions of ζ and η admit simple expressions (Fayolle et al., 1999, Section 3.1.2):

ζ(ω) = -ω + ω 1 + ω 2 , η(ω) = -ω + ω 1 + ω 2 + ω 3 , (9)
where, as the periods ω 1 and ω 2 , ω 3 ∈ (0, ω 2 ) is an elliptic integral (Fayolle et al., 1999, Lemma 3.3.3). Accordingly, the group is finite if and only if ω 2 /ω 3 ∈ Q, which provides a nice criterion in terms of elliptic integrals.

Continuation

While the generating function π(x) is defined through its power series in the unit disc, it is a priori unclear how to continue it to a larger domain. This is however crucial, since the curve M on which it satisfies a BVP (Lemma 2) is not included in the unit disc in general.

Theorem 5 The function π can be continued as a meromorphic function to

C \ [x 3 , x 4 ].
We notice that M does not intersect [x 3 , x 4 ] by (Fayolle et al., 1999, Theorem 5.3.3), so that Theorem 5 indeed provides a continuation of the generating function in the domain delimitated by M.

Proof: This result, stated as Theorem 3.2.3 in [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF], is a consequence of a continuation of the lifted generating functions (7) on the Riemann surface (or, better, on its universal covering -but we shall not go into these details here). The continuation on K uses the (lifted) functional equation ( 2) and the group of the walk ζ, η . 2

Conformal mapping

In the integral expression of Theorem 3, the conformal gluing function w is all-present, as it appears in the integrand and in f and g as well. The introduction of the Riemann surface K allows to derive an expression for this function (this is another major interest of introducing K). Let us recall that ω 1 and ω 2 are the periods of the elliptic functions of the parametrization (6), while ω 3 comes out in the lifted expression (9) of the automorphisms.

Theorem 6 The conformal gluing function w admits the expression:

w(x) = ℘(℘ -1 (x; ω 1 , ω 2 ); ω 1 , ω 3 ),
where for i ∈ {2, 3}, ℘( • ; ω 1 , ω i ) is the ℘-Weierstrass elliptic function with periods ω 1 and ω i .

Proof: While on the complex plane, M is a quartic curve, it becomes on K much simpler (typically, a segment). This remark (which again illustrates all the benefit of having introduced the Riemann surface) is used in (Fayolle et al., 1999, Section 5.5.2) so as to obtain the above expression for w. 2

Algebraic nature of the generating functions

Recall that a function of one variable is D-finite if it satisfies a linear differential equation with polynomial coefficients.

Theorem 7 If the group is finite, the generating functions π(x) and π(y) are D-finite.

Proof: This follows from manipulations on the Riemann surface, see (Fayolle et al., 1999, Chapter 4). The D-finiteness is proved on R; refined results (in the combinatorial context of the enumeration of paths) can be found in Bousquet-Mélou and Mishna ( 2010), where the D-finiteness is proved on Q. 2

The converse of Theorem 7 is not shown in full generality. It is true in combinatorics, see [START_REF] Kurkova | On the functions counting walks with small steps in the quarter plane[END_REF].

Reflected Brownian motion in the quadrant

Defining reflected Brownian motion in the quadrant

The object of study here is the reflected Brownian motion with drift in the quarter plane

R 2 + Z(t) = Z 0 + W (t) + µt + RL(t), ∀t 0, (10)
associated to the triplet (Σ, µ, R), composed of a non-singular covariance matrix, a drift and a reflection matrix, see Figure 3:

Σ = σ 11 σ 21 σ 12 σ 22 , µ = µ 1 µ 2 , R = (R 1 , R 2 ) = r 11 r 21 r 12 r 22 .
In Equation ( 10), Z 0 is any initial point in R 2 + , the process (W (t)) t 0 is an unconstrained planar Brownian motion starting from the origin, and for i ∈ {1, 2}, L i (t) is a continuous non-decreasing process, that increases only at time t such that Z i (t) = 0, namely t 0 1 {Z i (s) =0} dL i (s) = 0. The columns R 1 and R 2 represent the directions in which the Brownian motion is pushed when the axes are reached.

The reflected Brownian motion (Z(t)) t 0 associated with (Σ, µ, R) is well defined, see for instance [START_REF] Williams | Semimartingale reflecting Brownian motions in the orthant[END_REF]. Its stationary distribution exists and is unique if and only if the following (geometric flavored) conditions are satisfied (see, e.g., Harrison and Williams (1987a); [START_REF] Hobson | Recurrence and transience of reflecting Brownian motion in the quadrant[END_REF])

r 11 > 0, r 22 > 0, r 11 r 22 -r 12 r 21 > 0, r 22 µ 1 -r 12 µ 2 < 0, r 11 µ 2 -r 21 µ 1 < 0. ( 11
)
More that the Brownian motion in the quadrant, all results presented below concern the Brownian motion in two-dimensional cones (by a simple linear transformation of the cones). This is a major difference and interest of the continuous case, which also illustrates that the analytic approach is very well suited to that context.

Functional equation

Laplace transform of the stationary distribution

The continuous analogues of the generating functions are the Laplace transforms. As their discrete counterparts, they characterize the stationary distribution. Under assumption (11), that we shall do throughout the manuscript, the stationary distribution is absolutely continuous w.r.t. the Lebesgue measure, see Harrison and Williams (1987a); [START_REF] Dai | Steady-state analysis of reflected Brownian motions: Characterization, numerical methods and queueing applications[END_REF]. We denote its density by π(x) = π(x 1 , x 2 ). Let the Laplace transform of π be defined by

ϕ(θ) = E π [e θ|Z ] = R 2 + e θ|x π(x)dx.
We further define two finite boundary measures ν 1 and ν 2 with support on the axes, by mean of the formula

ν i (B) = E π 1 0 1 {Z(t)∈B} dL i (t) .
The measures ν i are continuous w.r.t. the Lebesgue measure by Harrison and Williams (1987a), and may be viewed as boundary invariant measures. We define their moment Laplace transform by

ϕ 2 (θ 1 ) = R+ e θ1x1 ν 2 (x 1 )dx 1 , ϕ 1 (θ 2 ) = R+ e θ2x2 ν 1 (x 2 )dx 2 .

Functional equation

There is a functional equation between the Laplace transforms ϕ, ϕ 1 and ϕ 2 , see ( 12), which is reminiscent of the discrete functional equation ( 2).

Lemma 8

The following key functional equation between the Laplace transforms holds

-γ(θ)ϕ(θ) = γ 1 (θ)ϕ 1 (θ 2 ) + γ 2 (θ)ϕ 2 (θ 1 ), ( 12
)
where

     γ(θ) = 1 2 θ|σθ + θ|µ = 1 2 (σ 11 θ 2 1 + σ 22 θ 2 2 + 2σ 12 θ 1 θ 2 ) + µ 1 θ 1 + µ 2 θ 2 , γ 1 (θ) = R 1 |θ = r 11 θ 1 + r 21 θ 2 , γ 2 (θ) = R 2 |θ = r 12 θ 1 + r 22 θ 2 .
By definition of the Laplace transforms, this equation holds at least for any θ = (θ 1 , θ 2 ) with θ 1 0 and θ 2 0. The polynomial γ in ( 12) is the kernel and is the continuous analogue of the kernel (1) in the discrete case. Polynomials γ 1 and γ 2 are the counterparts of k and k.

Proof: To show (12), the main idea is to use an identity called a basic adjoint relationship (first proved in Harrison and Williams (1987a) in some particular cases, then extended in [START_REF] Dai | Reflected Brownian motion in an orthant: numerical methods for steady-state analysis[END_REF]), which characterizes the stationary distribution. (It is the continuous analogue of the well-known equation πQ = 0, where π is the stationary distribution of a recurrent continuous-time Markov chain with infinitesimal generator Q.) This basic adjoint relationship connects the stationary distribution π and the corresponding boundary measures ν 1 and ν 2 . We refer to [START_REF] Foddy | Analysis of Brownian motion with drift, confined to a quadrant by oblique reflection (diffusions, Riemann-Hilbert problem)[END_REF]; [START_REF] Dai | Reflecting Brownian motion in two dimensions: Exact asymptotics for the stationary distribution[END_REF] for the details.

2

Elementary properties of the kernel

The kernel γ in ( 12) can be alternatively written as

γ(θ 1 , θ 2 ) = a(θ 2 )θ 2 1 + b(θ 2 )θ 1 + c(θ 2 ) = a(θ 1 )θ 2 2 + b(θ 1 )θ 2 + c(θ 1 ). ( 13
)
The equation γ(θ 1 , θ 2 ) = 0 defines a two-valued algebraic function Θ 1 (θ 2 ) by γ(Θ 1 (θ 2 ), θ 2 ) = 0, and similarly Θ 2 (θ 1 ) such that γ(θ 1 , Θ 2 (θ 1 )) = 0. Expressions of their branches are given by

Θ ± 2 (θ 1 ) = -b(θ 1 ) ± d(θ 1 ) 2a(θ 1 ) , where d(θ 1 ) = b 2 (θ 1 ) -4a(θ 1 )c(θ 1
) is the discriminant. The polynomial d has two zeros, real and of opposite signs; they are denoted by θ ± 1 and are branch points of the algebraic function Θ 2 . In the same way we define Θ ± 1 and its branch points θ ± 2 . Finally, notice that d is negative on R \ [θ - 1 , θ + 1 ]. Accordingly, the branches Θ ± 2 take complex conjugate values on this set.

Statement and resolution of the BVP

An important hyperbola

For further use, we need to introduce the curve

R = {θ 2 ∈ C : γ(θ 1 , θ 2 ) = 0 and θ 1 ∈ (-∞, θ - 1 )} = Θ ± 2 ((-∞, θ - 1 )). ( 14
)
It is the analogue of the curve M in Section 2.1. The curve R is symmetrical w.r.t. the real axis, see Figure 4 (this is a consequence of d being negative on (-∞, θ - 1 ), see above). Furthermore, it is a (branch of a) hyperbola by [START_REF] Baccelli | Analysis of models reducible to a class of diffusion processes in the positive quarter plane[END_REF]. We shall denote by G R the open domain of C bounded by R and containing 0, see Figure 4

. Obviously G R , the closure of G R , is equal to G R ∪ R.

BVP for orthogonal reflections

In the case of orthogonal reflections (see Figure 3), R is the identity matrix in (10), and we have γ 1 (θ 1 , θ 2 ) = θ 1 and γ 2 (θ 1 , θ 2 ) = θ 2 . We set

ψ 1 (θ 2 ) = 1 θ 2 ϕ 1 (θ 2 ), ψ 2 (θ 1 ) = 1 θ 1 ϕ 2 (θ 1 ). (15)
Fig. 4: Left: the discriminant d has two roots θ - 1 and θ + 1 of opposite signs. Right: the curve R in ( 14) is symmetric w.r.t. the horizontal axis and G R is the domain in blue Lemma 9 The function ψ 1 in (15) satisfies the following BVP:

(i) ψ 1 is meromorphic on G R with a single pole at 0, of order 1 and residue ϕ 1 (0), and vanishes at infinity;

(ii) ψ 1 is continuous on G R \ {0} and ψ 1 (θ 2 ) = ψ 1 (θ 2 ), ∀θ 2 ∈ R. (16)
Proof: The regularity condition of point (i) follows from Theorem 15, which provides a (maximal) meromorphic continuation of the function. Let us now consider (ii). Evaluating the (continued) functional equation ( 12) at (θ 1 , Θ ± 2 (θ 1 )), we obtain ψ 1 (Θ ± 2 (θ 1 )) + ψ 2 (θ 1 ) = 0, which immediately implies that

ψ 1 (Θ + 2 (θ 1 )) = ψ 1 (Θ - 2 (θ 1 )). (17)
Choosing θ 1 ∈ (-∞, θ - 1 ), the two quantities Θ + 2 (θ 1 ) and Θ - 2 (θ 1 ) are complex conjugate the one of the other, see Section 3.1. Equation ( 17) can then be reformulated as ( 16), using the definition (14) of the curve R.

2

The BVP stated in Lemma 9 is called a homogeneous BVP with shift (the shift stands here for the complex conjugation, but the theory applies to more general shifts, see [START_REF] Litvinchuk | Solvability Theory of Boundary Value Problems and Singular Integral Equations with Shift[END_REF]). It has a simpler form than the BVP in Lemma 2 for the discrete case, because there is no inhomogeneous term (as π 0,0 ) and also because in the coefficients in front of the unknowns there is no algebraic function (as Y 0 ) involved. Due to its particularly simple form, we can solve it in an explicit way, using the two following steps:

• Using a certain conformal mapping w (to be introduced below), we can construct a particular solution to the BVP of Lemma 9.

• The solution to the BVP of Lemma 9 is unique (see the invariant Lemma 2 in (Litvinchuk, 2000, Section 10.2)). In other words, two different solutions must coincide, and the explicit solution constructed above must be the function ψ 1 .

In [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF] it is explained that the above method may be viewed as a variation of Tutte's invariant approach, first introduced by Tutte for solving a functional equation arising in the enumeration of properly colored triangulations, see [START_REF] Tutte | Chromatic sums revisited[END_REF].

The function w glues together the upper and lower parts of the hyperbola R. There are at least two ways to find such a w. First, it turns out that in the literature there exist expressions for conformal gluing functions for relatively simple curves as hyperbolas, see (Baccelli and Fayolle, 1987, Equation (4.6)).

Here (based on Franceschi and Raschel ( 2016)), we use instead the Riemann sphere S, as we will see in Section 3.4. Indeed, many technical aspects (and in particular finding the conformal mapping) happen to be quite simpler on that surface.

We will deduce from Section 3.4 that function w can be expressed in terms of the generalized Chebyshev polynomial

T a (x) = cos(a arccos(x)) = 1 2 x + x 2 -1 a + x -x 2 -1 a as follows: w(θ 2 ) = T π β - 2θ 2 -(θ + 2 + θ - 2 ) θ + 2 -θ - 2 , (18)
where we have noted

β = arccos - σ 12 √ σ 11 σ 22 . (19)
In the case of orthogonal reflection, this methods leads to the main result of [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF], which is: Theorem 10 Let R be the identity matrix in (10). The Laplace transform ϕ 1 is equal to

ϕ 1 (θ 2 ) = -µ 1 w (0) w(θ 2 ) -w(0) θ 2 .

Statement of the BVP in the general case

We would like to close Section 3.2 by stating the BVP in the case of arbitrary reflections (non-necessarily orthogonal). Let us define for

θ 2 ∈ R G(θ 2 ) = γ 1 γ 2 (Θ - 1 (θ 2 ), θ 2 ) γ 2 γ 1 (Θ - 1 (θ 2 ), θ 2 ).
Similarly to Lemma 9, there is the following result: Lemma 11 The function ϕ 1 satisfies the following BVP:

(i) ϕ 1 is meromorphic on G R with at most one pole p of order 1 and is bounded at infinity;

(ii) ϕ 1 is continuous on G R \ {p} and ϕ 1 (θ 2 ) = G(θ 2 )ϕ 1 (θ 2 ), ∀θ 2 ∈ R. ( 20
)
Due to the presence of the function G = 1 in (20), this BVP (still homogeneous with shift) is more complicated than the one encountered in Lemma 9 and cannot be solved thanks to an invariant lemma. Instead, the resolution is less combinatorial and far more technical, and the solution should be expressed in terms of both Cauchy integrals and the conformal mapping w of Theorem 10. This will be achieved in a future work.

Asymptotics of the stationary probabilities

Overview

Let Π be a random vector that has the stationary distribution of the reflected Brownian motion. [START_REF] Dai | Reflecting Brownian motion in two dimensions: Exact asymptotics for the stationary distribution[END_REF] obtain the following asymptotic result: for a given directional vector c ∈ R 2 + they find (up to a multiplicative constant) a function f c (x) such that

lim x→∞ P[ c|Π x] f c (x) = 1.
In [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF] we solve a harder problem arisen in (Dai and Miyazawa, 2011, §8), namely computing the asymptotics of P[Π ∈ xc + B] as x → ∞, where c ∈ R 2 + is any directional vector and B ⊂ R 2 + any compact subset. Furthermore, we are able to find the full asymptotic expansion of the density π(x 1 , x 2 ) of Π as x 1 , x 2 → ∞ and x 2 /x 1 → tan(α), for any given angle α ∈ (0, π/2).

Main results

First we need to introduce some notations. The equation γ(θ) = 0 determines an ellipse E on R 2 passing through the origin, see Figure 5. Here we restrict ourselves to the case µ 1 < 0 and µ 2 < 0, although our methods can be applied without additional difficulty to other cases. For a given angle α ∈ [0, π/2], let us

The coordinates of θ(α) can be given explicitly. One can also construct θ(α) geometrically as on Figure 5. Secondly, consider the straight lines {γ 1 (θ) = 0} and {γ 2 (θ) = 0}, depending on the reflection matrix R only. They cross the ellipse E at the origin. The line {γ 1 (θ) = 0} (resp. {γ 2 (θ) = 0}) intersects the ellipse at a second point called θ * (resp. θ * * ). To present our results, we need to define the images on E of these points via the so-called Galois automorphisms ζ and η, to be introduced in Section 3.4. Namely, for the point θ * = (θ * 1 , θ * 2 ) ∈ E there exists a unique point ηθ * = (ηθ * 1 , θ * 2 ) ∈ E with the same second coordinate. Likewise, there exists a unique point ζθ * * = (θ * * 1 , ζθ * * 2 ) ∈ E with the same first coordinate as θ * * = (θ * * 1 , θ * * 2 ) ∈ E. Points θ * , θ * * , ηθ * and ζθ * * are pictured on Figure 5. Their coordinates can be made explicit.

Similarly to the discrete case, we introduce the set of parameters

Q --= ((Σ, µ, R), α) : γ 1 (ηθ(α)) < 0 and γ 2 (ζθ(α)) < 0 and Q +-, Q -+ and Q ++ accordingly.
The following theorem provides the main term in the asymptotic expansion of π(r cos α, r sin α).

Theorem 12 Let (x, y) = (r cos α, r sin α) with α ∈ (0, π/2). We assume that θ(α) ∈ R 2 + . Then as r → ∞ we have

π(re α ) = (1 + o(1)) •        C0 √ r e -r eα|θ(α) in Q --, C 1 e -r eα|ηθ * in Q +-, C 2 e -r eα|ζθ * * in Q -+ , C 1 e -r eα|ηθ * + C 2 e -r eα|ζθ * * in Q ++ , (22)
where C 0 , C 1 and C 2 are constants that can be expressed in terms of functions ϕ 1 and ϕ 2 and the parameters.

In [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF] the constants mentioned in Theorem 12 are specified in terms of functions ϕ 1 and ϕ 2 . But these functions are for now unknown. As we explained in Section 3.2, in a next work we are going to obtain ϕ 1 and ϕ 2 as solutions of BVP, thereby determining the constants in Theorem 12.

Proof of the key step of Theorem 12: Theorem 12 is proven in [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF]. The first step consists in continuing meromorphically the functions ϕ 1 and ϕ 2 on C \ [θ + 2 , ∞) or on the Riemann surface S, see Section 3.4. Then by the functional equation ( 12) and the inversion formula of Laplace transform (we refer to [START_REF] Doetsch | Introduction to the Theory and Application of the Laplace Transformation[END_REF] and [START_REF] Brychkov | Multidimensional Integral Transformations[END_REF]), the density π(x 1 , x 2 ) can be represented as a double integral. Using standard computations from complex analysis, we are able to reduce it to a sum of single integrals. We obtain the following (with the notation (13)):

π(x 1 , x 2 ) = -1 (2πi) 2 i∞ -i∞ i∞ -i∞ e -x1θ1-x2θ2 γ 1 (θ)ϕ 1 (θ 2 ) + γ 2 (θ)ϕ 2 (θ 1 ) γ(θ) dθ 1 dθ 2 = 1 2πi i∞ -i∞ ϕ 2 (θ 1 )γ 2 (θ 1 , Θ + 2 (θ 1 ))e -x1θ1-x2Θ + 2 (θ1) dθ 1 d(θ 1 ) + 1 2πi i∞ -i∞ ϕ 1 (θ 2 )γ 1 (Θ + 1 (θ 2 ), θ 2 )e -x1Θ + 1 (θ2)-x2θ2 dθ 2 d(θ 2 )
.

These integrals are typical to apply the saddle point method, see [START_REF] Fedoryuk | Asymptotic methods in analysis[END_REF]. The coordinates of the saddle point are the critical points of the functions cos(α)θ 1 + sin(α)Θ + 2 (θ 1 ) and cos(α)Θ + 1 (θ 2 ) + sin(α)θ 2 .

It is the point θ(α). Then we have to shift the integration contour up to new contours which coincide with the steepest-descent contour near the saddle point. When we shift the contours we have to take into account the poles of the integrands and their residues. The asymptotics will be determined by the pole if we cross a pole when we shift the contour and by the saddle point otherwise. 2

3.4 Riemann surface and related facts

Riemann surface

The Riemann surface

S = {(θ 1 , θ 2 ) ∈ C 2 : γ(θ 1 , θ 2 ) = 0}
may be viewed as the set of zeros of the kernel (equivalently, it is the Riemann surface of the algebraic functions Θ 2 and Θ 1 ). Due to the degree of γ, the surface S has genus 0 and is a Riemann sphere, i.e., homeomorphic to C ∪ {∞}, see [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF]. It admits a very useful rational parametrization, given by

θ 1 (s) = θ + 1 + θ - 1 2 + θ + 1 -θ - 1 4 s + 1 s , θ 2 (s) = θ + 2 + θ - 2 2 + θ + 2 -θ - 2 4 s e iβ + e iβ s , (23)
with β as in ( 19). The equation γ(θ 1 (s), θ 2 (s)) = 0 holds and S = {(θ 1 (s), θ 2 (s)) : s ∈ C ∪ {∞}}.

Group of the process

We finally introduce the notion of group of the model, similar to the notion of group of the walk in the discrete setting (see [START_REF] Malyšev | An analytic method in the theory of two-dimensional positive random walks[END_REF]; [START_REF] Fayolle | Random Walks in the Quarter-Plane[END_REF]; Bousquet-Mélou and Mishna ( 2010)). This group ζ, η is generated by ζ and η, given by (with the notation (13))

ζ(θ 1 , θ 2 ) = θ 1 , c(θ 1 ) a(θ 1 )

1 θ 2 , η(θ 1 , θ 2 ) = c(θ 2 ) a(θ 2 ) 1 θ 1 , θ 2 .
By construction, the generators satisfy γ(ζ(θ 1 , θ 2 )) = γ(η(θ 1 , θ 2 )) = 0 as soon as γ(θ 1 , θ 2 ) = 0. In other words, there are (covering) automorphisms of the surface S. Since ζ 2 = η 2 = 1, the group ζ, η is a dihedral group, which is finite if and only if the element ζ • η (or equivalently η • ζ) has finite order.

Algebraic nature of the Laplace transforms

With the above definition, it is not clear how to see if the group is finite, nor to see it its finiteness would have any implication on the problem. In fact, we have, with β defined in (19):

Lemma 13 The group ζ, η is finite if and only if π/β ∈ Q.

The proof of Lemma 13 is simple, once the elements ζ and η have been lifted and reformulated on the sphere S:

ζ(s) = 1 s , η(s) = e 2iβ s .
These transformations leave invariant θ 1 (s) and θ 2 (s), respectively, see (23). In particular, we have the following result (see [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF]), which connects the nature of the solution of the BVP to the finiteness of the group. Such a result holds for discrete walks, see our Theorem 7 and Bousquet-Mélou and Mishna (2010); [START_REF] Bernardi | Counting quadrant walks via Tutte's invariant method[END_REF].

Theorem 14 The solution ϕ 1 given in Theorem 10 and the conformal gluing function w in (18) are algebraic if and only if the group ζ, η is finite.

Conformal mapping

The conformal gluing function w introduced in Section 3.2 may be lifted on S. In fact its expression is even simpler using the parametrization of S. We show in [START_REF] Franceschi | Tutte's invariant approach for Brownian motion reflected in the quadrant[END_REF] that

where we make use of the principal determination of the logarithm.

Continuation of the Laplace transforms

To establish the BVP, we have stated a boundary condition for the functions ϕ 1 and ϕ 2 , on curves which lie outside their natural domains of definition (the half-plane with negative real-part), see Figure 4. In the same way, in the asymptotic study we use the steepest descent method on some curves outside of the initial domain of definition. We therefore need to extend the domain of definition of the Laplace transforms.

Theorem 15 The function ϕ 1 can be continued meromorphically on the cut plane C \ [θ + 2 , ∞).

Proof: The first step is to continue meromorphically ϕ 1 (θ 2 ) to the open and simply connected set {θ 2 ∈ C : θ 2 0 or Θ - 1 (θ 2 ) < 0}, by setting

ϕ 1 (θ 2 ) = γ 2 γ 1 (Θ - 1 (θ 2 ), θ 2 )ϕ 2 (Θ - 1 (θ 2 )).
This is immediate (see [START_REF] Franceschi | Asymptotic expansion of the stationary distribution for reflected Brownian motion in the quarter plane via analytic approach[END_REF] for the details). It is then possible to pursue the extension to the whole S using the invariance properties by the automorphisms ζ and η satisfied by the lifted Laplace transforms on S. 2

Fig. 1 :

 1 Fig. 1: Transition probabilities of the reflected random walk in the quarter plane, with four domains of spatial homogeneity

Fig. 2 :

 2 Fig. 2: Left: the polynomial d has three or four roots, denoted by y1, y2, y3, y4; exactly two of them are inside the unit disc. Right: the curve M = X([y1, y2]) is symmetrical w.r.t. the horizontal axis

Fig. 3 :

 3 Fig. 3: Drift µ and reflection vectors R 1 and R 2

5:

  Left: representation of the ellipse E, straight lines {γ1(θ) = 0}, {γ2(θ) = 0}, and points θ * , θ * * , ηθ * and ζθ * * . Right: geometric interpretation of the point θ(α) in (21) on E define the point θ(α) on the ellipse E by θ(α) = argmax θ∈E θ|e α ,where e α = (cos α, sin α).

  -s) + e -π β log(-s) ,

Acknowledgements

We acknowledge support from the projet MADACA of the Région Centre and from Simon Fraser University (British Columbia, Canada). We thank our three AofA referees for useful comments.


Martin Alain

Christine Guillemot

Dominique Thoreau

Philippe Guillotel

LEARNING CLUSTERING-BASED LINEAR MAPPINGS FOR QUANTIZATION NOISE REMOVAL

Keywords: clustering, linear mapping, regression, de-noising, HEVC

This paper describes a novel scheme to reduce the quantization noise of compressed videos and improve the overall coding performances. The proposed scheme first consists in clustering noisy patches of the compressed sequence. Then, at the encoder side, linear mappings are learned for each cluster between the noisy patches and the corresponding source patches. The linear mappings are then transmitted to the decoder where they can be applied to perform denoising. The method has been tested with the HEVC standard, leading to a bitrate saving of up to 9.63%.

INTRODUCTION

The main principles of modern video codecs, such as H.264 [START_REF] Wiegand | Overview of the H.264/AVC video coding standard[END_REF], HEVC [START_REF] Sullivan | Overview of the high efficiency video coding (HEVC) standard[END_REF] or VP9 [START_REF] Mukherjee | The latest open-source video codec VP9 -An overview and preliminary results[END_REF], first reside in the reduction of spatial and temporal redundancies, through prediction tools. The prediction residue is then transformed to further reduce the inner correlations of the signal, and the transformed coefficients are quantized to remove nonperceptive information. The quantized transformed coefficients are finally entropy coded to remove the remaining statistical redundancies. The improvements achieved over the years in video compression have been made by carefully optimizing each of these steps. For example, Intra prediction methods based on sparse representations have been proposed in [START_REF] Türkan | Sparse approximation with adaptive dictionary for image prediction[END_REF], and later based on Locally Linear Embedding (LLE) of patches [START_REF] Türkan | Image prediction based on neighbor embedding methods[END_REF]. The prediction methods based on LLE have then been extended to spatio-temporal prediction in [START_REF] Alain | Interprediction methods based on linear embedding for video compression[END_REF]. Recent works also investigate new transforms or quantization designs [START_REF] Arrufat | Mode-dependent transform competition for HEVC[END_REF] [START_REF] Valin | Perceptual vector quantization for video coding[END_REF]. Despite the obvious improvement brought by these techniques in terms of rate-distortion (RD) performances, the irreversible quantization step still results in a degradation of the reconstructed frames.

In this paper, we adopt a different perspective by proposing a novel and generic out-of-the-loop de-noising scheme, complementary to the aforementioned works. Contrary to most state-of-the-art de-noising methods, the proposed scheme does not rely on an additive white Gaussian noise (AWGN) model. In fact, it is known that the quantization noise fall outside these assumptions [START_REF] Guleryuz | Linear, worst-case estimators for denoising quantization noise in transform coded images[END_REF], and explicit modeling of this type of noise is difficult. Thus, even top performing de-noising techniques such as the block-matching 3D algorithm (BM3D) [START_REF] Dabov | Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering[END_REF] have a limited efficiency in our context.

The proposed approach first consists in clustering patches extracted from the reconstructed sequence. This clustering step can thus be performed at both the encoder and the decoder side. At the encoder side, we then learn for each cluster a projection function between the degraded patches and the corresponding source patches from the input sequence. These projection functions are then transmitted to the decoder, where they can be applied for each cluster in order to enhance the degraded patches quality.

The rest of this paper is organized as follow. Section 2 reviews related work on image de-noising. Section 3 describes the proposed clustering-based quantization noise removal scheme. Experimental results are then discussed in Section 4.

RELATED WORK

Video codecs such as H.264 [START_REF] Wiegand | Overview of the H.264/AVC video coding standard[END_REF] and HEVC [START_REF] Sullivan | Overview of the high efficiency video coding (HEVC) standard[END_REF] use in-loop deblocking filters to address blocking artifacts [START_REF] Norkin | HEVC Deblocking Filter[END_REF], which are the most visible artifacts due to the block transform and quantization. Compared to previous video compression standards, HEVC also introduced a Sample Adaptive Offset filter in addition to the deblocking filter. HEVC is thus a challenging codec for noise removal applications, and for this reason it was chosen to perform the experiments (see Section 4).

However, these in-loop filters are designed with a minimal complexity, and thus have a limited efficiency. Note that a more efficient Adaptive Loop Filter was also considered in [START_REF] Chen | The adaptive loop filtering techniques in the HEVC standard[END_REF]. The coefficients of a local Wiener filter are learned at the encoder in order to minimize the error between the decoded and the original pixels, and transmitted to the decoder. While this method performs a local pixel-based filtering, the technique we propose exploits instead non-local patchbased similarities. In fact, image de-noising has been a very active research topic, and more complex and efficient methods have been proposed. The BM3D is a two step algorithm [START_REF] Dabov | Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering[END_REF] and among the best de-noising methods [START_REF] Chatterjee | Is denoising dead?[END_REF], as it combines the strengths of several previous de-noising approaches. It thus exploits non-local similarities as in the Non-Local Mean (NLM) algorithm [START_REF] Buades | A non-local algorithm for image denoising[END_REF], as well as filtering in the transform domain, such as in [START_REF] Donoho | De-noising by soft-thresholding[END_REF] [START_REF] Ghael | Improved wavelet denoising via empirical Wiener filtering[END_REF].

For the sake of conciseness, we briefly describe here the first step of the BM3D algorithm, which shows the main principles used in this algorithm. Each block of pixels in the noisy image is stacked into a 3D group along with its nearest neighbors, in the sense of the Euclidean distance. A 3D transform is then applied on each group, and the coefficients are filtered using hard thresholding. The denoised estimates are then obtained by applying the inverse 3D transform. Advantageously, all the blocks of the 3D group are processed at once. Several estimates are thus obtained for each pixel, which are eventually aggregated to obtain the final estimate. In the second step of BM3D, the estimate from the first step is used to perform Wiener filtering in the 3D transform domain.

More sophisticated noise models than AWGN have been proposed, such as the signal dependent noise model [START_REF] Kuan | Adaptive Noise Smoothing Filter for Images with Signal-Dependent Noise[END_REF]. However, the main approaches to remove such noise eventually rely on techniques designed for AWGN. Such techniques in fact assume that the noisy image can be partitioned into segments or clusters for which the noise can be modeled by AWGN [START_REF] Liu | Automatic estimation and removal of noise from a single image[END_REF] [START_REF] Liu | Signal dependent noise removal from a single image[END_REF].

Recent efficient de-noising methods also rely on sparse representations [START_REF] Elad | Image denoising via sparse and redundant representation over learned dictionaries[END_REF][21] [START_REF] Sulam | Image denoising through multi-scale learnt dictionaries[END_REF] or the local learning of dictionaries through clustering [START_REF] Chatterjee | Clustering-based denoising with locally learned dictionaries[END_REF] [START_REF] Dong | Sparsity-based image denoising via dictionary learning and structural clustering[END_REF]. Nonetheless, it was shown that these methods are usually outperformed by the BM3D algorithm [START_REF] Chatterjee | Is denoising dead?[END_REF], which we thus choose to perform some experiments in Section 4.

CLUSTERING-BASED QUANTIZATION NOISE REMOVAL FOR COMPRESSION

In this section we describe the proposed compression scheme improvement, which uses clustering-based linear mappings learning for noise removal.

Proposed scheme for compression

We denote M d and Ms the matrices containing in their columns all the vectorized M × N patches (i.e. pixel blocks) extracted from the decoded sequence Y and the source sequence X respectively. The main idea of the proposed scheme is represented in Fig. 1, and consists in the following steps:

• At the encoder side: decode the linear mappings cluster the decoded sequence M d , as it is performed at the encoder side for each cluster c, apply the corresponding linear mapping Pc to the decoded sequence patches M c d to obtain the de-noised patches M c r :

M c r = PcM c d (1)
We can see that the proposed method does not rely on any assumption on the noise model. We can thus apply this scheme to different codecs, and expect it to be robust to new designs of transform or quantizations which could affect the quantization noise.

Furthermore, we can combine our method to existing de-noising techniques, which only partially remove the noise in our context. We thus propose an alternative two step scheme, where a blind denoising algorithm is first performed on the decoded frames, before applying the proposed method (see Fig. 1). The matrix M d then contains in its column the vectorized de-noised patches. We qualify this first pass de-noising as "blind" because no parameters should be transmitted between the encoder and the decoder side, so that there is no bit-rate increase. In practice, some parameters that have a negligible rate cost are sent.

We discuss more in details the different steps of the approach in the following sections.

Clustering

The clustering is performed on M × N non-overlapping patches extracted from the decoded frames. Any clustering algorithm is eligible for the proposed method, as long as it can be performed similarly at the encoder and the decoder side. We choose the popular K-means algorithm [START_REF] Macqueen | Some Methods for classification and Analysis of Multivariate Observations[END_REF]. Note that some tests performed with more advanced clustering methods, such as spectral clustering [START_REF] Uw | On spectral clustering: Analysis and an algorithm[END_REF], showed that the de-noising performances were in the end very similar, despite the increased complexity.

Selecting the number of clusters K is crucial for the overall RD performances, as increasing K improves the de-noising performances, but also the bit-rate associated to the linear mappings. Empirical results show that 10 is a good order of magnitude for K (see Section 4). However, a more specific tuning of K depending on the video content and the target bit-rate can improve the final RD performances. To address this problem, we present below an adaptive method for selecting K, based on an rate-distortion optimization (RDO) criterion.

Instead of explicitly partitioning the data into K clusters, we propose a recursive binary partitioning, which proceeds as follow:

• At initialization, the full data set is considered as one cluster.

• Each cluster is then recursively split into 2 clusters if the RDO criterion (explained below) is satisfied, creating a binary tree structure.

• The procedure stops when no RDO criterion is satisfied, or a pre-defined maximum tree depth is reached.

The RDO criterion used to decide to further split a cluster is designed to balance the de-noising performances and the linear mappings rate cost. Given a cluster Cn and its corresponding mapping Pn, we can allocate it a RD cost computed as:

Jn = Dn + λRn (2)
The distortion Dn is computed as the Sum of Squared Error (SSE) between the patches of Cn de-noised with Pn and the corresponding source patches. The rate Rn is estimated as the number of bits of the encoded mapping Pn. The Lagrangian parameter is computed as in the test model of HEVC:

λ = QP f actor × 2 ((QP -12)/3) ( 3
)
where QP is the quantization parameter and QP f actor is an adjustment factor, set to further divided in two clusters Cn+1 and Cn+2, we merely check the condition Jn+1 + Jn+2 < Jn. Note that this method is independent from the clustering algorithm. An example of binary tree structure obtained after recursive partitioning is given in Fig. 2.

To perform the exact same partitioning at the decoder, we then need to transmit the binary tree structure. In our experiment, the maximum binary tree depth is set to 4, which corresponds to a maximum of 16 clusters and requires a maximum of 15 bits to be encoded. Thus, the bit-rate allocated to binary tree structure is negligible compared to the bit-rate of the encoded sequence.

Linear mappings learning and encoding

We first describe here the process to learn the linear mappings at the encoder, which are then used to perform de-noising at the decoder side (see Eq. 1). Considering multivariate linear regression, the problem is of searching for each cluster c the function Pc minimizing:

E = ||M c s T -M c d T P T c || 2 (4)
which is of the form ||Y -XB|| 2 (corresponding to the linear regression model Y = XB + E). The minimization of Eq. ( 4) gives the least squares estimator

Pc = M c s M c d T (M c d M c d T ) -1 (5)
The linear mappings are thus obtained in a matrix form, which then need to be transmitted to the decoder. For that purpose, we consider the K linear mapping matrices of size (M N ) × (M N ) as a sequence of K images, and encode them using a video codec. An example of such matrices is given in Fig. 3, where we can observe redundancies within and between the different matrices.

The linear mappings matrices consist originally in floating points values, thus they need to be quantized to integer values before encoding. To avoid a loss in terms of de-noising performance due to the quantization, the matrices are quantized on 16 bits integer values. The matrices are then encoded with the Range Extension (RExt) of HEVC [START_REF] Sullivan | Standardized extensions of high efficiency video coding (HEVC)[END_REF] [28], which supports 16 bits input. The quantization parameter of HEVC was chosen very low (QP = -30) to limit the loss on the de-noising performances. Finally, in order to recover the floating point matrices, the minimum and maximum values of the original floating points matrices need to be sent to the decoder without loss (with a 32 bits precision).

Oracle clustering for de-noising based on linear mappings learning

Finally, we describe in this section a so-called oracle clustering, which goal is to maximize the de-noising performances of the proposed scheme for a fixed number of clusters K. In other words, we want to minimizes the distortion at the decoder side while keeping a steady bit-rate associated to the linear mappings. For that purpose, the oracle clustering jointly learns the partitioning and the associated linear mappings in order to minimize the overall reconstruction error. Note that this clustering is denoted "oracle", as it relies on the knowledge of the source patches, and can not be directly use in our scheme as it is not reproducible at the decoder side. Nonetheless, it allows to evaluate an upper bound on the coding performances.

Formally, the problem corresponding to the oracle clustering is formulated as follows:

min C i=1...K x d ,xs∈C i ||xs -Pix d || 2 (6)
where xs and x d are corresponding source and decoded patches respectively, and C = {Ci}i=1...K . The problem of Eq. ( 6) can be solved using a greedy algorithm, described in Algorithm 1.

The initial set of linear mappings P

1 . . . P

K is supposed known, e.g. obtained by first performing the K-means clustering on the coded/decoded patches and learning a linear mapping for each cluster.

Algorithm 1 Oracle clustering for linear mapping learning

Input: cluster number K, coded/decoded patches matrix M d = {x d }, corresponding source patches matrix Ms = {xs} Output: cluster set C = {C i , 1 ≤ i ≤ K}, corresponding mapping matrices P = {P i , 1 ≤ i ≤ K} Repeat the two following steps until convergence:

Assignment step: ∀i = 1 . . . K, C (t) i = {x d , xs : ||xs -P (t) i x d || 2 ≤ ||xs -P (t) j x d || 2 ∀j, 1 ≤ j ≤ K} Update step: ∀i = 1 . . . K, P (t+1) i = M i s M i d T (M i d M i d T ) -1

SIMULATIONS AND RESULTS

The sequences used in the experiment are presented in Table 1, and mainly consist in HEVC test sequences. The test sequences are processed per Group Of Pictures (GOP), where a GOP contains a number of frames equal to the frame rate. The sequences are encoded with the HEVC test model HM (ver 15.0) [29] using the Main profile in Random Access, with 4 values for the Quantization Parameter, QP = 22, 27, 32, 37. Note that the proposed method is applied on the luminance channel. The patch size is set to M × N = 4 × 4, which has been observed as offering a good trade-off between denoising performances and rate cost for transmitting the linear projection matrices.

We present below rate-distortion performances for several experiments, computed using the Bjontegaard metric [START_REF] Bjontegaard | Calculation of average PSNR differences between RD-curves[END_REF] with respect to the sequences encoded with HEVC. More detailed results are available at http://www.irisa.fr/temics/ demos/clusteringLinearMappingNoiseRemoval/ clusteringLinearMappingNoiseRemoval.html. We show in Table 2 the RD performances of our method (without first pass de-noising) for each sequence, first performed with a fixed number of clusters K = 10 (left column) and then compared to the adaptive selection of K (center column). The results are averaged over all GOPs. We can see that almost all sequences benefit from the adaptive K. We show the selected value of K as well, averaged over the four QP values, which varies between about 5 and 16 depending on the sequence. The number of clusters is generally reaching higher values for sequences with a higher resolution and/or frame rate. The K value depends on the bit-rate as well, and we have in average over the sequences K = 12.3, 11.1, 9.0, 8.0 for QP = 22, 27, 32, 37 respectively.

In the right column, we also give the performances of the oracle clustering with K = 10 clusters, which show that the proposed method could potentially reach about 30% of bit-rate reduction in average. Note that these results merely represent an upper bound on the performances, as the oracle clustering can not be performed independently at the decoder side.

We then show in Table 3 the RD performances of our method, first using the one step scheme as above (left column), then compared to the two step scheme (right column), where the videos are first denoised with the extension of BM3D to videos (VMB3D) [START_REF] Dabov | Video denoising by sparse 3D transform-domain collaborative filtering[END_REF]. The noise standard deviation parameter is estimated as the square root of the MSE computed on a GOP at the encoder side, and is then sent to the decoder to perform the exact same de-noising. We use the adaptive K selection for both schemes. Results are given for the first GOP of the sequences.

We can see that a substantial bit-rate reduction can be achieved when combining the proposed method with the VBM3D. In fact, the VBM3D only partially removes the quantization noise, which is further reduced when using the proposed method in the second step. All sequences benefit from the two step schemes, and in particular the People On Street sequence, which reaches almost 10% of bit-rate saving.

We also give the selected K values, and we can see that the number of clusters selected is higher when using the first pass de-noising. In fact, in the RDO decision, the improvement brought by the first pass de-noising reduces the distortion, and thus allows for a higher bit-rate cost, i.e. more linear mappings can be sent to the decoder.

It has been observed that the processing time increase brought by the proposed method (with a non-optimized Matlab implementation) compared to the HM software is around 110% at the encoder and around 4000% at the decoder side. In the current version, more than 95% of the proposed method running time is dedicated to the K-means clustering. We thus expect the complexity increase to be much more reduced with an optimized implementation, e.g. consid- ering fast algorithms [START_REF] Hamerly | Making k-means even faster[END_REF], or parallel implementation [START_REF] Zhang | The Study of Parallel K-Means Algorithm[END_REF], possibly on GPU [START_REF] Giuroiu | CUDA K-Means Clustering[END_REF].

CONCLUSION AND PERSPECTIVES

In this paper we introduced a novel generic scheme designed to de-noise compressed videos and improve the coding performances of existing codecs. The compressed video patches are first clustered, and linear mappings are then learned between the decoded and source patches at the encoder side. These linear mappings are then transmitted to the decoder to perform de-noising. We also propose an algorithm to optimally select the number of clusters in a RD sense. Experimental results thus show that our method can improve the coding performances of the HEVC standard. In addition, the proposed technique is also complementary to existing de-noising approaches, and by combining the different techniques, we can reach even higher bit-rate reduction, up to about 10%. Furthermore, we showed through the oracle clustering that much higher coding gains could be reached, which calls for the investigation of more specific clustering techniques instead of the K-means algorithm. Other perspectives could be explored as well thanks to the genericity of the proposed scheme. For example, an extension to color videos is quite straightforward, by vectorizing the three components of the patches instead of the sole luminance. Moreover, extension to scalable applications can be considered, such as spatial scalability, for which the proposed scheme could be used to perform joint de-noising and super-resolution. In fact, a similar method has been proposed recently for super-resolution applications in [START_REF] Zhang | Learning Multiple Linear Mappings for Efficient Single Image Super-Resolution[END_REF]. Finally, LDR/HDR or color gamut scalabilities could be also investigated [START_REF] Korshunov | Context-dependent JPEG backward-compatible high-dynamic range image compression[END_REF].

-

  cluster the decoded patches M d for each cluster c, learn a linear mapping Pc between the decoded patches M c d and the corresponding source patches M c s by minimizing ||M c s T -M c d T P T c || 2 encode the corresponding linear mappings (in matrix form) and transmit them to the decoder • At the decoder side:

Fig. 1 .

 1 Fig. 1. Proposed compression scheme using noise removal based on clustering and linear mappings learning.

Fig. 2 .

 2 Fig. 2. Example of a cluster partition based on a binary tree structure and an RDO criterion. The clusters are recursively split to reach 4 clusters. The tree structure is described by only 5 bits.

Fig. 3 .

 3 Fig. 3. Linear mapping matrices obtained on the first frame of the Kimono sequence encoded with HEVC at QP=22.

Table 1 .

 1 Test sequences

	Sequence	Resolution	Frames count	Frame rate
	City	1280 × 720	600	60
	Park Scene	1920 × 1080	240	24
	Tennis	1920 × 1080	240	24
	Kimono	1920 × 1080	240	24
	Cactus	1920 × 1080	500	50
	Terrace	1920 × 1080	600	60
	Basket	1920 × 1080	500	50
	Ducks	1920 × 1080	500	50
	People On Street	2560 × 1600	150	30
	Traffic	2560 × 1600	150	30

Table 2 .

 2 RD performances of the proposed scheme averaged over all GOPs (Bjontegaard bit-rate gain with respect to HEVC)

			K-means		Oracle clustering
		K = 10	Adaptive K	K = 10
		BD-rate	BD-rate K value	BD-rate
	City	-1.81	-2.24	8.9	-28.92
	Park Scene	0.38	-0.47	5.8	-21.66
	Tennis	0.08	-1.42	4.9	-38.20
	Kimono	0.43	-0.94	5.8	-33.60
	Cactus	-1.08	-1.15	10.7	-23.48
	Terrace	-8.28	-8.08	14.4	-39.50
	Basket	-0.65	-0.80	8.8	-28.99
	Ducks	-1.88	-1.96	15.8	-22.07
	People On Street	-2.56	-2.60	16.0	-23.26
	Traffic	-1.26	-1.43	10.0	-25.50
	Average	-1.66	-2.11	10.1	-28.52

Table 3 .

 3 RD performances of the proposed schemes for the first GOP (Bjontegaard bit-rate gain with respect to HEVC)

		Adaptive K	Two step scheme
		BD-rate	K value	BD-rate	K value
	City	-2.76	8.5	-3.24	12.3
	Park Scene	-0.69	5.8	-1.07	7.5
	Tennis	-0.90	6.0	-7.77	6.8
	Kimono	-1.32	5.8	-4.32	5.8
	Cactus	-1.12	12.0	-6.18	13.0
	Terrace	-5.25	15.8	-7.24	15.8
	Basket	-1.10	7.5	-3.80	9.0
	Ducks	-2.39	15.3	-4.13	16.0
	People On Street	-2.57	16.0	-9.63	15.5
	Traffic	-1.44	10.3	-5.43	11.5
	Average	-1.96	10.3	-5.28	11.3


We read with interest the article by Silverstein and colleagues [START_REF] Silverstein | P3b, consciousness, and complex unconscious processing[END_REF] who questioned the putative specificity of the P3b event-related potentials (ERP) component as a neural signature of conscious access to a visual representation. Prior to this new study, numerous empirical reports revealed that a brain response peaking ~300ms after stimulus onset and maximally distributed over parietal electrodes -the so called P3b -is closely related to subjective visibility [START_REF] Vogel | Electrophysiological evidence for a postperceptual locus of suppression during the attentional blink[END_REF][START_REF] Sergent | Timing of the brain events underlying access to consciousness during the attentional blink[END_REF]. These experimental findings provided the bases to develop neuronal and computational theories of consciousness such as the global workspace model [START_REF] Dehaene | Towards a cognitive neuroscience of consciousness: basic evidence and a workspace framework[END_REF][START_REF] Dehaene | Conscious, preconscious, and subliminal processing: a testable taxonomy[END_REF][START_REF] Dehaene | Experimental and theoretical approaches to conscious processing[END_REF]. Silverstein and colleagues used a 'passive attentive' version of a masked visual odd-ball paradigm while recording scalp ERPs. On each trial, subjects were presented with either the masked word 'LEFT' (in 80% or 20% of trials) or the masked word 'RIGHT' (in 20% or 80% of trials). Word frequency was balanced across subjects, who were asked to attend carefully to the masked sequence. Not only were they instructed that this sequence contained a masked word, but also that: "however implausible it might seem, our prior data suggested that the stimuli would nonetheless be unconsciously perceived and produce brain wave effects-but only if they maintained their attention". When contrasting ERPs elicited by rare and frequent masked words, Silverstein and colleagues identified a P3b ERP component followed by a late, and sustained, slow wave (LSW). Given that subjects subjectively reported the absence of conscious perception of words, and that they performed at chance-level in a stimulus detection task performed after the main experiment, Silverstein and colleagues concluded that a P3b can be observed during unconscious perception. If valid, their interpretation would then simply invalidate the P3b as a possible candidate neural signature of conscious access.

This original and provocative study, however, raises both methodological and conceptual concerns which need to be addressed before one can adopt Silverstein and colleagues interpretation.

A set of methodological problems

The P3b is part of a larger complex of positive deflections -the so-called P300. Of particular importance here, the P3a can be functionally distinguished from the P3b: it is known to occur in the absence of conscious perception (Muller-Gass, Macdonald, Schroger, Sculthorpe and Campbell 2007) and even in non-conscious patients [START_REF] Faugeras | Event related potentials elicited by violations of auditory regularities in patients with impaired consciousness[END_REF]). P3a and P3b are both positive deflections and occur in similar time windows, but they can nevertheless be separated based on their topographies. The P3b is maximally distributed over parietal electrodes, while the P3a is more frontally distributed. The spatial sampling of the EEG signal is therefore critical to separate these ERPs. Surprisingly, the authors only used 3 midline (Fz, Cz, Pz) electrodes referenced to linked ears, as well as 2 electrodes at the right eye to detect eye movement artifacts.

As expected for a P3b component, the effect reported by Silverstein and colleagues was maximal over Pz, but we simply don't have access to scalp topographies of the ERP effects reported in this work. We agree that there is no intrinsic relation between the number of electrodes and the quality of a result, but in the context of distinguishing P3b from P3a ERP components this limitation turns into a genuine problem.

In the same vein, one of the most reliable findings in the vast odd-ball literature, is the existence of a N2 and mismatch negativity (MMN) ERP effect occurring before the P3 complex [START_REF] Tiitinen | Attentive novelty detection in humans is governed by pre-attentive sensory memory[END_REF]. The apparent absence of such an effect (a small inverse difference is seen on Figure 2) confirms the necessity of sampling brain activity with a richer spatial resolution in order to describe reliably the observed effects.

Moreover, shortcomings in statistical analyses of the ERPs deserve further discussion. The visual inspection of 'effects' suggest that the effect size reported by Silverstein et al are not quite different from fluctuations within the baseline and from other periods of ERPs (see e.g. Figure 2C). Actually, the authors did not assess significant differences on the entire time course of the ERP but only on predefined time windows. Thus it is impossible to determine whether the reported effects are temporally and spatially precise and specific to the P3b. A better approach would consist in performing non selective sample-by-sample tests, and then identifying temporal clusters during which ERPs significantly differed.

More importantly, although the article by Silverstein et al. opens by asking the fundamental question "How can perceptual awareness be indexed in humans?", their experimental design is lacking the crucial comparison of unconscious 'P3b' with its conscious equivalent. Rather than using exclusively masked trials, the authors could have added unmasked trials, in order to compare the properties (latency, amplitude and effect size, duration, topography) of conscious and unconscious ERP effects. By doing so previous studies could identify specific components of conscious access [START_REF] Dehaene | Conscious, preconscious, and subliminal processing: a testable taxonomy[END_REF]. From a theoretical perspective, we previously mentioned and modeled the possibility for a masked stimulus to "evoke transient workspace activity of variable intensity and duration" (see also Figure 1 in [START_REF] Dehaene | Towards a cognitive neuroscience of consciousness: basic evidence and a workspace framework[END_REF]. Such transient and partial activation of the workspace could appear as brief and small patterns of activity distinct from a large and sustained P3b component. Therefore, without this crucial conscious contrast, it becomes almost impossible to precisely qualify the observed ERP effect.

It is noteworthy that according to our theory, conscious access associated with the P3b is also associated with other signatures [START_REF] Gaillard | Converging intracranial markers of conscious access[END_REF]) such as: long-range synchrony in theta-alpha-beta band, decrease of alpha power, and late increase of gamma power. None of these neural signatures, complementary to the P3b, are tested here and the nature of the observed ERP effects therefore remains unclear.

Additionally, the interesting use by Silverstein and al. of 'LEFT' and 'RIGHT' words opened the possibility of complementing the results by lateralized readiness potentials (LRPs) analyses. Such analyses proved to be very useful to explore both unconscious and conscious processing of masked primes (Dehaene, Naccache, Le Clec, Koechlin, Mueller, [START_REF] Dehaene | Imaging unconscious semantic priming[END_REF][START_REF] Eimer | Effects of masked stimuli on motor activation: behavioral and electrophysiological evidence[END_REF]. Unfortunately, the use of only 3 midline electrodes, and the absence C3/C4 electrodes precluded this interesting complementary approach.

Furthermore, from a Bayesian perspective, we consider the authors should have mentioned and discussed more extensively the large set of empirical evidence that their finding seems to contradict: numerous studies conducted in normal controls as well as in many clinical settings (e.g.: blindsight, visual neglect) support the P3b theory by reporting rich unconscious processing of visual stimuli without any late P3b signature (for a review see [START_REF] Dehaene | Experimental and theoretical approaches to conscious processing[END_REF]). This literature, acting here as a strong prior against Silverstein and colleagues interpretation, needs to be addressed.

Conscious metacognition of unconscious perceptual processes?

Beyond these notable methodological issues, this article also raises a more profound question. The major difference between this study and previous studies rests in the fact that subjects were told from the very beginning of the presence of masked stimuli, and were instructed to pay attention to them very carefully. Therefore, even if we discard the methodological issues we just raised, and consider that these results are correct, it may be the case that the P3b signature observed here between deviant and standard stimuli would correspond to a metacognitive effect, that is to say to a conscious access to the consequence of unconscious processing of masked primes. For instance, a motor effect induced by the processing of the rare 'LEFT' prime (or 'RIGHT' for other subjects) in motor areas may well lead to conscious access to a subjective confidence information that the prime was deviant or standard. By amplifying subjects' attention to monitor prime processing, this metacognitive interpretation may well explain the striking pattern of results reported here. Interestingly, a growing empirical evidence demonstrates that a large class of unconscious cognitive processes are strongly influenced by the conscious posture and endogenous attentional allocation [START_REF] Naccache | Unconscious masked priming depends on temporal attention[END_REF]. In addition to such an amplification, it might be the case that subject informed of the presence of subliminal stimuli could more easily introspect a form of surprise originating either from perceptual or from motor-related areas ('LEFT', 'RIGHT'). In other words, this study may illustrate conscious access to the downstream effects of an unconsciously perceived stimulus. Interestingly, a recent study using a visual masked priming paradigm reported that the conflict between masked prime and visible target stimuli modulated two ERP components (Desender, Van Opstal, Hughes and Van den Bussche in press): an early N2 component, as well as a late P3 complex. During this experiment, subjects had to perform two tasks on each trial: they first had to respond to the target, and then to introspect on the difficulty of the trial. Nicely, introspection of the prime-target conflict elicited by the unconscious processing of the prime was possible, and correlated only with the P3 component.

Similarly, in the study by Silverstein and colleague, one may suppose that the P3b component and the LSW they observed corresponds to a conscious introspection of processes elicited by the unconsciously perceived prime.

As a conclusion, if the results reported in Silverstein et al do correspond to a genuine P3b ERP component (but see our methodological concerns above), they may elegantly illustrate the complex relations prevailing between conscious and unconscious processes, and still not refute the relation prevailing between conscious access and the P3b ERP component.

  -HP, Groupe hospitalier Pitié-Salpêtrière, Department of Neurology, Paris, France 2-AP-HP, Groupe hospitalier Pitié-Salpêtrière, Department of Neurophysiology, Paris, France 3-INSERM, U 1127, F-75013, Paris, France 4-Institut du Cerveau et de la Moelle épinière, ICM, PICNIC Lab, F-75013, Paris, France 5-Sorbonne Universités, UPMC Univ Paris 06, Faculté de Médecine Pitié-Salpêtrière, Paris, France 6-Cognitive Neuroimaging Unit, CEA DSV/I2BM, INSERM, Université Paris-Saclay, NeuroSpin Center, 91191 Gif/Yvette, France 7-Ecole des Neurosciences de Paris Ile-de-France, 75006 Paris, France

Funding

This work has been supported by the Fondation pour la Recherche Médicale ('Equipe FRM 2015' grant to L.N. and by the 'Recovery of consciousness after severe brain injury Phase II' grant of the James S. McDonnell Foundation.