Vít Dolejší
email: dolejsi@karlin.mff.cuni.cz

Alexandre Ern
email: ern@cermics.enpc.fr

Martin Vohralík
email: vohralik@ann.jussieu.fr

A framework for robust a posteriori error control in unsteady nonlinear advection-diffusion problems *

Keywords: unsteady nonlinear advection-diffusion problem, a posteriori estimate, dual norm, flux reconstruction, flux equilibration, unified framework, robustness, discontinuous Galerkin method

We derive a framework for a posteriori error estimates in unsteady, nonlinear, possibly degenerate, advection-diffusion problems. Our estimators are based on a space-time equilibrated flux reconstruction and are locally computable. They are derived for the error measured in a space-time mesh-dependent dual norm stemming from the problem and meshes at hand augmented by a jump seminorm measuring possible nonconformities in space. Owing to this choice, a guaranteed and globally efficient upper bound is achieved, as well as robustness with respect to nonlinearities, advection dominance, domain size, final time, and absolute and relative size of space and time steps. Local-in-time and in-space efficiency is also shown for a localized upper bound of the error measure. In order to apply the framework to a given numerical method, two simple conditions, local space-time mass conservation and an approximation property of the reconstructed fluxes, need to be verified. We show how to do this for the interiorpenalty discontinuous Galerkin method in space and the Crank-Nicolson scheme in time. Numerical experiments illustrate the theory.

Introduction

We consider the unsteady nonlinear problem

∂ t u -∇•σ(u, ∇u) = f in Q := Ω × (0, t F ), (1.1a)
u = 0 on ∂Ω × (0, t F ), (1.1b) u(•, 0) = u 0 in Ω, (1.1c) with Ω ⊂ R d , d ≥ 2, a polygonal (polyhedral) domain, t F > 0 the final time, f the source term, and u 0 the initial datum. The function σ(u, ∇u) takes the form

σ(u, ∇u) := K(u)∇u -φ(u), (1.2)
where K(•) is a nonlinear, possibly degenerate, tensor-valued function associated with diffusive transport and φ(•) a nonlinear vector-valued function associated with advective transport. Precise assumptions on K(•) and φ(•) are specified in Section 2. [START_REF] Adams | FISHPACK -efficient FOR-TRAN subprograms for the solution of separable elliptic partial differential equations[END_REF]. In what follows, u is termed the potential and -σ(u, ∇u) the advective-diffusive flux. We assume that the problem (1.1a)-(1.1c) admits a unique weak solution u and that a fully discrete approximate potential, say u hτ , is available. The space discretization scheme can be nonconforming, so that u hτ exhibits jumps across the spatial mesh interfaces. Quite a broad literature has been devoted to the a posteriori error analysis for problems of the form (1.1a)-(1.1c). One of the main issues is to simultaneously prove global upper and (local) lower bounds for the error (reliability and (local) efficiency). For the linear heat equation, reliability is achieved in the energy norm by Picasso [START_REF] Picasso | Adaptive finite elements for a linear parabolic problem[END_REF] and Repin [START_REF] Repin | Estimates of deviations from exact solutions of initial-boundary value problem for the heat equation[END_REF], and in higher-order norms by Makridakis and Nochetto [START_REF] Makridakis | Elliptic reconstruction and a posteriori error estimates for parabolic problems[END_REF]. In Verfürth [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF] and Bergam et al. [START_REF] Bergam | A posteriori analysis of the finite element discretization of some parabolic equations[END_REF], both reliability and efficiency are proven by augmenting the energy norm by a dual norm of the time derivative. The lower bound is then local in time but global in space. A further analysis unifying various space discretization schemes has been recently given in [START_REF]A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF]. For nonlinear parabolic problems, reliability and local and global efficiency have been derived by Verfürth [START_REF] Verfürth | A posteriori error estimates for nonlinear problems. L r (0, T ; L ρ (Ω))-error estimates for finite element discretizations of parabolic equations[END_REF][START_REF]A posteriori error estimates for nonlinear problems: L r (0, T ; W 1,ρ (Ω))-error estimates for finite element discretizations of parabolic equations[END_REF] under a restriction on the relative size of space and time steps. This restriction has been lifted by Verfürth in [START_REF]A posteriori error estimates for non-linear parabolic equations[END_REF] for convex two-dimensional spatial domains thereby achieving global efficiency; the price to pay is a solution of a linear diffusion problem by the finite element method on each time step. Robustness with respect to advection dominance or nonlinearities is not addressed in [START_REF]A posteriori error estimates for non-linear parabolic equations[END_REF]. For unsteady linear advection-diffusion problems, Verfürth [START_REF]Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF] has proved robustness with respect to advection dominance while augmenting the energy norm by a dual norm of the material derivative. A solution of a finite element reaction-diffusion problem on each time step is, again, necessary. All the above works concern the nondegenerate case; still less results are available for degenerate parabolic problems. We cite in particular Nochetto et al. [START_REF] Nochetto | A posteriori error estimation and adaptivity for degenerate parabolic problems[END_REF], where only the error upper bound is derived.

The purpose of the present work is to derive guaranteed, (locally) efficient, and robust a posteriori error estimates for the problem (1.1a)- (1.1c). Here, guaranteed means that the estimates represent an explicitly computable upper bound on the error and robustness means that the efficiency is independent of the main model and discretization parameters, that is, nonlinearities, advection dominance, domain size, final time, and absolute and relative size of space and time steps. While these results are of independent theoretical interest, their practical motivation is on the one hand to certify the simulation error and on the other hand to provide criteria for space-time adaptive mesh refinement in view of efficient numerical computations.

Our key idea is the introduction of a space-time mesh-dependent dual norm to measure the error, see (2.7). This norm includes the full space-time nonlinear advection-diffusion operator, which stands in contrast to previous work where only parts of the differential operator at hand are included (i.e., the time derivative or the material derivative as in [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF]Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF]). Such approaches have been used recently by Chaillou and Suri [START_REF] Chaillou | Computable error estimators for the approximation of nonlinear problems by linearized models[END_REF] and in [START_REF] Alaoui | Guaranteed and robust a posteriori error estimates and balancing discretization and linearization errors for monotone nonlinear problems[END_REF] in the context of steady nonlinear diffusion problems. Moreover, evaluating the dual norm with respect to a specific mesh-dependent norm for test functions with bounded time and space derivatives in L 2 (Q), see (2.5a)-(2.5b), also allows one to simplify substantially the proof of the error lower bound, whereby space-time bubble functions can be considered instead of the more usual space bubble functions at fixed times. This point also presents the practically crucial advantage that an error lower bound can be achieved using locally computable estimators, in contrast to [START_REF]A posteriori error estimates for non-linear parabolic equations[END_REF][START_REF]Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF] where the solution of a global diffusion/reactiondiffusion problem on each time step is necessary to evaluate the estimators. The error measure (2.7) may seem rather weak at a first glance. However, this measure admits an easily and locally computable upper bound which consists of weighted L 2 -norms of the potential error and of the error in the nonlinear advectiondiffusion flux, see (2.10)- (2.11). Our efficiency results carry over to this norm and, moreover, can then be localized in space and in time. The numerical experiments of Section 8 actually show that our estimators generally provide efficient estimates for this weighted L 2 -norm.

Our results are derived in a unified framework where the actual numerical scheme used to obtain u hτ need not be specified. The error upper bound hinges on an advection-diffusion flux reconstruction and its local space-time equilibration, see Assumption 3.1, while the error lower bound requires a local approximation property on this flux, see Assumption 4.1. Applying the present framework to a given numerical scheme simply boils down to verifying these two assumptions. The reconstruction of the flux depends on the space discretization scheme, as discussed in [START_REF]A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF] for the linear heat equation. Our approach is thus in line with the developments relying for linear model problems on the Prager and Synge equality [START_REF] Prager | Approximations in elasticity based on the concept of function space[END_REF] and pursued later by, e.g., Ladevèze [START_REF] Ladevèze | Comparaison de modèles de milieux continus[END_REF], Luce and Wohlmuth [START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF], Braess et al. [START_REF] Braess | Equilibrated residual error estimates are p-robust[END_REF], or Ainsworth [START_REF] Ainsworth | A framework for obtaining guaranteed error bounds for finite element approximations[END_REF].

This paper is structured as follows. We present the continuous and discrete settings and define our error measure J u (u hτ ) in Section 2. This is a sum of the above-discussed space-time mesh-dependent dual norm J u,FR (u hτ ), see (2.7), and of a weighted jump term J u,NC (u hτ ) which measures possible nonconformities in space, see (2.12). We state our a posteriori error estimate, Theorem 3.3, in Section 3 yielding the error upper bound J u (u hτ ) ≤ η FR + η NC + η IC . Theorem 4.4 of Section 4 then establishes the error lower bound η FR + η NC J u (u hτ ) under the assumption that the quadrature error caused by the nonlinearity of the flux σ is small enough. Here, means up to a generic constant which is independent of the main model and discretization parameters as specified above (the possible dependencies of the constant are stated in Section 4). Moreover, Theorem 4.2 provides a space-time localized version in terms of the computable error upper bound (2.10)-(2.11) on the error measure J u,FR (u hτ ). We devote Sections 5 and 6 to the proofs of these results. To illustrate the developed abstract framework, we apply it in Section 7 to the interiorpenalty discontinuous Galerkin (IPDG) method in space and the Crank-Nicolson scheme in time. Numerical experiments, including nonlinear degenerate advection-diffusion problems, are presented in Section 8.

The setting

This section briefly describes the continuous and discrete settings and discusses the error measure.

Continuous setting and weak solution

We consider the function spaces

X := L 2 (0, t F ; H 1 0 (Ω)), (2.1a)
Y := {ϕ ∈ L 2 (0, t F ; H 1 0 (Ω)); ∂ t ϕ ∈ L 2 (Q); ϕ(•, t F ) = 0}, (2.1b) recalling that ϕ ∈ L 2 (0, t F ; H 1 0 (Ω)) and ∂ t ϕ ∈ L 2 (0, t F ; H -1 (Ω)) classically imply ϕ ∈ C 0 ([0, t F ]; L 2 (Ω)
). The weak solution u of (1.1a)-(1.1c) is sought in the space X, whereas the space Y is used as test space. Specifically, we assume that there exists a unique weak solution u ∈ X such that

tF 0 {(f, ϕ) + (u, ∂ t ϕ) -(σ(u, ∇u), ∇ϕ)}(t) dt + (u 0 , ϕ(•, 0)) = 0 ∀ϕ ∈ Y, (2.2)
where (•, •) denotes the inner product in L 2 (Ω) or [L 2 (Ω)] d . In order to ensure that all the terms in (2.2) are well-defined, we assume

f ∈ L 2 (Q), u 0 ∈ L 2 (Ω), and σ(u, ∇u) ∈ [L 2 (Q)] d which is satisfied, for u ∈ X ∩ L ∞ (Q), e.g., if K ∈ L ∞ loc (R; R d×d ) and φ ∈ C 1 (R; R d ).
In deriving our a posteriori error estimators in Section 3, we exploit the fact that (2.2) consists, except for the contribution of the initial condition, of space-time inner products in L 2 (Q).

The present framework also covers some particular cases of degenerate parabolic equations of the form

∂ t b(v) -∇•( K∇v -φ(v)) = f, (2.3)
where b(•) is an increasing function with locally Hölder regularity such that b(0) = 0. Problem (2.3) includes slow-diffusion-type problems, e.g., the porous media equation for which b(v) = v 1/κ , 1 < κ < +∞ (so that b ′ (0) = +∞), and fast-diffusion-type problems of elliptic-parabolic form, e.g., the Richards equation for which typically b ′ (0) = 0. Existence, uniqueness, and regularity results for problem (2.3) can be found in the work of Alt and Luckhaus [START_REF] Alt | Quasilinear elliptic-parabolic differential equations[END_REF] and Otto [START_REF] Otto | L 1 -contraction and uniqueness for quasilinear elliptic-parabolic equations[END_REF], leading to weak solutions

v ∈ X ∩ L ∞ (Q). Since the function b is increasing, equation (2.3) can be recast into the form (1.1a) by setting u := b(v), yielding K(u) := K(b -1 ) ′ (u) and φ(u) := φ(b -1 (u))
. Therefore, the present analysis can be applied to (2.3) under the assumption b(v) ∈ X. This assumption holds in the fast-diffusion regime, but not necessarily in the slow-diffusion regime where it is possible that v ∈ X but b(v) ∈ X. The Stefan problem, governed by ∂ t u -∆β(u) = 0 where β(•) is a nondecreasing Lipschitz function, is not covered by the present assumptions either, since, in this case, the weak solution u can exhibit jumps. Finally, we observe that we do not include zero-order terms in (1.1a).

Remark 2.1 (Physical units). In advection-diffusion problems, the physical unit of the components of the diffusion tensor K(•) is L 2 T -1 and that of the components of the transport velocity φ ′ (•) is LT -1 . Here, L stands for length and T for time.

Discrete setting and approximate solution

We consider an increasing sequence of discrete times {t n } 0≤n≤N such that t 0 = 0 and t N = t F . We set I n := (t n-1 , t n ] and τ n := t n -t n-1 for all 1 ≤ n ≤ N . We consider a time-sequence of matching simplicial meshes {T n } 0≤n≤N of the spatial domain Ω. The mesh T n is used to approximate the solution at the discrete time t n . We assume that T n is obtained from T n-1 by refining some elements and coarsening some other ones. For all 1 ≤ n ≤ N , we denote by T n-1,n the coarsest common refinement of T n-1 and T n and by T n-1,n the finest common coarsening of T n-1 and T n . Obviously, if the meshes are kept fixed, T n-1,n = T n-1,n = T n = T 0 . In the general case, the mesh T n-1,n is needed to handle nonconforming functions in space that are continuous in time, and the mesh T n-1,n is used to facilitate the flux equilibration. At this stage, we do not need to specify any numerical scheme with which to obtain the approximate solution u hτ . We simply assume that u hτ is in the space

X h := {ϕ ∈ L 2 (0, t F ; H 1 (T )); ∂ t ϕ ∈ L 2 (Q)}.
Here, L 2 (0, t F ; H 1 (T )) is spanned by those functions ϕ ∈ L 2 (Q) that, for all 1 ≤ n ≤ N and all T ∈ T n-1,n , satisfy ϕ| T ×In ∈ L 2 (I n ; H 1 (T )). Functions in X h can exhibit jumps across the interfaces of T n-1,n , but, owing to the assumption ∂ t ϕ ∈ L 2 (Q), they are continuous with respect to time. This latter assumption is needed to express the residual in terms of space-time inner products in L 2 (Q). Let 1 ≤ n ≤ N and let T ∈ T n or T ∈ T n-1,n or T ∈ T n-1,n . We use the notation h T for the diameter of T and • T for the norm in L 2 (T ). Similarly, the norm in

L 2 (T × I n ) is denoted by • T ×In .
The corresponding inner products are denoted by (•, •) T and (•, •) T ×In , respectively. We collect in F T all the faces of T and in F int T those that are subsets of Ω. The set of all faces of the mesh T n is denoted by F n , a generic mesh face by F , and we use a similar notation for norms and inner products as above. For an interface F , n F denotes its unit normal oriented in the sense the jump is evaluated, while, for a mesh element T , n T denotes its outward unit normal. Similarly, T T stands for all mesh elements sharing a face with the element T , whereas T F denotes the mesh elements sharing the face F . Finally, we denote by ∇ the broken space gradient.

The following inequality is instrumental in deriving our error upper bound.

Lemma 2.2 (Scaled space-time Poincaré inequality). Let 1 ≤ n ≤ N and let T ∈ T n-1,n . Let Π 0 denote the L 2 (Q)-orthogonal projection onto constants in each space-time element T × I n . Set C P := 1/π. Then, for all ϕ ∈ H 1 (T × I n ),

ϕ -Π 0 ϕ T ×In ≤ C P h 2 T ∇ϕ 2 T ×In + (τ n ) 2 ∂ t ϕ 2 T ×In 1 2 .
(2.4)

Proof. Let ϕ ∈ H 1 (T × I n ) and, for all t ∈ I n , set φ(t) := |T | -1 T ϕ(x, t) dx. Observing that (ϕ -φ) and ( φ -Π 0 ϕ) are L 2 (T × I n )-orthogonal, we infer ϕ -Π 0 ϕ 2 T ×In = ϕ -φ 2 T ×In + φ -Π 0 ϕ 2 T ×In .
For the first term on the right-hand side, the usual Poincaré inequality on T (which is convex) yields

ϕ -φ 2 T ×In ≤ In C 2 P h 2 T T |∇ϕ| 2 (x, t) dx dt = C 2 P h 2 T ∇ϕ 2 T ×In .
For the second term, observing that Π 0 ϕ = (τ n ) -1

In φ(t) dt, the one-dimensional Poincaré inequality on I n and the Cauchy-Schwarz inequality yield

φ -Π 0 ϕ 2 T ×In ≤ |T | -1 π -2 (τ n ) 2 In T ∂ t ϕ(x, t) dx 2 dt ≤ π -2 (τ n ) 2 ∂ t ϕ 2 T ×In .
Collecting the above bounds yields (2.4).

Error measure

The error measure described in this section combines a space-time mesh-dependent dual norm plus a nonconformity term.

Space-time mesh-dependent dual norm

Recalling definition (2.1b) of the space Y , we equip it with the norm

ϕ 2 Y,T ×In := C T,n h 2 T ∇ϕ 2 T ×In +(τ n ) 2 ∂ t ϕ 2 T ×In , ∀1 ≤ n ≤ N, T ∈ T n-1,n , (2.5a)
ϕ 2 Y := N n=1 T ∈T n-1,n ϕ 2 Y,T ×In . (2.5b)
• Y is indeed a norm since functions in Y vanish on Ω × {t F } and on ∂Ω × (0, t F ). We observe that this norm depends on the space-time meshes. The positive quantities C T,n are user-dependent weights, typically with physical unit T -1 , so that the error measure has the same unit as the classical energy norm. A possible example is

C T,n := (τ n ) -2 t F + h -2 T C φ,T,n + h -2 T C K,T,n , (2.6)
with a weight C φ,T,n proportional to a length scale times an advective velocity scale and a weight C K,T,n proportional to a diffusion scale. Since our estimators derived in Section 3 are also scaled by these weights, our robustness results are not influenced by the value assigned to these weights, see Remark 4.5.

The first building block of our error measure is the quantity

J u,FR (u hτ ) := sup ϕ∈Y, ϕ Y =1 tF 0 {(u hτ -u, ∂ t ϕ) + (σ(u, ∇u) -σ(u hτ , ∇u hτ ), ∇ϕ)}(t) dt. (2.7) Let v ∈ L 2 (0, t F ; H 1 (T )). Define the residual R(v) ∈ Y ′ such that, for all ϕ ∈ Y , R(v), ϕ Y ′ ,Y := tF 0 {(f, ϕ) + (v, ∂ t ϕ) -(σ(v, ∇v), ∇ϕ)}(t) dt + (u 0 , ϕ(•, 0)). (2.8)
Then we infer, owing to (2.2), that

J u,FR (u hτ ) = sup ϕ∈Y, ϕ Y =1 R(u hτ ), ϕ Y ′ ,Y , (2.9)
showing that J u,FR (u hτ ) is a dual norm of the residual of the approximate solution.

The error measure J u,FR (u hτ ) cannot be computed easily in practice (in the test cases where the exact solution u is available). Indeed, its evaluation requires solving the following (infinite-dimensional spacetime) problem: Find

ψ ∈ Y such that (ψ, ϕ) Y = R(u hτ ), ϕ Y ′ ,Y for all ϕ ∈ Y , where (•, •) Y denotes the inner product corresponding to the • Y -norm. Then, it is immediate that J u,FR (u hτ ) = ψ Y .
However, a computable upper bound on J u,FR (u hτ ) can be readily derived using the Cauchy-Schwarz inequality, leading to

J u,FR (u hτ ) ≤ sup ϕ∈Y, ϕ Y =1 N n=1 T ∈T n-1,n e n FR,T ϕ Y,T ×In ≤    N n=1 T ∈T n-1,n (e n FR,T ) 2    1 2
(2.10) that we denote by e FR , with, for all 1 ≤ n ≤ N and all T ∈ T n-1,n ,

e n FR,T := C -1 2 T,n (τ n ) -2 u hτ -u 2 T ×In + h -2 T σ(u, ∇u) -σ(u hτ , ∇u hτ ) 2 T ×In 1 2 .
(2.11)

Remark 2.3 (Dual norm of the residual and energy-type norms). Under appropriate assumptions, a functional framework can be introduced for the nonlinear differential operator in (1.1a), and the error between u and u hτ can be measured in the corresponding energy-type norms. For conforming approximations (u hτ ∈ X), following Verfürth (see, e.g., [30, Proposition 2.1]), the energy error can be bounded from above and from below by a dual norm of the residual. Such equivalence results hinge on suitable a priori bounds of the linearized differential operator, where it is in particular difficult to trace the influence of the size of the nonlinearities or of the advection dominance. Here we measure the error directly by (2.7) and avoid energy-type norms and the question of their equivalence with a dual norm of the residual.

Nonconformity

As we allow for nonconformities (X h is not a subspace of X), we need to introduce a second building block of our error measure,

J u,NC (u hτ ) :=    N n=1 T ∈T n-1,n F ∈FT C -1 T,n h -2 T C K,φ,T,F,n [[u -u hτ ]] 2 F ×In    1 2 , (2.12)
where [[•]] means the jump across interfaces and the actual value at boundary faces and, for

T ∈ T n-1,n , C T,n := C T ′ ,n where T ⊂ T ′ ∈ T n-1,n .
The measure (2.12) is inspired by the penalty terms used in IPDG methods, cf., e.g., Arnold [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF]. Owing to the flux approximation property in Assumption 4.1, the design of J u,NC (u hτ ) is linked to the way fluxes are equilibrated. Consequently, the positive weights

C K,φ,T,F,n (physical unit is L 3 T -2
) have to be chosen accordingly; we refer to (7.7) for an example. Moreover, the summation over the common refinement T n-1,n is linked to continuous-in-time flux reconstruction in the Crank-Nicolson scheme, whereas T n and piecewise constant-in-time flux reconstruction appear naturally in the backward Euler case. Furthermore, a straightforward consequence of the fact that the weak solution u is in X is the following lemma:

Lemma 2.4 (Jumps of the weak solution). For all 1 ≤ n ≤ N , all T ∈ T n-1,n , and all F ∈ F T ,

[[u]] = 0 in L 2 (F × I n ). Owing to Lemma 2.4, [[u -u hτ ]] can be replaced by [[u hτ ]] in (2.12). Thus, in contrast to J u,FR (u hτ ), J u,NC (u hτ ) is easily computable. Finally, note that J u,NC (u hτ ) = 0 if and only if u hτ ∈ X, that is, if and only if u hτ is X-conforming.

Error measure

Our error measure is the sum of (2.7) and (2.12), i.e. J u (u hτ ) := J u,FR (u hτ ) + J u,NC (u hτ ).

(

2.13)

There holds J u (u hτ ) = 0 if and only if u hτ = u. Indeed, if u hτ = u, J u (u hτ ) clearly equals 0. Conversely, if J u (u hτ ) = 0, then J u,FR (u hτ ) = J u,NC (u hτ ) = 0. Thus, u hτ is in X and, hence, since the weak solution is uniquely characterized by the property R(u), ϕ Y ′ ,Y = 0 for all ϕ ∈ Y , see (2.2) and (2.8), we infer u hτ = u.

A posteriori error estimate

This section collects the main results of this paper concerning the error upper bound. The approximation error is measured by (2.13), and the estimators are defined using an equilibrated flux reconstruction.

Equilibrated flux reconstruction

In order to proceed as generally as possible, in particular without the definition of any numerical scheme for approximating (1.1a)-(1.1c), we make the following assumption (recall that H(div, Ω) is spanned by vector fields in [L 2 (Ω)] d with weak divergence in L 2 (Ω), cf. Brezzi and Fortin [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]):

Assumption 3.1 (Space-time equilibrated flux reconstruction). There exists a flux reconstruction t hτ ∈ L 2 (0, t F ; H(div, Ω)), equilibrated in the sense that

(f -∂ t u hτ -∇•t hτ , 1) T ×In = 0 ∀1 ≤ n ≤ N, T ∈ T n-1,n . (3.1)
The function t hτ , which represents the advective-diffusive flux, has the opposite sign with respect to the function σ of (1.2), usually defined with the same sign as the gradient to yield monotonicity-type properties. Specific constructions of the flux t hτ for various spatial discretizations are discussed in [START_REF]A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF]. We present an example in the context of IPDG methods for problem (1.1a)-(1.1c) in Section 7.

Remark 3.2 (Local mass conservation). Equation (3.1) expresses local mass conservation over the spacetime element T × I n . A similar assumption has been made in [START_REF]A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF], see equation (3.4) therein. In [START_REF]A posteriori error estimation based on potential and flux reconstruction for the heat equation[END_REF], however, the local mass conservation had to be satisfied over a given mesh element T for all times t ∈ I n . The present assumption, being more general, allows for more flexibility. In particular, it allows us to use the scaled Poincaré inequality (2.4) in the proof of the error upper bound. This local mass conservation is supposed in (3.1) on the elements of the common coarsening of T n-1 and T n only, which in particular fits the Crank-Nicolson scheme of Section 7. Local mass conservation on the elements of T n appears naturally in the backward Euler case.

Guaranteed a posteriori error estimate

We are now in a position to state our main result concerning the error upper bound. For all 1 ≤ n ≤ N and all T ∈ T n-1,n , we define the residual, flux, and nonconformity estimators respectively

η n R,T := C -1 2 T,n C P f -∂ t u hτ -∇•t hτ T ×In , (3.2)
η n F,T := C -1 2 T,n h -1 T σ(u hτ , ∇u hτ ) + t hτ T ×In , (3.3)
η n NC,T :=    T ′ ∈T n-1,n , T ′ ⊂T F ∈F T ′ C -1 T ′ ,n h -2 T ′ C K,φ,T ′ ,F,n [[u hτ ]] 2 F ×In    1 2 , (3.4)
and set η n FR,T := η n F,T + η n R,T . These estimators are local-in-space and in-time. We define their global space-time versions as η

• := N n=1 T ∈T n-1,n η n •,T
η IC :=    N n=1 T ∈T n-1,n (η n IC,T ) 2    1 2 , η n IC,T := C -1 2 T,n (τ n ) -1 2 u 0 -u hτ (•, 0) T . (3.5)
Theorem 3.3 (Guaranteed a posteriori error estimate). Let u ∈ X be the weak solution given by (2.2) and let u hτ ∈ X h be arbitrary. Let Assumption 3.1 hold. Let η n R,T , η n F,T , η n NC,T , and η IC be defined by (3.2)-(3.5). Then,

J u (u hτ ) ≤ η FR + η NC + η IC . (3.6)
The proof is given in Section 5. We now present several remarks.

Remark 3.4 (Interpretation of the error estimators). The flux estimator η n F,T is related to the violation of the constitutive relation (1.2), the residual estimator η n R,T to the violation of the equilibrium condition (1.1a), while η n NC,T and η n IC,T are related to the violation of the constraints (u ∈ X and u(•, 0) = u 0 ) at the discrete level.

Remark 3.5 (Spatial and temporal errors). The decomposition of our estimators into spatial and temporal error estimators depends on the time-marching scheme. For the flux reconstruction t hτ within the Crank-Nicolson setting of Section 7, η n F,T and η n NC,T can be viewed as spatial error estimators and η n R,T as temporal error estimator. Remark 3.6 (Initial condition estimator). Whenever the weights C T,n are chosen according to (2.6), there holds C T,n ≥ t F (τ n ) -2 , whence we readily infer that η IC ≤ u hτ (•, 0) -u 0 Ω , that is, η IC is upper-bounded by the usual L 2 (Ω)-norm of the approximation error on the initial condition. Moreover, for very small time steps, η IC tends to u hτ (•, 0) -u 0 Ω , while for very large time steps, η IC becomes very small.

Efficiency and robustness

This section deals with the efficiency and robustness of our estimates. For simplicity, we assume that the source term f is a piecewise space-time polynomial on T n-1,n ; otherwise, a classical data oscillation term has to be included in the error lower bound. We likewise assume that both u hτ and t hτ are piecewise space-time polynomials on T n-1,n . We observe that, because of the nonlinear functions K(•) and φ(•), the flux σ(u hτ , ∇u hτ ) is not necessarily a piecewise polynomial, even if u hτ is. Let

σ(u hτ , ∇u hτ )| T := P T σ(u hτ , ∇u hτ )| T , ∀1 ≤ n ≤ N, T ∈ T n-1,n , (4.1)
where P T is a projection-type operator mapping onto piecewise space-time polynomials, see (7.5)-(7.6) below for an example. We define the quadrature error estimator

η n qd,T := C -1 2 T,n h -1 T σ(u hτ , ∇u hτ ) -σ(u hτ , ∇u hτ ) T ×In , (4.2)
together with its global space-time version η qd :=

N n=1 T ∈T n-1,n η n qd,T 2 1 2 .
Henceforth, A B means that there exists a constant C such that A ≤ CB, with C only depending on space dimension, maximal polynomial degree, shape-regularity of the meshes T n for all 1 ≤ n ≤ N , the maximal ratio h T /h T ′ over all 1 ≤ n ≤ N , T ∈ T n-1,n (the common coarsening of T n-1 and T n ) and T ′ ∈ T n-1,n (the common refinement of T n-1 and T n ), T ′ a subelement of T , and the maximal ratios

C T,n /C T ′ ,n
, for all T ∈ T n-1,n and T ′ ∈ T n-1,n sharing a face and all 1 ≤ n ≤ N . We thus in particular suppose that both refinement and coarsening are not too abrupt.

Approximation property of the flux reconstruction

Let, for all 1 ≤ n ≤ N and all T ∈ T n-1,n , η n clas,T be given by

h T f -∂ t u hτ + ∇• σ(u hτ , ∇u hτ ) T ×In +    F ∈F int T h F [[σ(u hτ , ∇u hτ )]]•n F 2 F ×In    1 2 + F ∈FT C K,φ,T,F,n [[u hτ ]] 2 F ×In 1 2 . ( 4.3)
The quantity η n clas,T can be viewed as a classical residual-based a posteriori error estimator. In order to carry the analysis without specifying a numerical scheme, we make the following assumption on the reconstructed flux: Assumption 4.1 (Flux approximation property). There holds

σ(u hτ , ∇u hτ ) + t hτ 2 T ×In T ′ ∈T n-1,n , T ′ ⊂T (η n clas,T ′ ) 2 ∀1 ≤ n ≤ N, T ∈ T n-1,n .
An example on how to achieve Assumption 4.1 in the context of IPDG methods in space and the Crank-Nicolson scheme in time is presented in Section 7.

Robust local-in-space and in-time efficiency

To control quadrature errors, we introduce coefficients 0 < γ n qd,T and require max

T ′ ∈TT η n qd,T ′ ≤ γ n qd,T (η n FR,T + η n NC,T ), ∀1 ≤ n ≤ N, T ∈ T n-1,n , (4.4)
recalling that T T collects the mesh elements of T n-1,n sharing a face with T .

Robust global efficiency

As above, in order to control quadrature errors, we introduce a coefficient 0 < γ qd and require

η qd ≤ γ qd (η FR + η NC ). (4.6)
Then we have the following full equivalence result between our estimators and J u (u hτ ).

Theorem 4.4 (Robust global efficiency). Let (4.6), with γ qd small enough, be satisfied. Let Assumption 4.1 hold. Then,

η FR + η NC J u (u hτ ). (4.7)
Remark 4.5 (Weights C T,n ). We stress that our robustness results are not influenced by the value assigned to the weights C T,n of (2.5a). Indeed, multiplying them by a positive factor λ scales both error measures J u,FR (u hτ ) and J u,NC (u hτ ) defined by (2.7) and (2.12) by the factor λ -1 2 , while the error estimators η n R,T , η n F,T , η n NC,T , and η n IC,T defined by (3.2)-(3.5) are scaled by the same factor. Hence, the ratio of error measure J u (u hτ ) to error indicators is independent of the scaling factor λ.

Remark 4.6 (Steady case). The present approach can be applied to the steady case leading to robustness with respect to advection dominance. Moreover, choosing

C T := h -2 T C φ,T + h -2 T C K,T with C φ,T := h T φ ′ (u hτ ) ∞,T , C K,T := K(u hτ ) ∞,T
, the steady version of the error measure e FR of (2.10) is multiplied by the square-root of the Péclet number with respect to the usual energy norm, while the estimators are divided by this factor with respect to the classical estimators (without cutoff factors), see, e.g., Verfürth [START_REF]Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF] and [START_REF] Ern | Guaranteed and robust discontinuous Galerkin a posteriori error estimates for convection-diffusion-reaction problems[END_REF]. Thus, a similarity to the approach of Sangalli [START_REF] Sangalli | Robust a-posteriori estimator for advection-diffusion-reaction problems[END_REF] can be observed.

Proof of the error upper bound

The goal of this section is to prove Theorem 3.3. The proof is decomposed in two steps. Lemma 5.1 (Bound on J u,FR (u hτ )). There holds J u,FR (u hτ ) ≤ η FR + η IC .

Proof. Let ϕ ∈ Y with ϕ Y = 1 be given. Using (2.8), subtracting the quantity tF 0 (t hτ , ∇ϕ)(t) dt + tF 0 (∇•t hτ , ϕ)(t) dt, which is equal to zero owing to the Green theorem since t hτ ∈ L 2 (0, t F ; H(div, Ω)) by Assumption 3.1, integrating by parts in time since ∂ t u hτ ∈ L 2 (Q) by assumption, and using that

ϕ(•, 0) = - tF 0 ∂ t ϕ(t) dt since ϕ(•, t F ) = 0 by assumption, we infer that R(u hτ ), ϕ Y ′ ,Y is equal to tF 0 {(f -∂ t u hτ -∇•t hτ , ϕ) + (u hτ (•, 0) -u 0 , ∂ t ϕ) -(σ(u hτ , ∇u hτ ) + t hτ , ∇ϕ)}(t) dt.
We now employ the second part of Assumption 3.1, namely (3.1), yielding

R(u hτ ), ϕ Y ′ ,Y = N n=1 T ∈T n-1,n (f -∂ t u hτ -∇•t hτ , ϕ -Π 0 ϕ) T ×In + (u hτ (•, 0) -u 0 , ∂ t ϕ) T ×In -(σ(u hτ , ∇u hτ ) + t hτ , ∇ϕ) T ×In ≤ N n=1 T ∈T n-1,n η n R,T + (η n F,T ) 2 + (η n IC,T ) 2 1 2 ϕ Y,T ×In ,
where we have used the Cauchy-Schwarz inequality, the scaled space-time Poincaré inequality (2.4), and the definition (2.5a). The assertion then follows from (η n F,T ) 2 + (η n IC,T ) 2 1 2 ≤ η n F,T + η n IC,T , the Cauchy-Schwarz inequality, the definition (2.5b) of the • Y -norm, the fact that ϕ Y = 1, and the relation (2.9). Lemma 5.2 (Bound on J u,NC (u hτ )). There holds J u,NC (u hτ ) = η NC .

Proof. Immediate owing to Lemma 2.4 and (3.4).

Proof of the error lower bound

We present here the proofs of Theorems 4.2 and 4.4. An important ingredient is the bubble function technique (see [START_REF]A posteriori error estimates for finite element discretizations of the heat equation[END_REF][START_REF]Robust a posteriori error estimates for nonstationary convection-diffusion equations[END_REF]) that we extend here to space-time bubbles.

Proof of Theorem 4.2

Proof. Let a time step 1 ≤ n ≤ N and a mesh element T ∈ T n-1,n be fixed. The first two steps of the proof are devoted to bounding the quantity η n clas,T defined by (4.3), while staying on the common refinement T n-1,n .

Step 1, bound on C

-1 2 T,n f -∂ t u hτ + ∇•σ(u hτ , ∇u hτ ) T ×In . Let us prove that C -1 2 T,n f -∂ t u hτ + ∇•σ(u hτ , ∇u hτ ) T ×In e n FR,T + η n qd,T , (6.1)
with e n FR,T defined by (2.11). Set v T,n := (f -∂ t u hτ + ∇•σ(u hτ , ∇u hτ ))| T ×In . Let ψ T,n be the space-time bubble function on T × I n given by the product of the barycentric coordinates on T and of the barycentric coordinates on I n . Note that both ψ T,n and v T,n are (space-time) polynomials since we are on a refinement of both T n-1 and T n . Then, by norm equivalence in finite-dimensional spaces, there holds

(v T,n , v T,n ) T ×In (v T,n , ψ T,n v T,n ) T ×In . (6.2)
Using an inverse inequality separately in space and in time, we obtain

h T ∇(ψ T,n v T,n ) T ×In ψ T,n v T,n T ×In , (6.3a)
τ n ∂ t (ψ T,n v T,n ) T ×In ψ T,n v T,n T ×In . (6.3b)
The norm • Y,T ×In defined by (2.5a) was precisely designed in order to use these inequalities at the present stage. Using (6.3), (2.5a), and ψ T,n ∞,T,n ≤ 1, we infer

C -1 T,n ψ T,n v T,n 2
Y,T ×In = h 2 T ∇(ψ T,n v T,n ) 2 T ×In + (τ n ) 2 ∂ t (ψ T,n v T,n ) 2 T ×In ψ T,n v T,n 2
T ×In ≤ v T,n 2
T ×In .

(6.4) Thus, using the definition of v T,n , (6.2), (2.2), and the Green theorem, yields, with the notation σ d := σ(u, ∇u)σ(u hτ , ∇u hτ ),

C -1 T,n v T,n 2 T ×In C -1 T,n (f -∂ t u hτ + ∇•σ(u hτ , ∇u hτ ), ψ T,n v T,n ) T ×In = C -1 T,n (u hτ -u, ∂ t (ψ T,n v T,n )) T ×In +(σ d , ∇(ψ T,n v T,n )) T ×In . (6.5)
Thus, owing to (6.4),

C -1 2 T,n v T,n T ×In (u hτ -u, ∂ t (ψ T,n v T,n )) T ×In + (σ d , ∇(ψ T,n v T,n )) T ×In ψ T,n v T,n Y,T ×In ,
whence (6.1) follows from the Cauchy-Schwarz inequality, the definitions (2.11) of e n FR,T and (4.2) of η n qd,T , and the triangle inequality.

Step 2, bound on C

-1 2 T,n h -1 2 T [[σ(u hτ , ∇u hτ )]]•n F F ×In for all F ∈ F int T .
Let F ∈ F int T be fixed. Recalling that T F denotes the simplices in T n-1,n sharing the face F , let us prove that F v F,n F ×In . Using these inequalities and proceeding as in step 1, we arrive at

C -1 2 T,n h -1 2 T [[σ(u hτ , ∇u hτ )]]•n F F ×In T ′ ∈TF (e n FR,T ′ + η n qd,T ′ ). ( 6
C -1 2 T,n h -1 2 T v F,n F ×In C -1 2 T,n f -∂ t u hτ + ∇•σ(u hτ , ∇u hτ ) TF ×In + (u -u hτ , ∂ t (ψ F,n v F,n )) TF ×In -(σ d , ∇(ψ F,n v F,n )) TF ×In ψ F,n v F,n Y,TF ,n .
Finally, the first term on the right-hand side is bounded using (6.1) for each T ′ ∈ T F , while proceeding as in step 1 for the second term yields (6.6).

Step 3, conclusion. Let now 1 ≤ n ≤ N and an element T from the common coarsening T n-1,n be fixed. Combining (6.1) and (6.6) yields

T ′ ∈T n-1,n , T ′ ⊂T (C -1 2 T ′ ,n h -1 T ′ η n clas,T ′ ) 2 T ′ ∈TT (e n FR,T ′ + η n qd,T ′ ) 2 + (η n NC,T ) 2 .
Using the triangle inequality and Assumption 4.1, we infer (η n

F,T ) 2 (η n qd,T ) 2 +C -1 T,n h -2 T T ′ ∈T n-1,n , T ′ ⊂T (η n clas,T ′ ) 2 .
Similarly, using in addition the inverse space inequality leads to (

η n R,T ) 2 C -1 T,n h -2 T T ′ ∈T n-1,n , T ′ ⊂T (η n clas,T ′ ) 2 .
Combining these inequalities we obtain

η n FR,T T ′ ∈TT (e n FR,T ′ + η n qd,T ′ ) + η n NC,T .
Finally, invoking the requirement (4.4) on quadrature errors to discard the quantities η n qd,T ′ from the righthand side completes the proof.

Proof of Theorem 4.4

Proof. We use the notations v T,n and v F,n from the previous proof. It is sufficient to show that

   N n=1 T ∈T n-1,n C -1 T,n v T,n 2 T ×In    1 2
J u (u hτ ) + η qd , (6.7a)

   N n=1 T ∈T n-1,n F ∈F int T C -1 T,n h -1 T v F,n 2 F ×In    1 2 J u (u hτ ) + η qd , (6.7b)
since choosing γ qd in (4.6) small enough then yields (4.7). We only show (6.7a); (6.7b) follows by similar arguments. Let 1 ≤ n ≤ N and let T ∈ T n-1,n . It follows from (6.5) and the Green theorem that

C -1 T,n v T,n 2 T ×In C -1 T,n {(f -∂ t u hτ , ψ T,n v T,n ) T ×In -(σ(u hτ , ∇u hτ ), ∇(ψ T,n v T,n )) T ×In +(σ(u hτ , ∇u hτ )-σ(u hτ , ∇u hτ ), ∇(ψ T,n v T,n )) T ×In }. Set λ| T ×In := C -1
T,n ψ T,n v T,n and observe that λ ∈ Y . Recall the notation R(u hτ ) from (2.8), J u,FR (u hτ ) from (2.7), and η n qd,T from (4.2). Summing the above inequality over all 1 ≤ n ≤ N and all T ∈ T n-1,n and using the Green theorem, the Cauchy-Schwarz inequality, (6.3a), and ψ T,n ∞,T,n ≤ 1 yields

N n=1 T ∈T n-1,n C -1 T,n v T,n 2 T ×In R(u hτ ), λ Y ′ ,Y + N n=1 T ∈T n-1,n η n qd,T C -1 2 T,n v T,n T ×In ≤ J u,FR (u hτ ) λ Y + η qd    N n=1 T ∈T n-1,n C -1 T,n v T,n 2 T ×In    1 2 . (6.8)
The definition (2.5b) of the • Y -norm, that of λ, the fact that we suppose that the ratio h T /h T ′ , T ∈ T n-1,n , T ′ ∈ T n-1,n , T ′ ⊂ T , is bounded, and (6.4) lead to

λ 2 Y = N n=1 T ∈T n-1,n C -2 T,n ψ T,n v T,n 2 Y,T ×In N n=1 T ∈T n-1,n C -1 T,n v T,n 2
T ×In . (6.9)

Combining (6.8) and (6.9) proves (6.7a).

Application to interior-penalty discontinuous Galerkin methods

We apply the framework of Sections 3 and 4 to IPDG methods as an example of nonconforming space discretization, with Crank-Nicolson time stepping. This consists in specifying the flux reconstruction t hτ and in verifying Assumptions 3.1 and 4.1.

The IPDG method

For all 0 ≤ n ≤ N , let V n h := P p (T n ), p ≥ 1, be spanned by piecewise polynomials on the mesh T n with total degree ≤ p. For each F ∈ F n , along with the jump operator [[•]], we consider the (arithmetic) average operator { {•} } (conventionally yielding the actual value at boundary faces).

The space-time approximation u hτ is continuous and piecewise affine in time, and is defined by its values u n h at t n for all 0 ≤ n ≤ N . We take u 0 h as the L 2 -orthogonal projection of u 0 onto V 0 h . Then, for all 1 ≤ n ≤ N , we look for

u n h ∈ V n h such that (∂ t u n h , v h ) + 1 2 n m=n-1 (σ(u m h , ∇u m h ), ∇v h ) + F ∈F m α m K,F h -1 F ([[u m h ]], [[v h ]]) F + F ∈F m (H F (u m h ), [[v h ]]) F - F ∈F m ({ {K(u m h )∇u m h } }•n F , [[v h ]]) F -θ F ∈F m ({ {K(u m h )∇v h } }•n F , [[u m h ]]) F -(f m , v h ) = 0 ∀v h ∈ V n h , (7.1)
where α n K,F is the penalty coefficient, θ ∈ {-1, 0, 1}, and

∂ t u n h := (τ n ) -1 (u n h -u n-1 h
); we suppose for simplicity that f is continuous and piecewise affine in time and denote f n := f (t n ). The advection term in (7.1) has been discretized using a numerical flux H F (u n h ) satisfying the following reasonable assumption:

Assumption 7.1 (Numerical flux for advection). For all 1 ≤ n ≤ N , all mesh faces F ∈ F n , and all v h ∈ V n h , there holds

H F (v h ) -{ {φ(v h )} }•n F F φ ′ (v h ) ∞,TF [[v h ]] F .
An example is the numerical flux of Lax-Friedrichs type, which consists of the centered flux { {φ(•)} }•n F supplemented by a stabilization term penalizing the interface jumps. In our numerical experiments, we employ a numerical upwinding flux, see Section 8.1.

Flux reconstruction

Let l ≥ 0. We construct t hτ continuous and piecewise affine in time, with

t n h := t hτ (t n ), 0 ≤ n ≤ N , in the Raviart-Thomas-Nédélec space RTN l (T n ) := v h ∈ H(div, Ω); v h | T ∈ RTN l (T ) for all T ∈ T n , where RTN l (T ) := [P l (T )] d + xP l (T ). In particular, v h ∈ RTN l (T n ) is such that, for all T ∈ T n , ∇•v h ∈ P l (T ),
and for all F ∈ F n , v h •n F ∈ P l (F ), with continuity across interfaces, cf. Brezzi and Fortin [START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF]. Following Kim [START_REF] Kim | A posteriori error estimators for locally conservative methods of nonlinear elliptic problems[END_REF] and [START_REF] Ern | An accurate H(div) flux reconstruction for discontinuous Galerkin approximations of elliptic problems[END_REF], we set: Definition 7.2 (Reconstructed flux t hτ ). Let l ≥ 0. Set w F := 1 2 for interfaces and w F := 1 for boundary faces. For all 0 ≤ n ≤ N , we specify t n h ∈ RTN l (T n ) by setting, for all T ∈ T n , all F ∈ F T , all q h ∈ P l (F ),

and all r h ∈ [P l-1 (T )] d , (t n h •n F , q h ) F = (-{ {K(u n h )∇u n h } }•n F + α n K,F h -1 F [[u n h ]] + H F (u n h ), q h ) F , (7.2)
(t n h , r h ) T = (-K(u n h )∇u n h + φ(u n h ), r h ) T + θ F ∈FT w F (K(u n h )r h •n F , [[u n h ]]) F . (7.3)
Remark 7.3 (Alternative construction). Instead of prescribing directly the degrees of freedom for t n h , it is also possible to reconstruct the flux by solving local Neumann problems by mixed finite elements, see [START_REF] Ern | Flux reconstruction and a posteriori error estimation for discontinuous Galerkin methods on general nonmatching grids[END_REF]. This approach can achieve a tighter relationship between the error and the estimator but is slightly more expensive.

Verification of Assumption

3.1 Lemma 7.4 (Local conservation). Let u n h , 1 ≤ n ≤ N , solve (7.1
) and let t hτ be defined by (7.2)-( 7.3) with l ≥ 0. Then, t hτ satisfies t hτ ∈ L 2 (0, t F ; H(div, Ω)) and, for all 1 ≤ n ≤ N and all T ∈ T n-1,n ,

(f -∂ t u hτ -∇•t hτ , v h ) T ×In = 0 ∀v h ∈ P min(p,l) (T ), (7.4)
so that, in particular, the local space-time conservation property (3.1) holds.

Proof. By construction, we have, for 0

≤ n ≤ N , t n h ∈ H(div, Ω), so that t hτ ∈ L 2 (0, t F ; H(div, Ω)). Let 1 ≤ n ≤ N , T ∈ T n-1
,n , and v h ∈ P min(p,l) (T ) be given. Using that u hτ , f , and t hτ are affine in time on I n , that T is from the common coarsening of the meshes T n-1 and T n , and the Green theorem, we obtain

(f -∂ t u hτ -∇•t hτ , v h ) T ×In = -τ n (∂ t u hτ , v h ) T + τ n 2 n m=n-1 (f m -∇•t m h , v h ) T =-τ n ∂ t u hτ - 1 2 n m=n-1 f m , v h T + τ n 2 n m=n-1 T ′ ∈T m , T ′ ⊂T {(t m h , ∇v h ) T ′ -(t m h •n T ′ , v h ) ∂T ′ }.
Let m = n or m = n -1. The facts that v h | F ∈ P l (F ) for all F ∈ F T ′ , (7.2), ∇v h ∈ [P l-1 (T ′ )] d for any T ′ , and (7.3) yield

(t m h •n T ′ , v h ) ∂T ′ = - F ∈F T ′ n T ′ •n F {(-{ {K(u m h )∇u m h } }•n F +α m K,F h -1 F [[u m h ]]+H F (u m h ), v h ) F }, (t m h , ∇v h ) T ′ = (-K(u m h )∇u m h + φ(u m h ), ∇v h ) T ′ + θ F ∈F T ′ w F (K(u m h )∇v h •n F , [[u m h ]]) F .
Extending v h by 0 outside T so that a function in V n h is obtained and using the above identities and the definition (7.1) of the scheme yields (7.4).

Verification of Assumption 4.1

To verify Assumption 4.1, we first specify the space-time piecewise polynomial σ(u n h , ∇u n h ). Let 0 ≤ n ≤ N and let T ∈ T n . We define σ(u n h , ∇u n h )| T ∈ RTN l (T ) such that, for all F ∈ F T , all q h ∈ P l (F ), and all

r h ∈ [P l-1 (T )] d , (σ(u n h , ∇u n h )| T •n F , q h ) F = (σ(u n h , ∇u n h )| T •n F , q h ) F , (7.5)
(σ(u n h , ∇u n h ), r h ) T = (σ(u n h , ∇u n h ), r h ) T . (7.6)
Here, l is the polynomial degree used for reconstructing the flux t n h in Section 7.2. We observe that, locally in each mesh element

T ∈ T n , σ(u n h , ∇u n h )| T belongs to RTN l (T ) (as t n h does), but, globally, σ(u n h , ∇u n h ) ∈ RTN l (T n ) because the normal component of σ(u n h , ∇u n h
) is in general discontinuous across interfaces. Finally, the space-time function σ(u hτ , ∇u hτ ) is taken to be continuous and piecewise affine in time, matching the values σ(u n h , ∇u n h ) at t n for all 0 ≤ n ≤ N .

Lemma 7.5 (Flux approximation). Let t hτ be defined by (7.2)-( 7.3) with l ≥ 0 and let σ(u hτ , ∇u hτ ) be defined by (7.5)- (7.6). For 1 ≤ n ≤ N and T ∈ T n-1,n , define the constants C K,φ,T,F,n of (2.12) by

C K,φ,T,F,n := (α n K,F ) 2 h -1 F + K(u n h ) 2 ∞,T h -1 F + φ ′ (u n h ) 2 ∞,TF h F . (7.7)
Then, Assumption 4.1 holds.

Proof. Let 0 ≤ n ≤ N and set

v h := σ(u n h , ∇u n h ) + t n h . Recall that σ(u n h , ∇u n h ) = K(u n h )∇u n h -φ(u n h ).
Let T ∈ T n and F ∈ F T . Using (7.2), (7.5) and (7.3), (7.6) yield, respectively,

(v h | T •n F , q h ) F = ((1 -w F )[[σ(u n h , ∇u n h )]]•n F + α n K,F h -1 F [[u n h ]] + H F (u n h ) -{ {φ(u n h )} }•n F , q h ) F , (v h , r h ) T = θ F ∈FT w F (K(u n h )r h •n F , [[u n h ]]) F ,
for all q h ∈ P l (F ) and all r h ∈ [P l-1 (T )] d . Thus, using Assumption 7.1 on the numerical flux for advection and employing an inverse inequality in space, we infer

v h •n F F (1 -w F ) [[σ(u n h , ∇u n h )]]•n F F + (α n K,F h -1 F + φ ′ (u n h ) ∞,TF ) [[u n h ]] F , (v h , r h ) T K(u n h ) ∞,T r h T F ∈FT h -1 2 F [[u n h ]] F .
Using these two estimates in the classical bound v h

2 T F ∈FT h F v h •n F 2 F + sup r h ∈[P l-1 (T )] d (v h ,r h )T r h T 2 ,
valid for all v h ∈ RTN l (T ), and owing to (7.7),

σ(u n h , ∇u n h ) + t n h 2 T F ∈F int T h F [[σ(u n h , ∇u n h )]]•n F 2 F + F ∈FT C K,φ,T,F,n [[u n h ]] 2 F . (7.8)
Let now 1 ≤ n ≤ N and T ∈ T n-

σ(u m h , ∇u m h ) + t m h 2 T = τ n n m=n-1 T ′ ∈T m , T ′ ⊂T σ(u m h , ∇u m h ) + t m h 2 T ′ T ′ ∈T n-1,n , T ′ ⊂T (η n clas,T ′ ) 2 ,
employing (7.8) on T ′ ∈ T m ; actually, only the jump terms of η n clas,T ′ appear.

Remark 7.6 (Choice of the reconstruction degree l). A typical choice for the polynomial degree in the flux reconstruction t hτ is l ∈ {p -1, p}. Choosing larger values for l ensures theoretically that the quadrature errors are small enough. In our numerical experiments, the choice l = p was sufficient to deliver quadrature errors one to two orders of magnitude smaller than the leading estimators.

Numerical experiments

In this section, we present numerical experiments illustrating our a posteriori error estimates.

Setting

We consider the problem (1.1a)-(1.1c) where, in some cases, we replace the homogeneous Dirichlet boundary condition (1.1b) by an inhomogeneous one. The error for nonpolynomial Dirichlet boundary data is neglected, since it is generally of higher order; for a numerical study, see, e.g., Luce and Wohlmuth [START_REF] Luce | A local a posteriori error estimator based on equilibrated fluxes[END_REF].

We employ the IPDG method (7.1) with θ = 0 and the upwind numerical flux such that, for all 1 ≤ n ≤ N and all

F ∈ F n , H F (u n h )| F := φ(u n h )| L F •n F if { {φ ′ (u n h )} } F •n F > 0 and H F (u n h )| F := φ(u n h )| R F •n F otherwise, where v n h | L F and v n h | R
F denote traces of a function v n h on F from the direction and the opposite direction of n F , respectively. The penalty parameter in (7.1) is chosen as α n K,F := 10p 2 K(u n h ) ∞,TF ; see Houston et al. [START_REF] Houston | Discontinuous Galerkin finite element approximation of quasilinear elliptic boundary value problems I: The scalar case[END_REF], [START_REF] Dolejší | Analysis and application of the IIPG method to quasilinear nonstationary convectiondiffusion problems[END_REF], and [START_REF] Ern | A discontinuous Galerkin method with weighted averages for advection-diffusion equations with locally small and anisotropic diffusivity[END_REF] for problems with internal layers caused by locally small diffusion. We consider polynomial degrees p ∈ {1, 2, 3}. For the linearization of (7.1), we use the Newton-like method of [START_REF] Dolejší | Efficient solution strategy for the semi-implicit discontinuous Galerkin discretization of the Navier-Stokes equations[END_REF] where the (approximate) construction of the Jacobian matrix is avoided via the idea of easy-to-evaluate flux matrix. The volume integrals are evaluated with the Dunavant quadrature rule of order (3p + 2), and the face integrals with the Gauss quadrature rule with (2p + 2) nodes.

We consider square domains Ω with uniform discretizations (consisting of right-angled triangles obtained by diagonal cuttings of squares). The time step and the meshes are kept fixed in time. We consider a family of uniformly refined space-time meshes characterized by the squares edge lengths and time steps (h m,p , τ m,p ), m ∈ {1, 2, 3}, p ∈ {1, 2, 3}. We choose (h 0 , τ 0 ) and then set h m,p := h 0 2 1-m , m, p ∈ {1, 2, 3}, and τ m,1 := τ 0 2 (1-m)/2 , τ m,2 := 0.5τ 0 2 (1-m) , τ m,3 := 0.2τ 0 2 3(1-m)/2 , m ∈ {1, 2, 3}. We evaluate experimental orders of convergence in the form EOC := log(Em,p/Em-1,p) log(hm,p/hm-1,p) , m ∈ {2, 3}, where E m,p is either an error or an error estimator on the space-time discretization (h m,p , τ m,p ).

The weights C T,n are chosen according to (2.6) with C φ,T,n := h Ω φ ′ (u hτ ) ∞,Q (here, h Ω denotes the diameter of Ω) and C K,T,n := (h Ω /h T ) K(u hτ ) ∞,Q . The first two addends on the right-hand side of (2.6) are of similar size if the Courant numbers τ n φ ′ (u hτ ) ∞,Q /h T are of order unity, as well as the ratios t F φ ′ (u hτ ) ∞,Q /h Ω (meaning that the final time allows particles to be advected across a relevant part of the domain). Moreover, the ratio of the second to the third addend is of the order of the Péclet numbers

h T φ ′ (u hτ ) ∞,Q / K(u hτ ) ∞,Q .
To evaluate the error measure J u,FR (u hτ ), we solve (approximately) the dual problem introduced in Section 2.3.1. This problem, which is posed on the space-time domain Q, is solved using the finite difference package Fishpack [START_REF] Adams | FISHPACK -efficient FOR-TRAN subprograms for the solution of separable elliptic partial differential equations[END_REF] on a structured space-time mesh obtained by typically three successive dyadic partitions of the space-time grid on which the approximate solution u hτ has been computed. Recalling that J u (u hτ ) = J u,FR (u hτ ) + J u,NC (u hτ ) where J u,NC (u hτ ) coincides with the nonconformity estimator η NC , we evaluate the effectivity index i e := η/J u (u hτ ), where η := η FR +η NC +η IC to illustrate Theorems 3.3 and 4.4. Moreover, to illustrate Theorem 4.2, we evaluate the effectivity index i e,FR := η/(e FR + J u,NC (u hτ )) where e FR is the locally computable upper bound on J u,FR (u hτ ), see (2.10)-(2.11).

Linear diffusion

We consider first a model advection-diffusion problem with Ω = (-1, 1) × (-1, 1), t F = 1, and linear diffusive but nonlinear advective flux. In particular, K(u) = εI (I being the identity tensor), ε > 0, and φ(u) = (u 2 /2, u 2 /2) T in (1.2). The initial and boundary conditions are such that the exact solution is u(x, y, t) = (1 + exp (ϑ(x, y, t))) -1 (8.1) with ϑ(x, y, t) := x+y+1-t 2ε , see [START_REF] Dolejší | An optimal L ∞ (L 2 )-error estimate of the discontinuous Galerkin method for a nonlinear nonstationary convection-diffusion problem[END_REF]. This problem exhibits an inner layer moving in the diagonal direction. The steepness of the layer increases with decreasing ε.

We first employ the initial space-time step (h 0 , τ 0 ) = (1/6, 0.05) and consider the value ε = 10 -2 . Results are reported in Table 1. The flux, residual, nonconforming, and initial condition estimators take relatively close values, while the quadrature estimator is lower. We observe that, for all values of m and p, the effectivity index i e takes stable values in a relatively narrow interval, independently of m and p. The same observation holds for the effectivity index i e,FR , with, however, values ≤ 1, since e FR is only an upper bound on J u,FR (u hτ ).

To assess robustness with respect to advection dominance and relative size of space and time steps, we consider some different choices for ε and (h 0 , τ 0 ). In Table 2, we compare the effectivity indices i e and i e,FR with those of Table 1 (repeated for convenience); for m = p = 3 and (h 0 , τ 0 ) = (1/6, 0.0125), the value of i e is not reported because of the cost of solving accurately the dual problem. We observe that all the considered changes in ε and (h 0 , τ 0 ) only mildly modify the effectivity indices.

Degenerate diffusion

We first consider a nonlinear degenerate advection-diffusion problem from Kačur [START_REF] Kačur | Solution of degenerate convection-diffusion problems by the method of characteristics[END_REF] with Ω = (0, 1)×(0, 1), t F = 1, K(u) = 2εuI, ε = 10 -2 , and φ(u) = 0.5(u 2 , 0) T . The initial and boundary conditions are such that the exact solution is

u(x, y, t) = 1 -exp v(x-vt-x0) 2ε for x ≤ vt + x 0 , 0 for x > vt + x 0 , (8.2)
where x 0 = 1/4 is the initial position of the front. If u = 0 (for x ≥ vt + x 0 ), the diffusive term degenerates. The initial space-time step is (h 0 , τ 0 ) = (1/8, 0.05). Results are summarized in Table 3. We observe again that the quadrature estimator is smaller than the other estimators and that the effectivity indices i e and i e,FR take relatively stable values, independently of m and p.

Our last example is the porous medium equation with Ω = (-6, 6) × (-6, 6), t F = 1, K(u) = κ|u| κ-1 I, κ > 1, and φ(u) = 0. The initial and boundary conditions are chosen so as to yield the so-called Barenblatt solution u(x, y, t) = 1 t + 1 1 -κ -1 4κ 2

x 2 + y 2 (t + 1) 1/κ κ κ-1

+ 1 κ , (8.3)
where [a] + = max(a, 0), a ∈ R. We consider the initial space-time steps (h 0 , τ 0 ) = (0.5, 0.02) and use here τ m,p := τ 0 2 (1-m)p/2 . We first set κ = 2. Results are summarized in Table 4. We can draw the same conclusions as above regarding the effectivity indices; we also observe that in this setting, the initial condition estimator is relatively important. To assess robustness with respect to the size of the nonlinearity, the last column of Table 4 reports the effectivity indices i e and i e,FR obtained for κ = 4; we observe that they take close values to those obtained for κ = 2. Finally, Figure 1 shows (for κ = 4) the distribution of the local error e n T := e n FR,T + η n NC,T and the local error estimator η n T := η n FR,T + η n NC,T at t = t F for m = 2 and p = 2. We observe a close agreement in the distributions of error and estimator, thereby indicating that the latter can be used to drive an adaptive mesh procedure. A similar agreement (not shown) was observed for κ = 2 and the other examples.

2 12

 2 for • ∈ {F, R, FR, NC}. Finally, we define the initial condition estimator as

Theorem 4 . 2 ( 2 +. 5 )

 4225 Robust local-in-space and in-time efficiency). Let a time step 1 ≤ n ≤ N and a mesh element T ∈ T n-1,n be fixed. Let (4.4), with γ n qd,T small enough, hold. Let Assumption 4.1 hold. Recall the definition (2.11) of e n FR,T and define J u,NC,T (u hτ ) as η n NC,T in (3.4), with [[u hτ ]] replaced by [[u -u hτ ]]. J u,NC,T (u hτ ). (4Remark 4.3 (Comment on Theorem 4.2). Estimate (4.5) says that our estimators represent a robust local lower bound for the error measures e n FR,T , whose Hilbertian sum provides an upper bound on the error measure J u,FR (u hτ ), see (2.10), augmented by the jump seminorm J u,NC (u hτ ).

Figure 1 :

 1 Figure 1: Example (8.3), κ = 4: distribution of the element error estimate η n T (left) and of the local error e n T (right) at t = t F for p = 2 and m = 2

  Usingthe same notation for the extension of the function v F,n onto T F by constant values, we also infer the estimate v F,n TF ×In h

	1
	2

.6) Set v F,n := [[σ(u hτ , ∇u hτ )]]•n F | F ×In . Let ψ F,

n be the space-time bubble function on F × I n given by the product of the barycentric coordinates with vertices in F and of the barycentric coordinates on I n . Then, by norm equivalence in finite-dimensional spaces, there holds (v F,n , v F,n ) F ×In (v F,n , ψ F,n v F,n ) F ×In .

  1,n . Using that both σ(u hτ , ∇u hτ ) and u hτ are piecewise affine in time, we have, cf. [17, Lemma 6.1], σ(u hτ , ∇u hτ ) + t hτ

		n
	2 T ×In	τ n
		m=n-1

Table 1 :

 1 Example (8.1), ε = 10 -2 , and (h 0 , τ 0 ) = (1/6, 0.05): errors, estimators, and effectivity indices

	m	p	Ju,FR(u hτ )		ηF	ηR		ηNC	ηIC			η qd	η	ie	ie,FR
	1	1	1.50E-02	1.11E-02	2.28E-02	4.11E-02	2.94E-02	3.82E-03	1.04E-01	1.85	1.15
	2	1	1.17E-02	8.30E-03	1.52E-02	2.29E-02	1.31E-02	1.92E-03	5.94E-02	1.71	1.35
			( 0.36)	( 0.43)	( 0.59)	( 0.84)	( 1.16)		( 0.99)	( 0.81)
	3	1	1.02E-02	5.16E-03	7.78E-03	1.16E-02	2.69E-03	7.49E-04	2.72E-02	1.25	1.36
			( 0.20)	( 0.69)	( 0.96)	( 0.98)	( 2.29)		( 1.36)	( 1.13)
	1	2	4.97E-03	3.78E-03	8.23E-03	1.23E-02	1.32E-02	9.38E-04	3.72E-02	2.15	1.01
	2	2	1.74E-03	1.36E-03	2.52E-03	4.02E-03	1.76E-03	2.34E-04	9.54E-03	1.65	0.94
			( 1.52)	( 1.47)	( 1.71)	( 1.61)	( 2.90)		( 2.00)	( 1.96)
	3	2	4.63E-04	4.00E-04	7.36E-04	1.26E-03	3.01E-04	3.97E-05	2.63E-03	1.53	1.08
			( 1.91)	( 1.77)	( 1.77)	( 1.67)	( 2.55)		( 2.56)	( 1.86)
	1	3	1.78E-03	9.11E-04	1.69E-03	3.41E-03	3.01E-03	2.20E-04	8.88E-03	1.71	0.59
	2	3	3.47E-04	1.57E-04	3.26E-04	6.06E-04	6.20E-04	2.50E-05	1.67E-03	1.75	0.73
			( 2.35)	( 2.54)	( 2.38)	( 2.49)	( 2.28)		( 3.14)	( 2.41)
	3	3	1.33E-05	1.80E-05	3.81E-05	6.97E-05	8.88E-05	1.64E-06	2.10E-04	2.54	0.97
			( 4.71)	( 3.12)	( 3.10)	( 3.12)	( 2.80)		( 3.93)	( 2.99)
			ε		10 -2	10 -2	10 -2		10 -4
			(h0, τ0)	(1/6, 0.05)	(1/6, 0.2)	(1/6, 0.0125)	(1/6, 0.05)
			m	p	ie	ie,FR	ie	ie,FR	ie	ie,FR	ie	ie,FR
			1	1	1.85	1.15	2.21	1.28	3.00	0.81	1.45	0.71
			2	1	1.71	1.35	2.38	1.12	2.45	1.03	1.68	1.06
			3	1	1.25	1.36	2.15	0.90	1.33	1.03	1.82	1.34
			1	2	2.15	1.01	3.13	1.71	3.69	0.67	1.38	0.62
			2	2	1.65	0.94	2.74	1.58	2.16	0.49	1.41	0.62
			3	2	1.53	1.08	2.38	1.52	1.83	0.58	1.54	0.69
			1	3	1.71	0.59	2.74	1.47	3.00	0.34	1.26	0.31
			2	3	1.75	0.73	2.63	1.67	3.15	0.46	1.13	0.21
			3	3	2.54	0.97	2.77	1.73	-	0.69	1.03	0.15

Table 2 :

 2 Example (8.1), effectivity indices for different choices of ε and (h 0 , τ 0 )

Table 3 :

 3 Example (8.2): errors, estimators, and effectivity indices

	κ = 2

Table 4 :

 4 Example (8.3), κ = 2: errors, estimators, and effectivity indices (also for κ = 4)

* This work was supported by the Groupement MoMaS (PACEN/CNRS, ANDRA, BRGM, CEA, EdF, IRSN) and by the ERT project "Enhanced oil recovery and geological sequestration of CO 2 : mesh adaptivity, a posteriori error control, and other advanced techniques" (LJLL/IFPEN). The research of V. Dolejší is supported by the Grant No. 201/08/0012 of the Czech Science Foundation.


Jerome Revaud
email: jerome.revaud@inria.fr

Matthijs Douze
email: matthis.douze@inria.fr

Cordelia Schmid
email: cordelia.schmid@inria.fr

Correlation-Based Burstiness for Logo Retrieval

Keywords: H.2.10 [Image processing and computer vision]: Scene analysis Image retrieval, burstiness, correlation

Detecting logos in photos is challenging. A reason is that logos locally resemble patterns frequently seen in random images. We propose to learn a statistical model for the distribution of incorrect detections output by an image matching algorithm. It results in a novel scoring criterion in which the weight of correlated keypoint matches is reduced, penalizing irrelevant logo detections. In experiments on two very different logo retrieval benchmarks, our approach largely improves over the standard matching criterion as well as other state-of-the-art approaches.

INTRODUCTION

Logo retrieval in real-world images is crucial to quantitatively measure the exposure of brands. In this context, a given logo must be detected and localized in a large collection of images. The preprocessing of the query logo may be costly, as long as the search in the images is fast. This is in contrast with content-based image search engines where the end-to-end search time must be minimized.

Recently, a number of works have tackled logo detection using Bag-of-Words (BoW) techniques [START_REF] Bagdanov | Trademark matching and retrieval in sports video databases[END_REF][START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF][START_REF] Kleban | Spatial pyramid mining for logo detection in natural scenes[END_REF][START_REF] Meng | Interactive visual object search through mutual information maximization[END_REF][START_REF] Romberg | Scalable logo recognition in real-world images[END_REF]. This makes sense since logo images typically contain sharp and contrasted regions that are well handled by keypoint-based representations. Also BoW, in combination with an inverted file, can quickly process large image collections.

The basic BoW model neglects the probabilistic dependence between local features [START_REF] Lewis | Naive (bayes) at forty: The independence assumption in information retrieval[END_REF]. Yet, our observations show that the independence model does not hold: some pairs, triples or n-uplets of keypoints often appear together in incorrect detections. As illustrated in Figure 1, this is due to the presence of local substructures in logo images that are frequently found in random images. For instance, the Adi-Figure 1: Locally, many logos are similar to patterns frequently found in random photos. The Adidas logo for instance, can partially match with building windows, leading to many detection ambiguities.

das logo is locally similar to printed text patterns such as "ll", "li" (with appropriate fonts) or building windows, see Figure 1. As a result, such patterns are often mistaken for the true logo.

In this paper, we show how to down-weight the score of those noisy detections by learning a dedicated burstiness model for the input logo. Section 2 describes how detection hypotheses are obtained by a state-of-the-art approach, then in Section 3 we describe a weighting scheme based on our spatial burstiness approach. Finally, we present experimental results on two logo datasets in Section 4.

DETECTING LOGO HYPOTHESES

Given a logo image L, we are interested in detecting all its occurrences in N test images. Some logos might come in different versions, e.g. the Apple logos shown in Figure 4. In this case, we apply the same detection and scoring procedure for each version, and retain the best score across the different versions.

A logo image is represented as a set of K keypoints {ki}i=1..K . Because the detection and scoring procedures are identical for each test image, without loss of generality, we consider a single test image with keypoints {k ′ i } i=1..K ′ .

Keypoint matching. In order to find other instances of the logo, we first compute a set of matches M = {(ki, k ′ j )} between the logo and the image keypoints, based on the quan- tization of their descriptors into visual words (the codebook size is set to 20,000). To further improve the matching, we compute an additional binary signature for each keypoint and exclude matches for which the Hamming distance between signatures is above 22 [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF]. This has been shown to significantly improve over basic quantization at a minor cost.

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
Spatial verification. Given a set of matches between a logo and an image, we verify the spatial consistency of the matches. As we use scale and rotation invariant keypoints, there exists a single similarity transformation that maps a keypoint to another [START_REF] Philbin | Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF]. We exhaustively examine the similarity transforms corresponding to each match (k, k ′ ) in M, as in [START_REF] David | Distinctive image features from scale-invariant keypoints[END_REF][START_REF] Philbin | Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF]. For each of those transforms, we count the number of inliers with a two-way transfer error [START_REF] Hartley | Multiple View Geometry in Computer Vision[END_REF]. Because logos are small and usually printed on flat surfaces, a similarity transform is an adequate approximation to the full homographic model and much faster to estimate. We set the threshold on the position error to 5 pixels plus a term that depends on the keypoint's scale. The procedure returns the 4 similarities yielding the maximum number of inliers, with non-maxima suppression of overlapping detections (Figure 2, left).

Descriptor calculation. For each returned hypothesis, we compute a Boolean descriptor x ∈ {0, 1} K , where xi = 1 iff. ki has a matching point k ′ j under the transformation hypothesis (Figure 2, right).

IMPROVED SCORING OF HYPOTHESES

We now explain how to score each hypothesis according to its descriptor x.

Baseline

Assuming that all components of x are independent, the optimal score in the Bayesian sense is the ratio of posterior probabilities for the presence of the logo L and its absence ¬L:

p(L|x) p(¬L|x) ∝ Y i p(xi|L) p(xi|¬L) . (1)
Assuming uniform distributions for p(xi|L) and p(xi|¬L), and taking the logarithm, we obtain the number of inliers:

s bsl (x) = log p(L|x) p(¬L|x) ∝ X i xi. (2)
This scoring criterion is widely used in RANSAC implementations [START_REF] Chum | Total recall: Automatic query expansion with a generative feature model for object retrieval[END_REF][START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF][START_REF] Philbin | Descriptor learning for efficient retrieval[END_REF].

Burstiness handling

In image retrieval, it has been observed that the independence assumption on p(xi|L) and p(xi|¬L) does not hold. Visual elements can appear more frequently in an image than a statistically independent model would predict. This phenomenon, initially observed in text retrieval, is called burstiness [START_REF] Jégou | On the burstiness of visual elements[END_REF][START_REF] Slava | Distribution of content words and phrases in text and language modelling[END_REF][START_REF] Rasmus | Modeling word burstiness using the Dirichlet distribution[END_REF]. Jégou et al. [START_REF] Jégou | On the burstiness of visual elements[END_REF] use the following improved score in image retrieval. They group matches with identical visual word indexes and apply a down-weighting scheme (i.e. square-rooting) on the score of each group:

s burst (x) = h X g=1 s X i xiGi,g. (3)
G ∈ {0, 1} K×h is the "grouping matrix": Gi,j = 1 if ki is quantized into visual word j (G's rows sum to one) 1 . This score works well in practice because square-rooting approximates a probabilistic model of the visual word dependencies [START_REF] Cinbis | Image categorization using non-iid image models[END_REF].

Proposed weighting

Instead of grouping keypoints according to their visual word index, we propose to learn a new grouping matrix H that better describes the burstiness in logo detection. The intuition is to group keypoints frequently matched together in incorrect detections (Figure 3), so as to down-weight their score as a whole. To learn their statistical distribution, we first measure the correlations between keypoints matched in incorrect hypotheses. Then, clustering correlated matches yields H.

Let X = [x1 . . . xL] ∈ {0, 1} K×L the matrix of L descriptors collected on a training dataset disjoint from the test set, from which the logo is known to be absent. All these descriptors correspond to incorrect hypotheses (Figure 1). Let C ∈ R K×K + be the corresponding correlation matrix 2 ,

1 In [START_REF] Jégou | On the burstiness of visual elements[END_REF], the formula is expressed differently but it is equivalent for binary scores. 2 Strictly speaking, the columns of X should be centered, but the mean of each column is approximately 0, so it makes little difference.

with cells defined as:

Cij = cov(Xi, Xj) σX i σX j = P ℓ x iℓ x jℓ r `Pℓ x 2 iℓ ´"P ℓ x 2 jℓ " (4)
Elements of C range in [0, 1] and represent how frequently two keypoints jointly appears in X. To extract groups, we perform hierarchical clustering [START_REF] Sibson | SLINK: An optimally efficient algorithm for the single-link cluster method[END_REF] in the space of logo keypoints using C as similarity metric, see Figure 3 for example clusters. This yields the grouping matrix H. In our experiments, we have set the minimum linkage threshold [START_REF] Sibson | SLINK: An optimally efficient algorithm for the single-link cluster method[END_REF] to 0 (yielding the minimum number of clusters), which produced the best results. Hereafter, we also investigate more general weighting schemes and soft-assigned clustering.

Generalized burstiness. We generalize the burstiness formulation of eq. ( 3) into

scorr(x) = h X g=1 ω ([x1H1,g, . . . , xK HK,g]) , (5)
where ω : {0, 1} K → R + is down-weighting function. Among many possibilities, we experiment square-rooting ω(x) = pP i xi as in eq. ( 3), ω(x) = log(1 + P i xi), and the family of norms

ω(x) = x γ (6)
with γ > 1. All those functions obey the desirable property that ∀a1 . . . an ∈ {0, 1} n ,

X i ω(ai) ≥ ω([a1 . . . an]), (7)
which guarantees that a small number of points spread into a large number of groups score higher than the same number of points belonging to a single group.

Soft-assigned groups. Because soft-assignment has been

shown to be beneficial for image retrieval [START_REF] Philbin | Lost in quantization: Improving particular object retrieval in large scale image databases[END_REF], we evaluate a more flexible grouping procedure. We extend the Boolean matrix H to a real-valued matrix S ∈ [0, 1] K×h . The scoring function [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF] generalizes smoothly if we choose S with rows summing to one. Also, we simply assume that the number of groups h is equal to K: there is one "soft-group" per logo keypoint ki.

A natural choice for the matrix S is to build it directly from the correlation matrix C. We use the following formulation:

Si,j = C α i,j P j C α i,j . ( 8
)
The normalization is required so that the rows sum to 1. We raise the elements of C to a power 0 < α ≤ 1, to compensate for the fact that the correlation values are underestimated due to the asymmetric repeatability noise in keypoint extraction and matching. Typically, α = 0.4 yields good results (see Section 4).

In the experiments, we evaluate both the hard and softassigned grouping H and S, respectively, in addition to different down-weighting schemes ω.

EXPERIMENTS

In this section, we present experimental results on two different datasets introduced to evaluate logo retrieval in real-world images: the BelgaLogos [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF] and the FlickrLogos [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF] datasets.

BelgaLogos dataset

The BelgaLogos dataset [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF] was created in collaboration with the Belgavox press agency. It is composed of 10,000 press photographs annotated for 26 logos. The set Qset1 is composed of 55 queries, each defined by an image from the database and the logo's bounding box in this image. The set Qset2 is composed of 26 thumbnails, representing the "ideal" logos, see top row in Figure 4.

As in [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF], we use SIFT keypoints [START_REF] David | Distinctive image features from scale-invariant keypoints[END_REF] and report mean Average-Precision (mAP) results for both sets. We introduce several optional improvements: DBL -doubling the image size prior to keypoint extraction, because the logo often appears at small scales in database images; INV -generating color-inverted versions of monochromatic logos like Adidas or Nike, that can be printed either black-on-white or white-on-black; MA -assigning multiple visual words to each query keypoint [START_REF] Jégou | Hamming embedding and weak geometric consistency for large scale image search[END_REF].

For our approach, we estimate the correlation matrix C from an external dataset composed of 15,000 images collected from Flickr. We manually verified that no logo appears in them, using the baseline method and checking the 10 top retrievals for each logo.

FlickrLogos dataset

Romberg et al. [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF] built the FlickrLogos dataset by downloading real-world images including one of 32 logos from Flickr. The dataset is very different from BelgaLogos, because logos are more textured, usually large, and come in several versions. Sample images are shown in Figure 4. The dataset is partitioned into 3 subsets: P1 contains 10 images per logo, chosen to contain little clutter and noise. Each of P2 and P3 contains 30 images per logo and 3,000 additional background images where no logo appears. Following the protocol of [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF], P1 serves as query set and results are reported in terms of average recall on P3, while P2 is used as training set to learn per-logo recognition thresholds (the target precision is set to 95% as in [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF]).

In order to estimate C and learn the recognition thresh- s bsl s burst scorr proposed method 0.577 0.614 0.726 state of the art [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF] 0.61 olds at the same time, we perform 10-fold cross-validation on P2.

Impact of the parameters

Table 1 compares different down-weighting schemes ω that were proposed in [START_REF] Jégou | On the burstiness of visual elements[END_REF] (#1, #2) and in this paper (#3), eq. ( 6). The proposed scheme (#3) depends upon parameters γ, plus α in the case of soft-assignment (eq. ( 8)). To evaluate their impact, we fix in turn γ and α to their optimal values and vary the other parameter. Overall, all examined values of γ and α yield good performance, with better results for larger values of γ and for α ∈ [0.3, 0.5]. In the following, we thus set γ = ∞ (i.e. ω(x) = x ∞ = max(x)) and α = 0.4.

In Table 1, we evaluate different weighting schemes ω for our scoring method scorr. The proposed weighting (#3) yields excellent performance: it counts the number of detected groups rather than the number of inlier points. In particular, it clearly outperforms the other ones for all datasets when combined with soft-assignment. It is also worth noting that the logarithmic down-weighting scheme (#2) only performs well for the case of hard-assignment. In the following, we thus use ω(x) = x ∞ and soft-assignment.

Quantitative results and conclusion

Table 2 shows quantitative results for the BelgaLogos dataset. The proposed scoring constantly improves over the baseline or the classical burstiness model of [START_REF] Jégou | On the burstiness of visual elements[END_REF]. The latter approach improves slightly by about 2%. In contrast, the proposed learned burstiness model outperforms the baseline by up to 9% and 10% on Qset1 and Qset2 for the best setting, i.e. with DBL, INV and MA. Due to the specificity of the BelgaLogos dataset, in which logos often appears at very small scales and with inverted colors (Figure 5), DBL, INV and MA are also important to obtain good performance. On FlickrLogos (Table 3), the proposed method largely improves over the baseline (+15%) and over the classical burstiness model [START_REF] Jégou | On the burstiness of visual elements[END_REF] (+3.7%).

Compared to the state of the art, our approach outperforms the a-contrario query expansion of Joly et al. [START_REF] Joly | Logo retrieval with a contrario visual query expansion[END_REF] by 7% on Qset1 and 22% on Qset2. In contrast to this method, our approach does not add any overhead at run-time and could, in fact, be combined with query expansion. On the Flick-erLogos dataset, we outperform the method of [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF] by 12%. Contrary to the cascaded triplets features of [START_REF] Romberg | Scalable logo recognition in real-world images[END_REF], our approach remains generic and can be added to existing BoWbased systems.

Figure 2 :

 2 Figure 2: Left: Detecting logo hypotheses. Right: construction of the Boolean descriptor x for a detection hypothesis.

Figure 3 :

 3 Figure 3: Groups of keypoints (shown as ellipses) obtained after correlation-based clustering (section 3.3).

Figure 4 :

 4 Figure 4: Sample logos from BelgaLogos (set #2) [7] (top row) and from Flickr-Logos [16] (bottom row).

Figure 5 :

 5 Figure 5: Examples of detections output by our system for the 3 queries on the left.

Table 1 :

 1 Comparison of the different downweighting schemes, for hard and soft assignment.

		ω(x)	Belga Qset1 Belga Qset2 Flickr
		None	0.328	0.383	0.577
	Hard H	#1 #2 log (1 + pP #3 x γ , γ = ∞ xi P xi)	0.387 0.385 0.388	0.440 0.440 0.428	0.658 0.656 0.724
	Soft S	#1 #2 log (1 + pP #3 x γ , γ = ∞ xi P xi)	left-out because (7) does not hold 0.333 0.389 0.566 0.414 0.481 0.726

Table 2 :

 2 mAP performance for various options on the BelgaLogos dataset (Section 4.1).

			Qset1			Qset2
	options	sbsl	s burst	scorr	sbsl	s burst	scorr
	default	0.205 0.216 0.250 0.182 0.215 0.250
	+DBL	0.286 0.302 0.367 0.206 0.226 0.283
	+INV	-		-	0.300 0.329 0.388
	+MA	0.328 0.344 0.414 0.383 0.404 0.481
	state of the art [7]		0.341			0.257

Table 3 :

 3 Recall performance on the FlickrLogos dataset (precision not shown but always ≥ 98%).

Acknowledgments

This work was realized as part of the Quaero Project, funded by OSEO, French State agency for innovation.


IBME Louis Mayaud
email: louis.mayaud@eng.ox.ac.uk

Marco Congedo
email: marco.congedo@gmail.com

Sabine +33476826252

Filipe

Guillaume Charvet
email: guillaume.charvet@cea.fr

Remy +33438780144

ROBUST VIRTUAL KEYBOARD FOR BRAIN-COMPUTER INTERFACE (ROBIK): AN HALFWAY UPDATE ON THE PROJECT

ABSTRACT

Introduction The principle of a Brain-Computer Interface or BCI is to control a device through the extraction and interpretation of signal features from electroencephalograms (EEG) collected either from the surface of the scalp or through invasive measurements. This late idea of communication technique [START_REF] Vidal | Toward direct brain-computer communication[END_REF], offers the advantage of bypassing the need for muscle activity in the control chain and is therefore presented as a promising alternative to restore communication and control in severely disabled patients (Wolpow, et al. 2002). However, the lack of robustness and ergonomics of both available software and EEG measurement techniques have delayed the transfer of this technology to patients bedsides.

Materials and methods The French Research Agency has funded a 3-year project gathering national leaders in microelectronics (CEA-Leti), EEG signal processing (Gipsa-Lab) and clinical management of severely disabled people (Raymond Poincar hospital). The aim of the project is the development and the clinical validation of a Brain-Computer Interface prototype for communication. As an initial step, a survey was carried out to assess patients' and users (family and caretakers) needs, which were translated into specifications, on the basis of which software and hardware were developed.

Results The survey (n=45) highlighted the need for easy-to-setup systems (installation time ¡ 15min), which stresses the importance of mechanical comfort and customization of application. The development of signal processing techniques has led to improvements of the P3Speller paradigm. A first prototype of a 32-channel EEG recording system is under development. To ease the EEG measurements and reduce installation time, the system has a reduced size. It includes the analog amplification and digital conversion of 32 channels sampled at 1 kHz, as well as the wireless data transmission to a computer. First in vivo validations were performed on small animals. This system will be optimized and connected to a headset specifically designed to provide a comfortable and handy interface with dry electrodes.

Discussion The present project will still run for one and a half years, ending with its clinical validation in a population of severely disabled patients, which will compare performances of the system with existing assistive technologies. At this stage, the proposed system yields very promising results, and outperforms the current state-of-theart. If such a system is shown to perform better than current users assistive technology, it could reach the commercial availability for severely disabled patients within the next 5 years.

KEY WORDS

Brain-Computer Interface (BCI), electroencephalography (EEG), Evoked Potentials, dry electrodes, assistive technology, P300 speller

Introduction

The principle of a Brain-Computer Interface (BCI) is to control a device through the extraction and interpretation of signal features from electroencephalograms (EEG) collected either from the surface of the scalp or through invasive measurements. This old idea of communication technique [START_REF] Vidal | Toward direct brain-computer communication[END_REF], offers the advantage of bypassing the need for muscle activity in the chain of control, and is therefore presented as a promising alternative to restore communication and control in severely disabled people diseases [START_REF] Birbaumer | Breaking the silence: brain-computer interfaces (bci) for communication and motor control[END_REF][START_REF] Birbaumer | Brain-computer interfaces: communication and restoration of movement in paralysis[END_REF].

The use of BCI in a severely disabled population is however less frequent than commonly believed. Firstly, a literature research on BCI shows that a third of published papers deal with algorithms aiming to provide higher and more robust overall performances. Unfortunately, only a few algorithms are made available on one of the two open-source platforms for BCI [START_REF] Schalk | Bci2000: a general-purpose brain-computer interface (bci) system[END_REF][START_REF] Renard | Openvibe: An open-source software platform to design, test, and use brain-computer interfaces in real and virtual environments[END_REF]. In fact, only two implementations of the well-known and broadly-published P3Speller [START_REF] Farwell | Talking off the top of your head: toward a mental prosthesis utilizing event-related brain potentials[END_REF] make use of these platforms. Secondly, nearly a quarter of the publications focus on neurophysiology and the use of EEG features with no major breakthrough in this direction, as BCI mainly makes use of late paradigms: Event Related Potentials (ERP) [START_REF] Halder | An auditory oddball brain-computer interface for binary choices[END_REF][START_REF] Citi | P300based bci mouse with genetically-optimized analogue control[END_REF], Steady State Visual Evoked Potentials (SSVEP) [START_REF] Muller-Putz | Control of an electrical prosthesis with an ssvep-based bci[END_REF], motor rhythms [START_REF] Bai | A high performance sensorimotor beta rhythm-based brain-computer interface associated with human natural motor behavior[END_REF] and Local Field Potentials (LFP) [START_REF] Moran | Evolution of brain-computer interface: action potentials, local field potentials and electrocorticograms[END_REF]. Finally, only a fifth and a tenth of publications deal with specific applications and evaluation in patients, respectively. Before 2005, the proportion of medical publications in the field was around 10%, a figure that nearly doubled since then.

The Clinical Investigation Centre for Innovative Technology (CIC-IT) at the Raymond Poincar hospital (Paris hospital trust -APHP) is the leading French multidisciplinary research group in the assessment of technology for use by disabled patients. BCI was lately identified as a potential assistive technology, despite being regarded with much skepticism by clinicians, due to the disadvantageous usability-performance ratio against existing assistive technologies. Further examination of the reasons why BCI are not being widely available led to the following conclusions:

• Lack of robustness of the application: no real adaptive signal processing tools taking into account the interand intra-individual variations of EEG features and high sensitivity to all sources of noise ;

• Lack of ergonomics: the setup procedures of EEG recording device is particularly tedious. Moreover, the software calibration procedures (primarily designed for researchers) are complex and are not integrated with end-user application such as emails or Internet.

A consortium was built and a project proposal developed whose aim is to investigate these issues, while working on hardware and software simultaneously, in close relation with patients and medical teams specifications. Finally, the French Research Agency funded a 3year project regrouping national leaders in microelectronics (CEA-Leti), EEG signal processing (Gipsa-Lab) and clinical management of severely disabled people (Raymond Poincar hospital) together with an industrial leader in invasive brain signal monitoring (DIXImicrotechniques). The project aims to develop and validate a BCI prototype for communication. As an initial step, a survey was carried out to assess patients and users (family and caretakers) needs, which were translated into specifications, subsequently used as guidelines for both software and hardware development.

Materials and methods

Patients need and technological locks

A questionnaire was written by occupational therapists specialized in providing effective assistive technology to a The aim of the clinical study is to assess the feasibility of a state-of-the-art BCI for communication described in [START_REF] Krusienski | Toward enhanced p300 speller performance[END_REF]. 24 cupule-electrodes were fitted to 10/20 standard locations by an EEG technician and signals were recorded at 256 Hz with a Porti32 T M from TMSi (Twente, Netherlands). The Porti32 T M was connected to the open-source platform OpenViBE. Signals were filtered using the xDAWN spatial filter [START_REF] Rivet | xdawn algorithm to enhance evoked potentials: application to brain-computer interface[END_REF] prior to SVM voting classifier [START_REF] Rakotomamonjy | Bci competition iii: Dataset ii-ensemble of svms for bci p300 speller[END_REF].

New P300 paradigm for daily use

Our cohort consists of 20 healthy volunteers who gave their informed consent. For each of them, 32-EEG signals were recorded using a Nexus32 T M from MindMedia (Roarmond-Herten, Netherlands). Recordings took place in a room specifically labeled by the national health and safety regulatory committee (AFSSAPS). Subjects were given instructions for a standard P3Speller experiment [START_REF] Krusienski | Toward enhanced p300 speller performance[END_REF]. Stimulation and signal processing were handled by Open-ViBE. Signals were analyzed with the aim of optimizing:

• the number and locations of sensors ;

• pre-processing and classification parameters.

Wireless EEG headset for daily use at home

A wireless multichannel data acquisition electronic system is being designed for electroencephalograms (EEG) recording. This system is based on the EleCtro encephalogramme Recording Interface for Neuroactivity Studies (ECRINS) prototype which is shown on the Figure 1. This electronic system is consists of a custom integrated circuit (ASIC) for signal conditioning, amplification and digitization, as well as a commercial module for RF transmission. The RF transmission module is capable of handling a 32-channel EEG recording sampled at 1 kHz with 12-bit resolution. This integration is the first step towards a lightweight EEG headset for Brain Computer Interface (BCI) studies.

In parallel, an extensive literature review for EEG sensors was carried out. Three dry-electrodes solutions for easy-to-setup BCI headsets were identified: micro-array electrodes [START_REF] Hutton | Diamond microarray electrodes[END_REF], comb electrodes [START_REF] Pradeep | Dry electrodes for electroencephalography[END_REF][START_REF] Zander | Frontiers: A dry eeg-system for scientific research and braincomputer interfaces[END_REF] and hybrid electrodes [START_REF] Estepp | Validation of a dry electrode system for eeg[END_REF]. Commercially available solutions are currently evaluated together with a prototype developed by the projects industrial partner. Eventually the selected candidate will be interfaced with an EEG headset designed to comply with specifications. To date, several headset designs were suggested and are currently prototyped.

Results

Patients need and technological locks

The survey (n=40) highlighted the need for easy-to-setup systems (installation time shorter than 15 min in 82% of responses) and revealed that the expected daily time of use would be higher than 2 hours. This, in particular, stresses the importance of mechanical comfort (selected as first priority for 72% of users and 60% of users) and customization of the application (second priority for 79% of users and 62% of medical staff). In terms of mechanical comfort, the main point was the compatibility with potential headrests and most cited applications were, in order of importance: voice synthesizer, web surfing, emails and domotics interfaces.

At the time of writing 4 ICU patients have been included in the feasibility study (2 male, 2 female, age range: 22 to 63 y.o.). Two patients presented with quadriplegia (C1, C4-C5), one with Guillain-Barre syndrome, and one was a post-surgery patient. Three patients were mechanically ventilated through translaryngeal tracheotomy (75%). Even though the BCI spelling system was calibrated at a spelling rate of 1.8 characters per minute, post-analysis showed that two patients (50%) had sufficient Signal-to-Noise-Ratio (SNR) to achieved spelling rate comparable to those found in literature on healthy subjects. The other two patients withdrew from the study, one because of respiratory distress and the other because of an impossibility to set up occipital electrodes, respectively.

Robust P300 paradigm

Identification of optimal subset of sensors was carried out by comparing P300 identification performance from three different group of sensors: the complete EEG headset (32 electrodes), 8 electrodes individually selected to provide best performance, and common subgroup of 8 electrodes providing the best average performance in the population of 20 healthy volunteers [START_REF] Cecotti | Suboptimal sensor subset evaluation in a p300 brain-computer interface[END_REF][START_REF] Rivet | Eeg sensor selection by sparse spatial filtering in p300 speller brain-computer interface[END_REF][START_REF] Cecotti | A robust sensor-selection method for p300 brain-computer interfaces[END_REF]. This study, illustrated in Figure 2, shows that an individually customized subsets of 8 leads provides similar accuracy (Acc) than the complete headset (32 electrodes) in 85% of subjects. The subset of 8 electrodes (out of 32) common to the 20 subjects also provides equivalent performances, except for Optimization of the time window for P300 detection [START_REF] Cecotti | Impact of the time segment analysis for p300 detection with spatial filtering[END_REF] and spatio-temporal feature selection [START_REF] Jrad | Spatio-temporal feature extraction and classification of event-related potentials[END_REF] led to an average accuracy greater than 95% amongst the 20 subjects. In particular, partial use of training data provided an average accuracy of 80% with as little as 2 min signal for calibration [START_REF] Cecotti | Reducing calibration time for the p300 brain-computer interface speller[END_REF]. This particularly meets the needs expressed by clinicians and patients for a short setup time.

Finally, the data set was used to understand better the implications of close target stimuli. Analysis of such events revealed that in order to elicit an unaltered P300, target stimuli should be separated by at least 450 ms. A novel pseudo-random paradigm preventing such events, led to an estimated 30% relative increase in bit-rate [START_REF] Cecotti | A pilot study for improving the graphical user interface of p300 based bcis[END_REF].

Wireless EEG headset for daily use at home

A prototype of the electronic board was designed. Thanks to the experimental setup based on the ECRINS board shown on Figure 3, the input referred noise of the recording system (based on a dedicated ASIC) and preliminary throughput tests have been performed. These first tests confirmed the feasibility of transmitting 32 channels with a 12-bit resolution sampled at 1 kHz. An input-referred RMS noise in BW1 [0,5Hz; 300Hz] of about 1µV for each channel was measured. A software for wireless recording and visualization of EEG signals was developed, and the complete system was validated in vivo on a rodent model, as described in Figure 1. The experimental procedures and animal care were carried out in compliance with the European Community Council Directive of 24th November 1986 (86/609/EEC).

Literature on dry electrode technology discarded micro-array technology for lack of sufficient evidence, and health and safety issues associated with minimally invasive modalities over time. Similarly, hybrid capacitive technology was discarded because of the necessary additional electronic stages, leaving comb designs the preferable technology for BCI.

Conlusion

The current project will run for another year and a half, and should end with the clinical validation of the prototype in a population of severely disabled patients. The study will compare performances of the system with existing assistive technology. Major remaining steps involve the break down of line-columns stimuli in the P3Speller paradigm, development of an optimized version of the electronic board. The board will be connected to a headset specifically designed to provide comfortable and practical interface with dry electrodes.

The major advantage of the approach described here is the central role played by the clinical partner: from the original project design phase to the clinical validation, developments are led by pragmatic recommendations from practitioners and patients, allowing the know-how of engineers and researchers to benefit patients directly. Such a collaboration allowed us to address the issues encountered during the feasibility study at an early stage, and should remain the driving force behind this project.

In addition to the hardware and algorithms (available on the OpenViBE forge) developed in the scope of the project, benefits for the BCI community involve a first evaluation of patients' needs as well as patients' data acquired in Intensive Care Unit during a P300 speller protocol. Such an evaluation of patients' need and in-situ data are extremely rare in the BCI community that works mainly with healthy volunteers and is therefore highly valuable material to guide future research in the field.

Finally, if the system is proven to perform better than current users assistive technology in the clinical validation stage, commercial availability for severely disabled patients could be achieved within the next 5 years.

Figure 1 .

 1 Figure 1. Picture of the ECRINS electronic board with its RF antenna used for the RoBIK project.

Figure 2 .

 2 Figure 2. Description of classification accuracy for n=20 subjects + average (right) for all 32 sensors, an optimized subset of 8 sensors for all subjects and an individually optimized subset of sensors.

Figure 3 .

 3 Figure 3. Description of the technical setting for the in vivo validation of the electronic card and the wireless link.

Acknowledgement

We would like to thank the National Research Agency (ANR) and the General Directorate for Armament (DGA) for funding this project (Project ANR-09-TECS-013-01-RoBIK). We would like to thank the medical staff who have been (and still are) involved on this project, and in particular: Aurelien Van Laghenhove, Michelle Lejaille, Marjorie Figère, Marine Garguillo from the CICIT; Jean-Marie, Gilles and Caroline together with Dr. Azabou and Prof. Lofaso and Hérault from the Functional Exploration Unit, Justine Bouteille and Samuel Pouplin from the New Technologies Platform (PFNT), as well as many others involved on the project.