Xavier Blanc
email: blanc@ann.jussieu.fr

Benjamin-Edouard Peigney

Homogenization of heat diffusion in a cracked medium

HAL is

Dans cette Note, on traite le cas de la diffusion linéaire dans un milieu périodique fracturé. Nous montrons que l'effet des fractures sur l'énergie du système peut être modélisé par un terme source en volume au sein du milieu homogénéisé. L'analyse asymptotique du problème de départ (1)-( 2) conduit à la formulation du modèle homogénéisé (10)-(11) pour lequel il est nécessaire d'introduire des conditions de transmission à l'interface correspondant à la singularité en pointe des fissures. Nous montrons en particulier que la solution u, et son gradient ∂ n u, y subissent un saut.

Nous mettons au point une approche par point fixe consistant à résoudre successivement (13-gauche) dans le sous-domaine Ω + intact , et (13-droite) dans le sous-domaine Ω -contenant les fractures (voir figure [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF]). Nous montrons ainsi que le problème (10)-(11) est bien posé. Cela nous donne au passage une méthode de construction de la solution, que nous exploiterons à des fins numériques.

Nous énonçons alors la proposition 2.2 où l'on établit que la solution du problème exact (1)-( 2) défini dans la géométrie fracturée Ω ε converge faiblement vers la solution de (10)-(11) défini dans le domaine homogénéisé Ω, ne contenant plus la description des fractures (voir figure [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF]). La preuve rigoureuse de la proposition 2.2 sera exposée dans une publication ultérieure [START_REF] Blanc | Homogenization of heat diffusion in a cracked medium[END_REF].

Enfin, nous appliquons la méthode du point fixe précédemment décrite pour résoudre numériquement (10)-(11). Nous développons une approche équivalente en écrivant la formulation faible du problème écrite dans tout le domaine homogénéisé conduisant à l'équation (14), caractérisée par la présence d'une masse de Dirac localisée en pointe de fissure. Ces deux méthodes donnent des résultats cohérents avec le calcul direct de la fracture.

Motivation and setting of the problem

We consider the propagation of heat through a cracked medium, exposed to an incoming energy flux. Physically, the exchange surface between the medium and the source may be greatly modified by the fractures. This may have a significant impact on the energy balance of the considered system. In many situations, the geometry of the cracked media is too intricate to be described precisely. Thus, we cannot model the surface of the cracked medium directly. Besides, the shape of the fractures may have a stochastic feature and it may involve many spatial scales. Full numerical simulations of such multi-scaled media become hence infeasible.

That is why we have been looking for an average approach, to capture the effects of cracks in a homogenized medium. The model presented here is simple enough to be coupled to standard Finite Element codes. The physical idea behind the method developed in this Note, called "MOSAIC" (Method Of Sinks Averaging Inhomogeneous behavior of Cracked media), is to treat the flux enhancement induced by the crack as a volume source term in the homogenized energy equation. We will show that this can be justified rigorously by homogenization theory.

For the sake of simplicity, we shall assume a linear behavior law of the material 1 . The linear diffusion problem can thus be modeled by:

                   ∂ t u ε -∆u ε = 0 in Ω ε , ∂ n u ε = 0 on Γ 0 ε ∂ n u ε = 1 on Γ 1 ε , ∂ n u ε = α -β 2 ε on Γ α ε , ∂ n u ε = β α on Γ β ε . (1)
We also need an initial condition:

u ε (x, y, t = 0) = u 0 (x, y), (2)
so that problem (1)-( 2) is well-posed. We impose u ε (x, y, t) to be periodic of period ε with respect to the variable y. The domain Ω ε , as well as the boundaries Γ 0 ε , Γ 1 ε , Γ α ε , Γ β ε are defined on Figure 1. As suggested by the red lines, the left part of the figure is a zoom in the y variable of the small shaded area in Ω.

The period ε > 0 is supposed to be small and will tend to 0, whereas α ∈ [0, 1) is a fixed parameter related to the width of the crack. The parameter β ∈ [0, α) measures the portion of the flux which, coming through the segment {x = -1, -αε/2 < y < αε/2}, reaches the bottom Γ β ε of the crack. The remaining part of the incoming flux is distributed on the horizontal part of the boundary, namely Γ α ε . The parameter β is supposed to be fixed. The boundary conditions in [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF] are defined in such a way that the total incoming flux is exactly equal to 1, which is the value we impose on the left boundary in the case α = 0 (no crack). A space-time dependence of the flux applied on the boundaries Γ α ε may be introduced but it does not affect the homogenization process that we describe here.

Γ ε 1 Γ ε 1 Γ α ε Γ α ε Γ ε β εα Γ β Γ 1 Γ 0 ε Ω Γ ε 0 Ω - Ω + (0,0) x y ε 1
Similar homogeneization problems have already been tackled in [START_REF] Achdou | Effective boundary conditions for laminar flows over periodic rough boundaries[END_REF], [START_REF] Amirat | Effective boundary condition for Stokes flow over a very rough surface[END_REF] or [START_REF] Dalibard | Effective boundary condition at a rough surface starting from a slip condition[END_REF] but the geometry, the equation type as well as the boudary conditions were not the same as the one considered here.

Asymptotic expansion and homogenized equation

To carry out an asymptotic expansion of the solution u ε of (1) in powers of ε, we "scale" the variable y, in the spirit of [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF]. Actually, 2 scales describe the model: the variable y is the macroscopic one, whereas y ε represents the "microscopic geometry". Thus, we define:

u ε (x, y, t) = v ε x, y ε , t ,
so that v ε is periodic of period 1 in y. Applying the change of variable in (1), we can write the system satisfied by v ε :

                     -∂ 2 x v ε - 1 ε 2 ∂ 2 y v ε + ∂ t v ε = 0 in Ω 1 , ∂ n v ε = 0 on Γ 0 1 , ∂ n v ε = 1 on Γ 1 1 , 1 ε 2 ∂ n v ε = α -β 2 on Γ α 1 , ∂ n v ε = β α on Γ β 1 . (3)
Firstly, we notice that in the system (3) the domain Ω 1 does not depend on ε anymore. We have to study an equation depending on ε in a fixed domain. Secondly, the parameter ε appears in the equation (3) only as ε 2 , which means that ε 2 is a good parameter for an asymptotic expansion. Thus, it seems natural to look for v ε as follows:

v ε (x, y, t) = v 0 (x, y, t) + ε 2 v 1 (x, y, t) + ε 4 v 2 (x, y, t) + . . . (4)
Hence, we insert the ansatz (4) into the system (3) and identify the different powers of ε 2 . We obtain:

-At the order ε -2 : ∂ 2 y v 0 = 0, and the condition of periodicity in y verified by v: v 0 (x, 1/2, t) = v 0 (x, -1/2, t), implies that v 0 does not depend on y :

v 0 (x, y, t) = v 0 (x, t). (5)
-At the order ε 0 :

-∂ 2 x v 0 -∂ 2 y v 1 + ∂ t v 0 = 0. ( 6
)
The boundary conditions on v 1 give

∂ n v 1 = 0 on Γ 0 1 , ∂ n v 1 = 0 on Γ 1 1 , ∂ n v 1 = α -β 2 on Γ α 1 , ∂ n v 1 = 0 on Γ β 1 .
Integrating ( 6) with respect to y and using the boundary value for ∂ n v 1 as well as the periodicity in y, we get:

       -∂ 2 x v 0 + ∂ t v 0 = α -β 1 -α in Ω 1 ∩ {x < 0}, ∂ n v 0 = 1 on Γ 1 1 , ∂ n v 0 = 0 on Γ α 1 ,        -∂ 2 x v 0 + ∂ t v 0 = 0 in Ω 1 ∩ {x > 0}, ∂ n v 0 = 0 on Γ 0 1 , ∂ n v 0 = β α on Γ β 1 . (7)
This system is not well-posed, since boundary conditions are missing at the interface {x = 0}. Actually, we can show that the flux ∂ x v 0 is not continuous accross the interface {x = 0}, so that we need to introduce a boundary layer at this interface. This is done by changing variables once again and defining v ε (x, y, t) = w ε x ε , y, t . This function w ε is thus defined on the set

Ω ε := (x, y) ∈ - 1 ε , 1 ε × - 1 2 , 1 2 , |y| > α 2 if x < 0 .
On this set, w ε satisfies the equation

ε 2 ∂ t w ε -∂ 2 x w ε -∂ 2 y w ε = 0,
with the boundary conditions ∂ x w ε = -εβ α on {x = 0} ∩ {|y| < α/2}. Integrating this equation in the domain Ω ε ∩ {|x| < δ} , letting ε, then δ, go to zero, we find that lim

ε→0 {x=0 + } ∂ x w ε = lim ε→0 {x=0 -,α/2<|y|<1/2} ∂ x w ε -εβ .
Recalling the link between w ε and v ε , this leads to the so-called transmission conditions (8):

v 0 (x = 0 -) = v 0 (x = 0 + ), (1 -α)∂ x v 0 (x = 0 -) = ∂ x v 0 (x = 0 + ) + β. ( 8
) on {x = 0} \ Γ β 1 = {x = 0} ∩ {|y| > α/2}. Making use of: v y ε * -⇀ 1/2 -1/2 v(y)dy, (9)
in L ∞ for any periodic function v, together with the fact that

u ε (x, y, t) ≈ v 0 x, y ε , t ,
we find that, in the limit ε → 0, u ε is well approximated (in the weak sense) by the solution u of

     -∆u + ∂ t u = α -β in Ω -, ∂ n u = 1 -α on Γ 1 , u is 1 -periodic in y.      -∆u + ∂ t u = 0 in Ω + , ∂ n u = 0 on Γ 0 , u is 1 -periodic in y. ( 10
)
We also have the corresponding transmission condition inherited from (8):

u(x = 0 -) = (1 -α)u(x = 0 + ), ∂ x u(x = 0 -) = ∂ x u(x = 0 + ) + β. ( 11
)
u evolves in the homogenized domain Ω = Ω + ∪ Ω -with the boundaries Γ 1 , Γ β , Γ 0 defined on Figure 1.

Remark 1 Of course, if u 0 does not depend on y, the solution of (10)-( 11) is constant with respect to y. However, it is also possible to apply an incoming flux wich depends on y, that is, impose ∂ n u ε = F (y) on Γ 1 ε , in system (1). In such a case, the above analysis is still valid, and the solution a priori depends on y, even if u 0 does not. Proposition 2.1 Assume that α ∈ (0, 1). Then for any T > 0, problem (10)-( 11) has a unique solution (u -, u + ) ∈ X, where

X = C 0 [0, T ], H 1 (Ω -) ∩ C 1 [0, T ], L 2 (Ω -) × C 0 [0, T ], H 1 (Ω + ) ∩ C 1 [0, T ], L 2 (Ω + ) .
Proposition 2.2 Let u ε be the unique solution to [START_REF] Bensoussan | Asymptotic analysis of periodic structures[END_REF]. We extend it by 0 outside Ω ε , and assume that the initial data u ε (t = 0) is such that u ε (t = 0) -⇀ u 0 in L 2 (Ω). Then, for any T > 0, we have

u ε -⇀ ε→0 u in L 2 (Ω × [0, T ]), ( 12
)
where u is the unique solution to (10)-(11).

Building of the solution

Fixed-point approach

To compute the solution of the homogenized problem (10)-(11) , we firstly devise a fixed-point approach connecting the sub-domains Ω + and Ω -:

               -∆u + ∂ t u = 0 in Ω + , ∂ n u = 0 on Γ 0 , ∂ n u = F on Γ β , u is 1 -periodic in y, u(t = 0) = u 0 in Ω + .                -∆u + ∂ t u = α -β in Ω -, ∂ n u = 1 -α on Γ 1 , u = g on Γ β , u is 1 -periodic in y, u(t = 0) = u 0 in Ω -. (13)
We start with an initial guess for the flux F 0 ∈ L 2 (Γ β ), to which we associate the solution u +,0 of (13-left). Then, define g 0 as the trace of (1 -α)u +,0 on Γ β , and solve (13-right) with data g = g 0 : this defines u -,0 , and a new flux F 1 = β -∂ n u -,0 in L 2 (Γ β ). Repeating this procedure, we build a sequence converging towards the solution of (10)-(11). This scheme is actually used to prove Proposition 2.1. Figure 2. Comparison of the direct and homogenized approaches for α = 0.1.The red curve is the direct calculation of the crack, i.e the limit solution of (1) as ε → 0. The two other curves correspond to the solution of the homogenized problem (10) computed by two approaches: in black, it is the fixed-point method whereas the blue curve represents the numerical solution of the weak formulation (14) involving a Dirac mass at x = 0. Computations presented here are performed using a P 1 finite element approximation on triangular meshes. It has been implemented using the software FreeFem++ [START_REF] Hecht | FreeFem++ (manual[END_REF].

Weak formulation approach

Another possible way to solve (10)-( 11) is to compute the weak formulation of the problem in the whole domain Ω corresponding to equation ( 14

A Dirac mass at {x = 0} has appeared in (14). This represents the singularity at the bottom of the crack.

Those two methods are illustrated on figure (2) and compared with the direct calculation. We can note that the fixed-point method appears to be more accurate than the solution of the weak formulation (14) especially as α increases. This is due to the fact that (14) involves a Dirac mass at x = 0, proportional to the width of the crack α. This term is treated approximately in our finite element simulation and leads to more significant errors for greater α. In some way, the fixed-point method amounts to treat the Dirac mass at x = 0 exactly.

1 (Figure 1 .

 11 Figure 1. The cracked domain Ωε and the homogenized domain Ω = Ω + ∪ Ω -.

  Spatial profile at time t 1 = 0.5.

-

  ∆u + ∂ t u = (α -β)1 {x<0} -α∂ x u(x = 0 + )δ x=0 + βδ x=0 in Ω, ∂ n u = 1 -α on Γ 1 , ∂ n u = 0on Γ 0 , u est 1 -periodic in y.


Axel Hutt
email: axel.hutt@inria.fr

Meysam Hashemi
email: meysam.hashemi@inria.fr

Peter Beim Graben

Peter Beim Graben
email: peter.beim.graben@hu-berlin.de

How to render Neural Fields more realistic

Complex systems are omnipresent in nature. They exhibit multiple spatial and temporal scales whose understanding and control is one of the most challenging problems of research for the last centuries. The notion of complexity is not welldefined [1], but typically a system is said to be complex if, in addition to the multiple scales, there exist self-organised sub-units in the system which are generated by smaller units. In many biological systems, climbing this hierarchy from smaller units to larger units implies increasing the spatial and temporal scale.

Before discussing neural systems, we would like to clarify a major scientific problem in this context: the choice of the description level. Let us assume the task to describe a water wave mathematically, i.e., derive a mathematical model whose solution describes the spatiotemporal phenomenon of a wave. Some researchers may ask: we want to understand how water works, and especially, how a wave is generated. Hence, a straightforward simple and fundamental approach is to take a closer look at the building blocks of water (remark : this is very similar to the way how biologists work in neuroscience today). The researchers find H 2 O-molecules, study the properties of the atoms and the inter-molecule interactions like hydrogen bridges and van-der-Waals bounds. Now, to describe a water wave, one has to study millions and millions of these molecules and their interactions.

take the reasonable approach to start simulating two, three, then and even hundreds of these molecules. The numerical analysis is demanding, but still there is no description of the water waves, since the system appears to be so complex taking into account all the different interactions of molecules and the research in this hard task got stuck in some way. Of course it is well-known that it is not necessary to study single molecules to describe water waves, today we take the famous Navier-Stokes equation (NSE) which includes the solution for this phenomenon. This equation involves mean fluid properties, e.g., inner friction and viscosity [START_REF] Landau | Fluid Mechanics[END_REF]. Hence the NSE considers average properties of interacting single molecules. In other words, it does not know single molecules, but provides a powerful description of a large ensemble or mass of molecules. It allows to describe several different rather complex fluid phenomena, but of course no phenomena related to single molecules. Hence the NSE equation provides a very good mathematical description of the system at a macroscopic description level only. Going back to the intended description of water waves by a single molecule study, this approach is not reasonable, probably it will not lead to a good description level of macroscopic phenomena and hence it is not constructive.

In todays' neuroscience, the approach of linking single neuron activity to macroscopic phenomena is attractive, e.g. in the context of cognition [START_REF] Gross | [END_REF][4][5], sleep [6] or anaesthesia [7,8]. These studies state a relationship between the single neuron activity (microscopic scale) and behavioral phenomenon (macroscopic scale). However, to our best knowledge it is not understood how the different experimental findings on different scales are linked to each other, i.e., the link between the two scales is not understood and no model linking the scales has been developed yet. This situation resembles closely the water wave task described above: it is clear that there is a relation between small sub-units (molecules or neurons) to large complex systems units (water wave or cognition), since the dynamics of the large units is generated by the sub-units, but a link appears to be too complex. Consequently, learning from physics and the NSE, it is necessary to consider more abstract, intermediate models whose elements are based on small sub-units properties but which allow to model large unit phenomena. In other words, it is much more effective to consider mesoscopic population models involving average properties of interacting neurons and which allow to describe macroscopic experimental phenomena, such as Local Field Potentials (LFP), encephalographic activity (EEG/MEG) or even behaviour. Promising candidates for such models are neural mass or neural field models. The present book chapter discusses recent advances in these models rendering the standard neural population models more realistic.

The subsequent sections do not give a complete overview over the recent advances in the field, a recent excellent review article already provides this information [9]. The present chapter first introduces briefly to two types of neural field models. Then it gives some details of few selected extensions, always introducing the neuroscientific problem by experimental data before presenting a mathematical description of the phenomena described.

2 Two classes of neural field models

Amari model

The first neural field models developed by Wilson and Cowan [10] and Amari [11] are continuum limits of large-scale neural networks. Typically, their dynamic variables describe either mean voltage [11] or mean firing rate [10,12] of a population element of neural tissue, see also the excellent review article of Bressloff [9]. In some subsequent sections we consider the paradigmatic Amari equation [11] describing the spatiotemporal dynamics of mean potential V (x,t) over a cortical d-

dimensional manifold Ω ⊂ R d : τ ∂V (x,t) ∂t = -V (x,t) + Ω K(x, y)S[V (y,t)] dy + I(x,t) . (1)
with the spatial synaptic kernel K(x, y) which defines the connectivity between site y ∈ Ω and site x ∈ Ω . The transfer function S is nonlinear and typically of sigmoidal shape. This model considers external inputs I(x,t), e.g., originating from other extracortical populations and from external stimulation. The model ( 1) takes into account a single synaptic time scale τ assuming an exponential synaptic response function.

However we point out that re-scaling of time allows to set τ = 1.

In general, the connectivity kernel K(x, y) fully depends on both sites x and y reflecting spatial heterogeneity. If the connectivity solely depends on the difference between x and y, i.e. K(x, y) = K(xy), then the neural field activity does not depend on specific spatial locations and hence is translational invariant. This case is called spatially homogeneity [11]. If the connectivity even depends on the distance between x and y only, i.e. K(x, y) = K(||x -y||), with ||x|| as some norm in Ω , then the neural field is spatially homogeneous and isotropic [13].

Several extensions of the Amari model (1) are possible, such as the consideration of finite axonal transmission speeds [14,15], constant feedback delays [16,17](see also section 3), heterogeneity [18,19] (see also section 4), spike-frequency adaption [20], statistical properties of single neurons [21], the combination of several brain areas [22], electromagnetic fields [23,[START_REF] Beim Graben | Neural Fields: Theory and Applications[END_REF] and many more [9].

Mathematically, Eq. ( 1) is an integro-differential equation. Spatially homogeneous (respectively isotropic) neural fields have been intensively studied in the literature due to their nice analytical properties [14,15,[START_REF] Coombes | [END_REF]26]. Moreover, these models may be transformed to derive partial differential wave equations [START_REF] Coombes | [END_REF]27,28] for certain classes of synaptic kernels.

Robinson model

As mentioned at the end of the previous paragraph, under certain conditions integrodifferential equations may be transformed to partial differential equations. In the 1990's James Wright and David Liley started developing partial differential equation models for neural population activity [29]. Their work has inspired other teams, e.g., Peter Robinson and colleagues, who developed a similar neural model which has been proven to be very successful. This type of neural field model [30][31][32] is based on a population-level model of a single thalamo-cortical module consisting of excitatory (E) and inhibitory (I) cortical population, thalamic relay neurons (S), and thalamic reticular neurons (R). The average soma membrane potential is modeled by

V a (t) = ∑ b=E,I,R,S h(t) ⊗ K a,b φ b (t -τ a,b ), a = E, I, S, R (2)
where ⊗ denotes temporal convolution and φ b is the pulse firing rate of the population b. The constants K a,b are the strengths of the connections from population of type b to population of type a. The delay term, τ a,b is zero for intra-cortical and intra-thalamic connections and non-zero for thalamocortical or corticothalamic connections [31].

The model assumes that only axons of excitatory cortical neurons are long enough to emit axonal propagating pulses. Moreover, φ E obeys the damped oscillator equation

Dφ E = S(V a ), (3)
with the operator

D D = 1 γ ∂ ∂t + 1 2 . (4)
In Eq. ( 2) h(t) denotes the mean synaptic response function

h(t) = αβ β -α e -αt -e -βt , (5)
where α and β are the synaptic decay and rise rate of synaptic response function, respectively.

Delayed nonlocal feedback between populations

In neural fields, one might include delays in several ways. The finite axonal transmission delay is proportional to the fraction of distance between two spatial locations and transmission speed and takes into account the finite propagation speed of action potentials along axonal branches [33], or in more general terms, it considers the finite-time interaction between two elements in a spatially extended system [28].

In addition, one could argue that delayed interactions happen on a single-neuron scale between neurons and it is more reasonable to treat these inter-neuron delays as a kind of effective delay [17,34]. This latter type of delay is constant. In addition to these two delay types, the nonlocal feedback delay takes into account the finite axonal transmission speed along axonal pathways between two brain areas. Since this axonal pathway has a finite defined length, the transmission delay is fixed and hence also constant [16]. Variations of all these delay types may be considered by distributed transmission speeds and/or distributed delays [15,35]. The subsequent paragraphs consider constant delays reflecting the finite transmission speed along axonal branches between brain areas.

The primary sensory area in weakly-electric fish

To examine the neural decoding in weakly electric fish, Doiron et al. [36] had performed an experimental in vivo stimulation study. A dipole was placed near the skin of a fish to stimulate only a part of the receptive field. A rather simple population model considers the primary sensory areas in the electro-sensory system of weakly electric fish [37], but similar configurations can also be found in parts of the vertebrate brain [38]. The model sketched in Fig. 2 [39] is made up of the ELL, a layer of pyramidal cells driven by the primary receptors that receive an external stimulus, and the higher area N p. These areas are spatially coupled via a delayed topographic feedback with connectivity kernels K en (x) and K ne (x) which reflect connections from the N p(n) to the ELL(e) and vice versa, respectively, see Fig. 2. The neurons in both populations have insignificant direct couplings and the coupling from the ELL to N p is excitatory and delayed in time by τ 1 , the coupling back to the ELL is inhibitory with delay τ 2 . Moreover, according to the experimental setup in [36], the model considers excitatory spatiotemporal stimuli I(x,t) to the ELL. The aim of the study is to learn more about the mechanism how to change the principal oscillation frequency by the properties of an external stimulus. The model considers spatiotemporal noise with a well-defined and adjustable spatial correlation length. The population model [39] describes two coupled neural fields whose activities are strongly related to experimentally observable local field potentials [START_REF] Nunez | Neocortical dynamics and human EEG rhythms[END_REF].

Then the theoretical power spectrum in the ELL reads

P(ν) = ∞ -∞ R(ν, l) C(l)dl (6)
where R(ν, l) is the spectral response function and C(l) is the scaled Fourier transform of the input correlation function. It turns out, that the power spectrum (6) does not depend on the spatial scale of the feedback loop σ f and the input correlation scale σ i independently, but just depends on their ratio, called η. This finding reflects the coupling of the spatial scale of the external input to the intrinsic spatial scale of the system. Figure 3 shows the resulting power spectra for two values of the spatial scale ratio η. We observe that a small ratio η = σ f /σ i ≪ 1 generates a spectral peak at about 20Hz, whereas the large ratio η ≫ 1 generates a power peak at about 0Hz. Hence retaining the topographic feedback but decreasing the input correlation function from large values of σ i = σ f /η (global noise) to small values of σ i (local noise) switches the spectral peak, similar as observed experimentally by Doiron et al. [36]. To understand this, Figure 4 shows the response function R(ν, l) and the integrand R(ν l ) C(l) in the definition of the power spectrum (6). The response function R and R C have a single maximum at about 20Hz for global noise (Fig. 4(a,b)), whereas R and R C have two local maxima at 20Hz and 0Hz (Fig. 4(c,d)). Since the peak of R C at 0Hz is broader than the peak about 20Hz and the power spectrum is the integral over R( l), cf. Eq. ( 6), the contribution of R C to the power at 0Hz exceeds the contribution at 20Hz yielding a strong peak at 0Hz. This result reveals that the spectral peak in Fig. 3 at 0Hz results from the selection of one mode out of two possible modes at frequencies of 0Hz and 20Hz, whereas the spectral peak at 20Hz is the only oscillation mode present in the system. The switch between these two configurations depends on the spatial correlation of the stimulus noise. These modes reflect activity subnetworks from which only one is engaged.

General anaesthesia

General anaesthesia is an important medical application in today's hospital surgery, but its underlying neural interactions is still a mystery. In the last decades, general anaesthesia has attracted theoretical researchers [START_REF] Steyn-Ross | [END_REF][42][43][START_REF] Ching | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Mccarthy | [END_REF]. Most theoretical studies aim to explain signal features of electroencephalographic data (EEG) observed during anaesthesia, such as the attenuation or enhancement of α-activity accompanied by a subsequent enhancement of δ -activity while increasing anaesthetic concentration [START_REF] Cimenser | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Murphy | [END_REF], cf. Fig. 5. The subsequent paragraphs show how neural field models Fig. 5 The power spectra measured in frontal EEG electrodes in the absence of anaesthesia (blue line) and during propofol anaesthesia (red line) in a group of subjects (a) and for a single subject (b). Taken from [30] by permission. may explain on the power enhancement and the frequency shift of maximum power while increasing the anaesthetic concentration.

To this end, we consider a derivative of the Robinson model introduced in section 2.2 and introduce a new sigmoid function derived from properties of type-I neurons [48] S(V a ) = F(V a , 0) -F(V a , γ),

with

F(V a , γ) = S max 2 1 + erf V a -θ -γσ 2 √ 2σ e -γ(V a -θ )+γ 2 σ 2 /2 , (8)
in which the parameter γ < ∞ takes into account the properties of type I-neurons, S max is the maximum population firing rate, θ is the mean firing threshold, and σ is related to the standard deviation of firing thresholds in the populations.

External input to the system can be considered as a non-specific input to relay neurons as

φ N = φ N + √ 2κξ (t), (9)
where φ N is its mean value and ξ (t) is Gaussian white noise of strength κ with zero mean.

It is well-known that general anaesthetics (GA) bind directly to sensitive target receptors. For instance, a large number of studies link the effect of many GAs to altered function of the GABA A receptors. Recent clinical findings have revealed several sites of elicited anaesthetic action of the anaesthetic propofol in the human brain, e.g., propofol suppresses field potentials in the rat thalamus and cortex [49].

In detail, the anaesthetic propofol increases the decay time constant of synaptic GABA A receptors, and hence increases the total charge transfer in these synapses but not that of excitatory synapses [50]. To integrate physiological observations into neural models such as the Robinson model detailed in section 2.2, the anaesthetic action on synaptic receptors is modelled by α → α/p with p ≥ 1 leading to a decrease of the decay rate constant α, or equivalently, an increase in decay time constant of GABA A receptors [51]. The factor p = 1 reflects absent anaesthetic action, i.e, the baseline condition. The model considers inhibitory synapses at excitatory neurons (factor p 1 ), and at inhibitory neurons (factor p 2 ) and at thalamic relay neurons (factor p 3 ).

To compute the power spectrum of the system, we consider the stationary state of Eq. ( 2), which obeys dV a (t)/dt = 0 where V a is taken from Eq. ( 2). Increasing the anaesthetic concentration, i.e. increasing the three factors p 1 , p 2 and p 3 changes the stationary states dependent on the relation of these three factors, cf. Fig. 6. In (a) the two lower stationary states collide to a single state whereas in (b) the two upper states collide. This difference indicates two fundamentally different mechanisms which may yield different power spectra.

The power spectrum characterises small fluctuations about this stationary state. Assuming that excitatory activity generates the EEG, and by virtue of the specific choice of external input to relay neurons, the power spectrum of the EEG depends just on one matrix component of the Greens function [51] by

P E (ω) = 2κ √ 2π G1,3 (ω) 2 . ( 10
)
in which G1,3 (ω) is a matrix element of the 4 × 4 matrix

G(ω) = 1 √ 2π     L D K 11 -K 1 -K 2 -K 3 e -iωτ 0 -K 4 L -K5 -K 6 0 -K 7 e -iωτ 0 L -K 8 -K 9 e -iωτ 0 -K 10 L     -1 (11)
with L = (1 + iω/α)(1 + iω/β ), D = (1 + iω/γ) 2 , and the constants K i , i = 1, 2, ..., 11 are proportional to the synaptic strengths and the nonlinear gains ∂ S/∂V computed at the stationary state of the system. Figure 7 shows the theoretical power spectrum P E in the baseline condition and after the administration of propofol for two different relations of the anaesthetic factors p 1 , p 2 and p 3 . At first we note that the spectra resemble well the spectrum obtained from experimental observations: increasing p 2 in a specific relation to the p 1 and p 3 yields increases in delta and theta power as well as more pronounced alpha oscillation with increased peak-frequency. The dynamical analysis of the model [30] reveals three resonances in the baseline condition, including an oscillatory resonance corresponding to the peak in the alpha-band and a pair of zero-frequency resonances. Increasing the anaesthetic concentration diminishes the damping rate of alpha resonances (and hence increases its magnitude) while its frequency increases. Hence increases in alpha power and its peak-frequency results from the approach of the system of an oscillatory instability [52]. Moreover, the two zero-frequency resonances collide and gradually increases in frequency leading to a magnitude increases in delta and theta power (solid red line in Fig. 7.) Summarizing, previous studies [30,51] have revealed the differential role of synaptic inhibition at GABAergic receptors located on the dendrites of different neural types. Moreover, synaptic potentiation within local cortical inhibitory neurons suffices to reproduce experimental observation in EEG during anesthesia [30,52].

Heterogeneous neural fields

Previous experimental studies on the structure of biological neural networks and their connections, e.g., as the exciting work of Hellwig [53], reveals a large rather simple mesoscopic structure at the spatial scale of several hundreds of micrometers with an underlying more complex structure at smaller scales. For instance, [53] reveals a rather homogeneuous Gaussian or exponential connectivity distribution of single neurons in layers II and III in rat visual cortex on a scale of ∼ 500µm with an overlying complex connectivity structure at much shorter scales reflecting heterogeneous connections. Hence, a homogeneous model serves as a first approximation at larger spatial scales, whereas heterogeneous connections need to be considered in more realistic models.

Several previous theoretical studies have investigated the effect of heterogeneous connections in neural fields, primarily in the context of pattern formation [54][55][56][57][58]. Most of the previous studies consider the heterogeneous fields as small heterogeneous perturbations of homogeneous states. The subsequent paragraphs show a new approach taking into account the heterogeneous structure not as a limit of the homogeneous case.

Starting from the Amari equation ( 1), we expand the integral into a power series [19] ∂V (x,t) ∂t +V (x,t)

=

Ω K 1 (x, y)V (y,t) dy + Ω Ω K 2 (x, y, z)V (y,t)V (z,t) dy dz ,
where the kernels K 1 and K 2 have to be reconstructed from experimental data.

Event-related potentials

In order to illustrate this, we combine a recently developed technique for nonlinear data analysis [18] with our new method for training heterogeneous neural fields [19] [19]. Figure 8 shows the grand average event-related potentials (ERP) U n (t), n = 1 . . . , N of N = 25 recording electrodes from a language processing experiment of subject-object ambiguities in German [59]. The first column displays the control condition with a canonical subject-verb-object order: "Die Rednerin hat den Berater beim Kongress gesucht" ("The speaker has sought the advisor at the congress"), while the second column (Fig. 8(e -h)) shows ERPs for a non-canonical, yet in German grammatical, object-verb-subject order: "Die Rednerin hat der Berater beim Kongress gesucht" ("The speaker has been sought by the advisor at the congress"). The panels (a) and (b) in the first row of Fig. 8 present the ERPs averaged over 14 subjects elicited by the disambiguating article "den" vs. "der" at time zero. Each trace shows the voltage of one of N = 25 EEG electrodes. Comparing panels (a) and (b) in Fig. 8, one recognises a difference between conditions at about 600ms after onset of the critical stimulus, known as the P600 ERP component.

Heteroclinic orbits

To analyze these spatiotemporal ERP patterns we apply our recently developed recurrence domain segmentation technique [18] whose results are shown in Fig. 8(c) and (d). This method computes the recurrence plots from the N-dimensional trajectories of each condition using a ball-size ε > 0 as parameter. Here, we optimized ε according to a reasonable heuristics for obtaining a relatively low number of recurrence domains with relatively long dwell times and good (visual) discrimination of conditions which worked quite well for ε = 1.9µV . Following [18], the recurrence plots are interpreted as rewriting grammars to replace large time indices by lower, recurrent, ones. These transformed time indices are plotted in Fig. 8(c) and (d) as different colors. The processing differences between conditions are clearly seen as different segments, indicating that ERP components, such as the P600, can be regarded as recurrence domains, or, in first approximation, as saddle nodes [60,61]. These saddles are estimated in Fig. 8(c) and (g) as follows: For each segment k of Fig. 8(c) and (d) we compute its center of gravity, i.e. its temporal average, to obtain a spatial EEG activity distribution on the scalp, namely a spatial voltage distribution U k (z), k = 1 . . . , P with number of segments P and where z is the spatial location on the scalp. The patterns U k (z) are activity patterns extrapolated in space on the basis of the average ERP components which are discrete in space. They are shown in Fig. 8(e) and (f). These distributions are assumed to originate from underlying neural patterns embed-ded in neural populations, e.g. sub-networks in neural populations. These underlying corresponding neural patterns may be called V k (x). Let us call V (x,t) the neural population activity generating the EEG signal. Then the Moore-Penrose pseudoinverses V k (x) + of V k (x) form a bi-orthogonal basis with the patterns V k (x), such that the projections

ξ k (t) = Ω V k (x) + V (x,t)dx
serve as order parameters in a decomposition

V (x,t) = ∑ k ξ k (t)V k (x) .
The ξ k (t) are plotted in Figs. 8(g) and (h) where each trace depicts the projection onto one segment from the recurrence domain analysis in Figs. 8(c) and (d).

The dynamics of the order parameters {ξ k } is of high temporal complexity, indicating that a large number of time scales is involved in the ERP language processing dynamics. However, in a first approximation we assume that the recurrence domains are in fact isolated saddles that are connected along stable heteroclinic sequences (SHS). This can be modeled through winnerless competition in a network of Lotka-Volterra populations [62,63]. This assumption leads to the order parameter dynamics shown in 8(g,h). From the growth rates σ k > 0, and interaction weights ρ k j > 0, ρ kk = 1, of the Lotka-Volterra dynamics and from the bi-orthogonal pattern systems V j (x),V k (x) + , one can construct the kernels of the neural field equation ( 12) as

K 1 (x, y) = ∑ k (σ k + 1)V + k (y)V k (x) K 2 (x, y, z) = ∑ k j ρ k j σ j V + k (y)V + j (z)V k (x) ,
We observe that the kernel K 1 (x, y) describes Hebbian synapses between sites y and x trained with pattern sequence V k [19]. Interestingly, this memory storage mechanism resembles well the storage of patterns in a Bidirectional Associative Memory (BAM) [64]. Moreover, to our best knowledge the three-point kernel K 2 (x, y, z) has not been studied yet in the context of neural fields. It further generalizes Hebbian learning to interactions between three sites x, y, z ∈ Ω . Interestingly, one can write

K 1 (x, y) = ∑ k (σ k + 1)K 1,k (x, y) K 2 (x, y, z) = ∑ j ∑ k ρ k j K 1,k (x, y) σ j V + j (z) .
with K 1,k (x, y) = V + k (y)V k (x). According to Eq. ( 13) the three-point kernel can be see as a linear superposition of two-point kernels weighted by the interaction weights ρ k j and the underlying stored patterns V + j (z) and its corresponding growth rates σ k .

The additional dependence of the spatial location z is reasonable in heterogeneous fields: the action at spatial location x does not only depend on the activity in spatial location y (as assumed in homogeneous systems), but may also depend on the path between x and y which in turn depends on the underlying memory structure. This argument motivates the presence of an additional spatial variable z, but does neither explain the terms in K 2 (x, y, z) nor its dependence on the specific choice of the heteroclinic dynamics (here a Lotka-Volterra system). Future work will attempt to elaborate more details on the presence of three-point kernels.

Conclusion

This chapter has presented some theoretical neural field studies describing the power spectra in homogeneous neural fields and the nonlinear dynamics in heterogeneous neural fields. In this context, we want to mention recent extensions of standard neural field models by incorporating single neuron properties [48,[65][66][67].

Figure 1

 1 sketches the experimental setup, the spike autocorrelation A(t) and the interspike interval (ISI) histograms of a typical pyramidal cell response to a local and global stimulus. The stimulus was temporal noise evoking spatially weakly correlated sensory receptor activity for local stimulation and strong spatial correlations in the receptor dynamics in case of the experimentally global stimulation. The global random stimulus evokes bursting in the neuron activity whereas the local stimulus evokes a single principal firing mode. This experiment raises the question how the spatial correlations in the input stimulus interact with those imposed by the physiological system.

Fig. 1

 1 Fig.1Experimental setup and firing statistics of stimulation experiment[36]. In (a), the electric skin stimulation is local (left panel) and induces a single main oscillation as seen in the spike autocorrelation function A(τ) (center panel) and the corresponding histogram (right panel). In (b) the stimulus was global inducing an additional oscillation mode. Taken from[36] by permission.

Fig. 2

 2 Fig. 2 Topography of the delayed feedback model. The plus and minus signs indicate the excitatory and inhibitory connections, respectively.

Fig. 3

 3 Fig. 3 Theoretical power spectrum computed for η = 40 (solid line) and η = 1/40 (dashed line).

Fig. 4

 4 Fig. 4 The response function R(ν, l) and the integrand of the power spectrum integral R(ν, l) C(l) for (a,b) σ i = 40σ f (global noise), (c,d) σ i = σ f /40 (local noise). Taken from [39] by permission.

Fig. 6

 6 Fig. 6 The stationary states and the nonlinear gain dS/dV computed at the lowest stationary state of pyramidal neurons V E subjected to the factor p 2 . (a), (b) p 1 = p 3 = 1 + 0.3(p 2 -1), (c), (d) p 1 = p 2 = p 3 . We observe three states in (a) and (c) for p 2 = 1 where (a) the two lower states collide and (c) the two upper states collide. The center branch (dashed red) is linearly unstable, whereas the other branches are linearly stable. For clarification, the lower branches are shown in the insets. Parameters are S max = 250Hz, V th = 15mV, γ = 0.08mV and σ = 10mV.

Fig. 7

 7 Fig. 7 Theoretical power spectrum in the baseline condition (p 1 = p 2 = p 3 = 1.0 encoded in blue) and in the anesthesia condition (red line). (a) p 1 = p 3 = 1 + 0.3(p 2 -1) and p 2 = 1.15 (red line) and (b) p 1 = p 2 = p 3 = 1.15 (red line).

Fig. 8

 8 Fig.8Experimental event-related potentials (ERP)[59] in two experimental conditions (left and right column) and their interpretation as heteroclinic orbits. (a,b) ERPs of N = 25 single EEG time series, (c,d) temporal sequence of extracted recurrence domains, (e,f) the spatial EEG activity distributions on the scalp (ERP components, seen from above, nose on top, back of head on bottom of the maps) averaged over the time windows of recurrence domains shown in panels (c,d), projections of the multivariate EEG signal (shown in panels (a,b)) on the corresponding ERP component maps (shown in panels (e,f), (g,h) sequence of heteroclinic saddles (SHS) modeled by the Lotka-Volterra model.

Acknowledgements

The authors thank Stefan Frisch and Heiner Drenhaus for conducting the ERP experiment. AH and MH acknowledge funding from the European Research Council for support under the European Union's Seventh Framework Programme (FP7/2007-2013) ERC grant agreement No.257253. PbG acknowledges financial support through a Heisenberg Fellowship of the German Research Foundation DFG (GR 3711/1-2).


Nathanaël Fijalkow

Hugo Gimbert

Florian Horn

Youssouf Oualhadj

Two Recursively Inseparable Problems for Probabilistic Automata ⋆

This paper introduces and investigates decision problems for numberless probabilistic automata, i.e. probabilistic automata where the support of each probabilistic transitions is specified, but the exact values of the probabilities are not. A numberless probabilistic automaton can be instantiated into a probabilistic automaton by specifying the exact values of the non-zero probabilistic transitions. We show that the two following properties of numberless probabilistic automata are recursively inseparable:

• all instances of the numberless automaton have value 1,

• no instance of the numberless automaton has value 1.

Introduction

In 1963 Rabin [START_REF] Michael | Probabilistic automata[END_REF] introduced the notion of probabilistic automata, which are finite automata able to randomise over transitions. A probabilistic automaton has a finite set of control states Q, and processes finite words; each transition consists in updating the control state according to a given probabilistic distribution determined by the current state and the input letter. This powerful model has been widely studied and has applications in many fields like software verification [START_REF] Chatterjee | Algorithms for omega-regular games of incomplete information[END_REF], image processing [START_REF] Culik | Digital images and formal languages[END_REF], computational biology [START_REF] Durbin | Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids[END_REF] and speech processing [START_REF] Mohri | Finite-state transducers in language and speech processing[END_REF].

Several algorithmic properties of probabilistic automata have been considered in the literature, sometimes leading to efficient algorithms. For instance, functional equivalence is decidable in polynomial time [START_REF] Schützenberger | On the definition of a family of automata[END_REF][START_REF] Tzeng | A polynomial-time algorithm for the equivalence of probabilistic automata[END_REF], and even faster with randomised algorithms, which led to applications in software verification [START_REF] Kiefer | Language equivalence for probabilistic automata[END_REF].

However, many natural decision problems are undecidable, and part of the literature on probabilistic automata is about intractability results. For example the emptiness, the threshold isolation and the value 1 problems are undecidable [START_REF] Paz | Introduction to probabilistic automata[END_REF][START_REF] Bertoni | Some recursive unsolvable problems relating to isolated cutpoints in probabilistic automata[END_REF][START_REF] Gimbert | Probabilistic automata on finite words: Decidable and undecidable problems[END_REF].

A striking result due to Condon and Lipton [START_REF] Condon | On the complexity of space bounded interactive proofs (extended abstract)[END_REF] states that, for every ǫ > 0, the two following problems are recursively inseparable: given a probabilistic automaton A, ⋆ The research leading to these results has received funding from the European Union's Seventh Framework Programme (FP7/2007-2013) under grant agreement n o 259454 (GALE) and from the French Agence Nationale de la Recherche projects EQINOCS (ANR-11-BS02-004) and STOCH-MC (ANR-13-BS02-0011-01).

• does A accept some word with probability greater than 1ε?

• does A accept every word with probability less than ε?

In the present paper we focus on numberless probabilistic automata, i.e. probabilistic automata whose non-zero probabilistic transitions are specified but the exact values of the probabilities are not. A numberless probabilistic automaton can be instantiated into a probabilistic automaton by specifying the exact values of the non-zero probabilistic transitions (see Section 2 for formal definitions).

The notion of numberless probabilistic automaton is motivated by the following example. Assume we are given a digital chip modelled as a finite state machine controlled by external inputs. The internal transition structure of the chip is known but the transitions themselves are not observable. We want to compute an initialisation input sequence that puts the chip in a particular initial state. In case some of the chip components have failure probabilities, this can be reformulated as a value 1 problem for the underlying probabilistic automaton: is there an input sequence whose acceptance probability is arbitrarily close to 1? Assume now that the failure probabilities are not fixed a priori but we can tune the quality of our components and choose the failure probabilities, for instance by investing in better components. Then we are dealing with a numberless probabilistic automaton and we would like to determine whether it can be instantiated into a probabilistic automaton with value 1, in other words we want to solve an existential value 1 problem for the numberless probabilistic automaton. If the failure probabilities are unknown then we are facing a second kind of problem called the universal value 1 problem: determine whether all instances of the automaton have value 1. We also consider variants where the freedom to choose transition probabilities is restricted to some intervals, that we call noisy value 1 problems.

One may think that relaxing the constraints on the exact transition probabilities makes things algorithmically much easier. However this is not the case, and we prove that the existential and universal value 1 problems are recursively inseparable: given a numberless probabilistic automaton C,

• do all instances of C have value 1?

• does no instance of C have value 1?

This result is actually a corollary of a generic construction which constitutes the technical core of the paper and has the following properties. For every simple probabilistic automaton A, we construct a numberless probabilistic automaton C such that the three following properties are equivalent: (i) A has value 1, (ii) one of the instances of C has value 1, (iii) all instances of C have value 1.

The definitions are given in Section 2. The main technical result appears in Section 3, where we give the generic construction whose properties are described above. In Section 4, we discuss the implications of our results, first for the noisy value 1 problems, and second for probabilistic Büchi automata [START_REF] Baier | Probabilistic ω-automata[END_REF].

Definitions

Let A be a finite alphabet. A (finite) word u is a (possibly empty) sequence of letters u = a 0 a 1 • • • a n-1 ; the set of finite words is denoted by A * .

A probability distribution over a finite set Q is a function δ : Q → Q ≥0 such that q∈Q δ(q) = 1; we denote by 1 3 • q + 2 3 • q ′ the distribution that picks q with probability 1 3 and q ′ with probability 2 3 , and by q the trivial distribution picking q with probability 1. The support of a distribution δ is the set of states picked with positive probability, i.e., supp(δ) = {q ∈ Q | δ(q) > 0}. Finally, the set of probability distributions over Q is D(Q).

Definition 1 (Probabilistic automaton). A probabilistic automaton (PA) is a tuple A = (Q, A, q 0 , ∆, F ), where Q is a finite set of states, A is the finite input alphabet, q 0 ∈ Q is the initial state, ∆ : Q × A → D(Q) is the probabilistic transition function, and F ⊆ Q is the set of accepting states.

For convenience, we also use P A (s u -→ t) to denote the probability of going from the state s to the state t reading u, P A (s u -→ S) to denote the probability of going from the state s to a state in S reading u, and P A (u), the acceptance probability of a word u ∈ A * , to denote P A (q 0 u -→ F ). We often consider the case of simple probabilistic automata, where the transition probabilities can only be 0, 1 /2, or 1.

Definition 2 (Value).

The value of a PA A, denoted val(A), is the supremum acceptance probability over all input words val(A) = sup

u∈A * P A (u) .
The value 1 problem asks, given a (simple) PA A as input, whether val(A) = 1.

Theorem 1 ([8]

). The value 1 problem is undecidable for simple PA.

Definition 3 (Numberless probabilistic automaton). A numberless probabilistic automaton (NPA) is a tuple

A = (Q, A, q 0 , T, F ), where Q is a finite set of states, A is the finite input alphabet, q 0 ∈ Q is the initial state, T ⊆ Q × A× Q is the numberless transition function, and F ⊆ Q is the set of accepting states.
The numberless transition function T is an abstraction of probabilistic transition functions. We say that ∆ is consistent with T if for all letters a and states s and t, ∆(s, a, t) > 0 if, and only if (s, a, t) ∈ T .

A numberless probabilistic automaton is an equivalence class of probabilistic automata, which share the same set of states, input alphabet, initial and accepting states, and whose transition functions have the same support.

A NPA A = (Q, A, q 0 , T, F ) together with a probabilistic transition function

∆ consistent with T defines a PA A[∆] = (Q, A, q 0 , ∆, F ). Conversely, a PA A = (Q, A, q 0 , ∆, F ) induces an underlying NPA [A] = (Q, A, q 0 , T, F ), where T ⊆ Q × A × Q is defined by (q, a, p) ∈ T if ∆(q, a)(p) > 0.
We consider two decision problems for NPA:

-The existential value 1 problem: given a NPA A, determine whether there exists ∆ such that val(A[∆]) = 1. -The universal value 1 problem: given a NPA A, determine whether for all ∆, we have val(A[∆]) = 1.

Proposition 1. There exists a NPA such that:

-there exists ∆ such that val(A[∆]) = 1, -there exists ∆ ′ such that val(A[∆ ′ ]) < 1. C1 C2 L1 L2 R1 R2 a f i, 1
3 Recursive inseparability for numberless value 1 problems

In this section, we prove the following theorem:

Theorem 2. The two following problems for numberless probabilistic automata are recursively inseparable:

all instances have value 1, -no instance has value 1.

Recall that two decision problems A and B are recursively inseparable if their languages L A and L B of accepted inputs are disjoint and there exists no recursive language L such that L A ⊆ L and L ∩ L B = ∅.

Note that it implies that both A and B are undecidable. Equivalently, this means there exists no terminating algorithm which has the following behaviour on input x:

if x ∈ L A , then the algorithm answers "YES" -if x ∈ L B , then the algorithm answers "NO".

On an input that belongs neither to L A nor to L B , the algorithm's answer can be either "YES" or "NO".

Overall construction

Lemma 1. There exists an effective construction which takes as input a simple PA A and constructs a NPA C such that

val(A) = 1 ⇐⇒ ∀∆, val(C[∆]) = 1 ⇐⇒ ∃∆, val(C[∆]) = 1 .
We first explain how Lemma 1 implies Theorem 2. Assume towards contradiction that the problems "all instances have value 1" and "no instance has value 1" are recursively separable. Then there exists an algorithm A taking a NPA as input and such that: if all instances have value 1, then it answers "YES", and if no instance has value 1, then it answers "NO". We show using Lemma 1 that this would imply that the value 1 problem is decidable for simple PA, contradicting Theorem 1. Indeed, let A be a simple PA, applying the construction yields a NPA C such that

val(A) = 1 ⇐⇒ ∀∆, val(C[∆]) = 1 ⇐⇒ ∃∆, val(C[∆]) = 1 .
In particular, either all instances of C have value 1, or no instance of C has value 1. Hence, if it answers "YES" then val(A) = 1 and it if answers "NO" then val(A) < 1, allowing to decide whether A has value 1 or not. This concludes the proof of Theorem 2 assuming Lemma 1.

The construction follows two steps.

The first step is to build from A a family of PA's B λ whose transitions are all of the form q a -→ (λ • r, (1λ) • s) in such a way that B λ has the same value as A for any value of λ. Note that, while all the B λ 's belong to the same NPA, they are not, in general, a NPA: for example, if A were the simple version of the automaton of Figure 1, the B λ 's would be the instances where x = y = λ, while the underlying NPA would also include the cases where x = y.

The second step is to build from the B λ 's a NPA C such that, for each probabilistic transition function ∆, there is a λ such that C[∆] has value 1 if, and only if B λ has value 1.

It follows that:

∃∆, val(C[∆]) = 1 =⇒ ∃λ, val(B λ ) = 1 =⇒ val(A) = 1 =⇒ ∀λ, val(B λ ) = 1 =⇒ ∀∆, val(C[∆]) = 1 .

The fair coin construction

Let A = (Q, A, q 0 , ∆, F ) be a simple PA over the alphabet A. We construct a family of PAs (B λ ) λ∈]0,1[ over the alphabet B = A ∪ {♯}, whose transitions have probabilities 0, λ, 1λ or 1, as follows.

The automaton B λ is a copy of A where each transition of A is replaced by the gadget illustrated in Figure 2 (for simplicity, we assume that all the transitions of A are probabilistic). The initial and final states are the same in A and B λ .

The left hand side shows part of automaton A: a probabilistic transition from q reading a, leading to r or s each with probability half. The right hand side shows how this behaviour is simulated by B λ : the letter a leads to an intermediate state q a , from which we can read a new letter ♯. Each time a pair of ♯'s is read, the automaton B λ goes to r with probability λ•(1-λ), goes to s with probability (1λ) • λ, and stays in q a with probability λ 2 + (1λ) 2 . Reading a letter other than ♯ while the automata is still in one of the new states leads to the new sink state ⊥, which is not accepting. Thus, the probability of going to r is equal to the probability of going to s, and we can make the probability of a simulation failure as low as we want by increasing the number of ♯'s between two successive "real" letters.

Let u be a word of A * . We denote by [u] k the word of B * where each letter a ∈ A of u is replaced by a • ♯ 2k . Conversely, if w is a word of B * , we denote by w the word obtained from w by removing all occurrences of the letter ♯.

Intuitively, a run of A on the word u is simulated by a run of B λ on the word [u] k . Whenever there is a transition in A, B λ makes k attempts to simulate it through the gadget of Figure 2, and each attempt succeeds with probability (1 -2λ • (1λ)), so each transition fails with probability:

A λ,k = 1 -(1 -2λ • (1 -λ)) k .
Proposition 2. The probabilistic automaton B λ satisfies: 1. For all q, r ∈ Q, a ∈ A, and k ∈ N, P B λ (q

q r s a, 1 2 a, 1 2 q qa qa, L qa, R r s ⊥ ♯ a ♯, λ ♯, 1 -λ ♯, 1 -λ ♯, λ ♯ ♯ A A A A ∪ {♯}
[a] k --→ r) = A λ,k • P A (q a -→ r) ,
2. For all q, r ∈ Q, u ∈ A * , and k ∈ N, P B λ (q

[u] k --→ r) = A |u| λ,k • P A (q u -→ r) .
3. For all q, r ∈ Q, w ∈ (A ∪ {♯}) * , and k ∈ N, P B λ (q w -→ r) ≤ P A (q w -→ r) .

It follows from Proposition 2 that for any λ, the value of B λ is equal to the value of A.

The simulation construction

All the B λ 's induce the same NPA, that we denote by B. The problem is that there are many other instances of B, whose values may be higher than the value of A (recall the example of Figure 1, where the B λ have value 1 /2, while there are instances of B with value 1). In this subsection, we construct a NPA C (over an extended alphabet C) whose instances simulate all the B λ 's, but only them.

The idea is that the new NPA should only have one probabilistic transition. An instance of this transition translates to a value for λ. Figure 3 describes a first attempt at this (notice that our convention is that an non-drawn transition means a loop, rather than a transition to a sink state). apply(b, q) Fig. 3. Naive fusion of the probabilistic transitions.

In this automaton, that we call B ′ , there are two copies of the set of states, and a single shared probabilistic transition q R $ -→ (λ•s 0 , (1-λ)•s 1 ) between them. In order to make all the probabilistic transitions of B happen in this center area, we use new letters to detect where the runs come from before the probabilistic transition, and where it should go afterwards.

For each pair of a letter b in B and a state q in Q, we introduce two new letters check(b, q) and apply(b, q). The letter check(b, q) loops over each state except the left copy of q, which goes to q R . The letter apply(b, q) loops over each state except s 0 , from where it goes to the λ-valued successor of (q, b) in B λ , and s 1 , from where it goes to the (1λ)-valued successor of (q, b) in B λ . The new letter next transition sends the run back to the left part once each possible state has been tested.

Thus, if we define the morphism by its action on letters:

b = check(b, q 0 )•$•apply(b, q 0 ) • • • check(b, q n-1 )•$•apply(b, q n-1 )•next transition ,
where the q i 's are the states of B λ , we get for any word u on the alphabet B,

P B λ (u) = P B ′ λ ( u).
The problem with this automaton is that one can "cheat", either by not testing an unwelcome state, or by changing state and letter between a check and the subsequent apply. In order to avoid this kind of behaviour, we change the automaton so that it is necessary to win an arbitrarily large number of successive times in order to approach the value 1, and we can test whether the word is fair after the first successful attempt. A side effect is that the simulation only works for the value 1: for other values, it might be better to cheat. The resulting automaton is described in Figure 4.

The structure of the automaton of Figure 3 is still there, but it has been augmented with an extra layer of scrutiny: each time we use the probabilistic transition, there is now a positive probability (1θ) to go to a new wait state. There is also a new letter next word which has the following effect:

if the run is in an accepting state, it goes to the initial state of the fairness checker D; -if the run is in a non-accepting state, it goes to the non-accepting sink state ⊥ of D; -if the run is in the wait state, it goes back to the initial state.

The fairness checker D is a deterministic automaton which accepts the language

{ u • next word | u ∈ B * } * . Its final state is the only final state in all of C.
Intuitively, a run can still cheat before the first next word letter, but the benefits of doing so are limited: the probability that C[λ, θ] accepts a word at that point is at most θ (except if the empty word is accepting, but that case is trivial). After that point, cheating is risky: if the run already reached D, a false move will send it to the sink.

A more formal proof follows. A simple inspection of the construction of C yields Proposition 3:

D p q0 q q ′ ∈ F q ′′ / ∈ F wait s qR s0 s1 p p0 p1 ⊥ q D 0 next transition check(b, q) next word next word $, λ • θ $, (1 -λ) • θ $, (1 -θ)
apply(b, q)

apply(b, q)

apply(b, q) Fig. 4. The Numberless Probabilistic Automaton C.

Proposition 3. Let u be a word of B * of length k. We have:

P C[λ,θ] ( u • next word) = θ k • P B λ (u) P C[λ,θ] (( u • next word) ℓ ) = (1 -(1 -θ k ) ℓ ) • P B λ (u) .
It follows from Proposition 3 that the value of C[λ, θ] is at least the value of B λ . Proposition 4. Let u be a word of (C \ {next word}) + . We have:

P C[λ,θ] (u) ≤ θ .
Proposition 5 formalises the fact that there is no point in cheating after the first next word letter: Proposition 5. Let u 1 , . . . , u k be k words of (C \ {next word}) * and w be the word

u 1 • next word • • • u k • next word. Then, for any 1 ≤ i ≤ k, if u i / ∈ B * , we have: P C[λ,θ] (w) ≤ P C[λ,θ] (u i • next word • • • u k • next word) .
Proof. After reading u 1 • • • u i-1 •next word, a run must be in one of the following three states: q 0 , q D 0 , and ⊥. As u i / ∈ B * , reading it from q D 0 will lead to ⊥. Thus,

P C[λ,θ] (w) = P C[λ,θ] (q 0 u1•••ui-1•next word ------------→ q 0 ) • P C[λ,θ] (u i • • • u k • next
, val(C[∆]) = 1 =⇒ ∃λ, val(B λ ) = 1 =⇒ val(A) = 1 =⇒ ∀λ, val(B λ ) = 1 =⇒ ∀∆, val(C[∆]) = 1 .
Theorem 2 follows.

Consequences

In this section, we show several consequences of the recursive inseparability results from Theorem 2 and of the construction from Lemma 1. The first is a series of undecidability results for variants of the value 1 problem. The second is about probabilistic Büchi automata with probable semantics, as introduced in [1].

The noisy value 1 problems

Observe that Theorem 2 implies the following corollary:

Corollary 1. Both the universal and the existential value 1 problems are undecidable.

We can go further. Note that the universal and existential value 1 problems quantify over all possible probabilistic transition functions. Here we define two more realistic problems for NPA, where the quantification is restricted to probabilistic transition functions that are ε-close to a given probabilistic transition function:

-The noisy existential value 1 problem: given a NPA A, a probabilistic transition function ∆ and ε > 0, determine whether there exists ∆ ′ such that |∆ ′ -∆| ≤ ε and val(A[∆ ′ ]) = 1.

-The noisy universal value 1 problem: given a NPA A, a probabilistic transition function ∆ and ε > 0, determine whether for all ∆ ′ such that |∆ ′ -∆| ≤ ε, we have val(A[∆ ′ ]) = 1.

It follows from Lemma 1 that both problems are undecidable:

Corollary 2. Both the noisy universal and the noisy existential value 1 problems are undecidable.

Indeed, we argue that the construction from Lemma 1 implies a reduction from either of these problems to the value 1 problem for simple PA, hence the undecidability. Let A be a simple PA, the construction yields a NPA C such that: This completes the reduction.

Probabilistic Büchi automata with probable semantics

We consider PA over infinite words, as introduced in [START_REF] Baier | Probabilistic ω-automata[END_REF]. A probabilistic Büchi automaton (PBA) A can be equipped with the so-called probable semantics, defining the language (over infinite words):

L >0 (A) = {w ∈ A ω | P A (w) > 0} .

It was observed in [START_REF] Baier | Probabilistic ω-automata[END_REF] that the value 1 problem for PA (over finite words) easily reduces to the emptiness problem for PBA with probable semantics (over infinite words).

Informally, from a PA A, construct a PBA A ′ by adding a transition from every final state to the initial state labelled with a new letter ♯. (From a nonfinal state, the letter ♯ leads to a rejecting sink.) As explained in [START_REF] Baier | Probabilistic ω-automata[END_REF], this simple construction ensures that A has value 1 if and only if A ′ is non-empty, equipped with the probable semantics.

This simple reduction, together with Theorem 2, implies the following corollary:

Corollary 3. The two following problems for numberless PBA with probable semantics are recursively inseparable:

all instances have a non-empty language, -no instance has a non-empty language.

Fig. 2 .

 2 Fig. 2. The fair coin gadget.

  val(A) = 1 ⇐⇒ ∀∆, val(C[∆]) = 1 ⇐⇒ ∃∆, val(C[∆]) = 1 .It follows that for any probabilistic transition function ∆ and any ε > 0, we have:val(A) = 1 ⇐⇒ (∀∆ ′ , |∆ ′ -∆| ≤ ε =⇒ val(C[∆ ′ ]) = 1) ⇐⇒ (∃∆ ′ , |∆ ′ -∆| ≤ ε ∧ val(C[∆ ′ ]) = 1) .

  word) , and Proposition 5 follows. Finally, Proposition 6 shows that C[λ, θ] cannot have value 1 if B λ does not. Proof. Let us write w = u 1 • next word • • • u k • next word with u 1 , . . . , u k ∈ (C \ {next word}) * . By Proposition 4, k > 1, and by Proposition 5 we can assume that u 2 , . . . , u k belong to B * . Let v 2 , . . . , v k be the words of B * such thatu i = v i .The C[λ, θ]-value of w can be seen as a weighted average of 1 (the initial cheat, with a weight of θ) and the B λ values of the v i 's (the weight of P B λ (v i ) is the probability that the run enters D while reading u i ). It follows that at least one of the v i 's has a B λ -value greater than the C[λ, θ]-value w Thus, for each λ and θ, the value of C[λ, θ] is 1 if and only if the value of B λ is 1. As all the B λ 's have the same value, which is equal to the value of A, we get: ∃∆

	Proposition 6. For all words w ∈ C * such that P C[λ,θ] (w) > θ, there exists a word v ∈ B * such that
	P B λ (u) ≥	P C[λ,θ] (w) -θ 1 -θ	.

a, 1 -x f a, x i, a, f i, 1 2 a, 1 -y f a, y i, a, f Fig. 1. This NPA has value 1 if and only if x > yIn this automaton, adapted from[START_REF] Gimbert | Probabilistic automata on finite words: Decidable and undecidable problems[END_REF][START_REF] Fijalkow | Deciding the value 1 problem for probabilistic leaktight automata[END_REF], the shortest word that can be accepted is i • f , as i goes from C 1 to L 1 , and f goes from L 1 to L 2 . However, there are as much chances to go to R 2 , so the value of i• f is 1 /2.If x is strictly less than y, one can tip the scales to the left by adding a's between the i and f : each time, the run will have more chances to stay left than to stay right. After reading i • a n • f , the probability of reaching L 2 is equal to x n , while the probability of reaching R 2 is only y n . There is also a very high chance that the run went back to C 1 , but from there we can simply repeat our word an arbitrary number of times.Let x, y, and ε be three real numbers such that 0 ≤ y < x ≤ 1, and 0 < ε ≤ 1. There is an integer n such that x n /(x n +y n ) is greater than 1ε /2, and an integer m such that (1x ny n ) m is less than ε /2. The word (i • a n • f ) m is accepted with probability greater than 1ε.On the other hand, if x ≤ y, there is no word with value higher than 1 /2.

Acknowledgments

We would like to thank the referees for their helpful comments.


Aurélien Schutz
email: aurelien.schutz@ims-bordeaux.fr

Lionel Bombrun
email: lionel.bombrun@ims-bordeaux.fr

Yannick Berthoumieu
email: yannick.berthoumieu@ims-bordeaux.fr

Intrinsic prior for Bayesian classification of texture images

Keywords: Information geometry, intrinsic prior, texture classification

This paper introduces an intrinsic prior distribution for supervised classification of texture images. First, we introduce the intrinsic prior distribution as the normal law on a Riemannian manifold. Next, based on this definition, we derive the estimation and classification schemes. Finally, we propose an application for the classification of texture images. Experiments on the VisTex texture database are conducted and demonstrate the interest of the proposed intrinsic classification algorithm.

INTRODUCTION

In the framework of natural texture image processing, parametric stochastic models have been well studied. When parametric stochastic modeling is associated with scale-space representation, very effective algorithms can be designed, obtaining successful results for a large class of natural textures. Providing an unified view of model estimation, classification and synthesis, multiscale scheme enables us to prototype stochastic models which have a non-Gaussian distribution.

However, based on scale-space decomposition and stochastic modeling, none of the existing works have addressed the problem of the intra-class diversity characterization to increase the classification performance. In this paper, we employ the Bayesian framework to take into account the intra-class diversity. From the concept of intrinsic prior, we develop a new supervised parametric classification algorithm.

Parametric classification is closely linked with the theory of statistical manifolds, which aims at providing a Riemannian structure to the parameters space of probability density functions (pdf). Many works have proposed extrinsic tools by embedding the manifold in an Euclidean space [START_REF] Han | Isometric Embedding of Riemannian Manifolds in Euclidean Spaces[END_REF]. Nevertheless, to compare two observations in a parameter space, the geometry of the Riemannian manifold should be considered which is the main objective of information geometry theory [START_REF] Amari | Translations of Mathematical monographs[END_REF]. In this context, the notion of intrinsic tools is the key point. For example, in an estimation problem, the estimate should be invariant under any reparametrization of the parameter space. This notion of intrinsic has hence been considered in a wide range of disciplines such as in the definition of intrinsic loss function [START_REF] Robert | Intrinsic loss functions[END_REF], intrinsic discrepancy in the context of Bayesian estimation [START_REF] Bernardo | Intrinsic Estimation, Bayesian statistics[END_REF] and detection [START_REF] Bernardo | Integrated Objective Bayesian Estimation and Hypothesis Testing, Bayesian statistics[END_REF], and also in intrinsic version of Cramér-Rao bound [START_REF] Smith | Covariance, subspace, and intrinsic Cramér-Rao bounds[END_REF]. In this paper, our contributions are threefold. First, in the Bayesian framework, we introduce the concept of intrinsic prior distribution as the normal law on a Riemannian manifold when the Jeffrey divergence is considered and derive an intrinsic estimation and classification scheme. Second, based on the intrinsic prior, we show that the optimal decision characterizing the classification procedure leads to a decision directly on the parameter space (Riemannian manifold). Third, we propose to validate the proposed methodology in a texture based image retrieval experiments.

The paper is structured as follows. Section 2 gives a stateof-the-art about intrinsic distribution and introduces the proposed intrinsic prior. Based on this definition, Section 3 derives the estimation and classification scheme. Section 4 introduces an application for the classification of texture images, and some experiments results are displayed on the Vis-Tex database. Conclusions and future works are finally reported in Section 5.

INTRINSIC PRIOR DISTRIBUTION

Let χ = (χ 1 , . . . , χ K ) be K independent sets of N i independent and identically distributed random variables (vectors) x according to a parametric model p(x|θ). Let ( θ1 , . . . , θK ) ∈ Θ be the K maximum likelihood estimates computed on these sets (χ 1 , . . . , χ K ). This collection of K parametric vectors can be described by its first and second order characteristics: the most central element (i.e. the centroid θ) and the standard deviation σ around this central point. Based on the definition of the entropic prior and the normal law on a Riemannian manifold, we introduce the notion of intrinsic prior distribution p(θ| θ, σ).

Entropic prior

Given a parametric model p(•|θ), the entropic prior on θ is given by [START_REF] Rodríguez | Entropic[END_REF][START_REF] Snoussi | Information Geometry and Prior Selection[END_REF]:

p (θ|θ 0 , α) ∝ 1 |G(θ)| -1 2 exp {-αI(p(•|θ 0 ), p(•|θ))} (1)
where α is a positive scalar parameter, |G(θ)| is the determinant of the Fisher information matrix computed at point θ and I(p(•|θ 0 ), p(•|θ)) is a divergence between the probability measures p(•|θ 0 ) and p(•|θ).

As observed, this entropic prior is invariant under any change of coordinate in the parameter space. This prior is hence intrinsic and assigns to θ a probability which decreases exponentially with the divergence. In (1), the parameter α controls the sensitivity to changes in distance. Note that when α tends toward 0, the entropic prior reduces to the non-informative Jeffrey prior.

Normal law on a manifold

In [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF], Pennec introduces the concept of normal law on a Riemannian manifold M knowing the mean value θ and covariance matrix Σ as:

p(θ| θ, Γ) = k exp - - → θθ T Γ - → θθ 2 ( 2
)
where k is the normalizing constant and Γ the concentration matrix linked to the covariance matrix by:

Σ = k M - → θθ - → θθ T exp - - → θθ T Γ - → θθ 2 dM(θ) (3)
The model defined in ( 2) approximates the normal model by the usual Gaussian distribution in the tangent space TθM at the mean value θ. In this definition, the projection from the Riemannian manifold to the tangent space is given by the exponential map. Here, the geodesic distance (GD) induced by the Riemannian metric, derived from the Fisher information matrix is considered to compute the proximity between two

observations, i.e. GD(p(•| θ), p(•|θ)) = - → θθ .
This model has notably been considered in [START_REF] Lenglet | Statistics on the manifold of multivariate Normal distributions: Theory and application to diffusion tensor MRI processing[END_REF] for the segmentation of magnetic resonance images.

Intrinsic prior

Inspired from the definitions of the entropic prior (1) and the normal law on a Riemannian manifold (2), we introduce the definition of the intrinsic prior as the normal law on a manifold when the Jeffrey divergence is considered instead of the Riemannian metric as:

p θ| θ, σ = 1 (2π) d/2 σ d |G(θ)| -1 2 exp - 1 2σ 2 J(p(•|θ), p(•| θ)) , ( 4
)
where d is the dimension of the parameter space Θ and

J(p(•|θ), p(•| θ))
is the Jeffrey divergence computed between the parametric models p(•|θ) and p(•| θ).

INTRINSIC ESTIMATION AND CLASSIFICATION 3.1. Intrinsic estimation

By considering the notation introduced in Section 2, the maximum likelihood estimators of the centroid θ and standard deviation σ are obtained as solution of:

λ = arg max λ p(χ|λ), (5)
where λ = { θ, σ} is the set of hyperparameters, and

p(χ|λ) = K i=1 p(χ i |λ) = K i=1 Θ p(χ i |θ i ) p(θ i |λ) dM(θ i ), (6)
where dM

(θ i ) = |G(θ i )| 1
2 dθ i is the volume element on a Riemannian manifold. In [START_REF] Smith | Covariance, subspace, and intrinsic Cramér-Rao bounds[END_REF], the main difficulty relies on the computation of the integral. Many works have been dedicated to this problem including Monte Carlo integration, variational based approaches or Laplace approximation [START_REF] Kass | Approximate bayesian inference in conditionally independent hierarchical models (parametric empirical bayes models)[END_REF][START_REF] Miyata | Fully exponential laplace approximations using asymptotic modes[END_REF]. This latter has been successfully validated by many authors in applied mathematics (see [START_REF] Johnson | An asymptotic expansion for posterior distributions[END_REF] for instance). Next, after developping p(χ i |θ i ) as Ni j=1 p(x j |θ i ) and introducing the Laplace approximation and the expression of the proposed intrinsic prior (4) in [START_REF] Smith | Covariance, subspace, and intrinsic Cramér-Rao bounds[END_REF], one can derive the maximum likelihood estimates of the centroid θ and standard deviation σ. After some cumbersome computations and removing the terms independent of θ and σ, it yields:

θ = arg min θ 1 K K i=1 J(p(•| θi ), p(•| θ)) (7) σ2 = 1 dK K i=1 J(p(•| θi ), p(•| θ)) (8)
First, it can be noticed that the maximum likelihood estimators of θ and σ do not depend on the set χ. They are directly expressed as a function of the maximum likelihood estimators θi . Note also that, equations [START_REF] Rodríguez | Entropic[END_REF] and ( 8) coincide with the maximum likelihood estimators of p(θ|λ) where θ = { θ1 , . . . , θK }.

To estimate θ from (7), a stochastic gradient descent algorithm can be considered [START_REF] Amari | Translations of Mathematical monographs[END_REF][START_REF] Amari | Natural gradient works efficiently in learning[END_REF][START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]. θ is obtained as the fixed point solution of

θ t+1 = θ t -η t C(θ t )∇l(θ t ) (9)
where η t is the step size which may depend on iteration t, C(θ t ) is a positive definite matrix and ∇l(θ t ) is the gradient of the cost function defined in [START_REF] Rodríguez | Entropic[END_REF]. When C(θ t ) is equal to G -1 (θ t ), the inverse of the Fisher information matrix, [START_REF] Pennec | Intrinsic statistics on riemannian manifolds: Basic tools for geometric measurements[END_REF] corresponds to the natural gradient descent. This latter is intrinsic since it does not depend on the chosen parametrization θ of the pdf [START_REF] Amari | Translations of Mathematical monographs[END_REF][START_REF] Amari | Natural gradient works efficiently in learning[END_REF][START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF].

Classification

Let χ t be a set of N t independent and identically distributed random vectors x. Let θt be the maximum likelihood estimate computed on χ t . Let λ1 , . . . , λC be a collection of C hyperparameters, corresponding to C classes, estimated according to [START_REF] Rodríguez | Entropic[END_REF] and [START_REF] Snoussi | Information Geometry and Prior Selection[END_REF]. The sample χ t is classified to the class

Equations ( 12) and ( 11) can respectively be interpreted as linear (resp. quadratic) discriminant analysis on a Riemannian manifold. Note that those decision rules are decisions on the parameter space Θ (not on the original space χ), hence reducing the computational complexity.

In the next section, we propose an application of this intrinsic classification scheme for the recognition of texture images.

INTRINSIC TEXTURE CLASSIFICATION 4.1. Context

Many works in texture image recognition have shown that the wavelet representation is a well-adapted domain to characterize the texture, yielding to a multiscale analysis scheme which consists in modeling each wavelet subband. Let I be a texture image. Let N o and N s be respectively the number of orientation and scale of a multi-scale decomposition. I is hence decomposed into N o × N s sub-bands. Let us consider the parametric vector θ s,o of the pdf associated to each sub-band. The collection T I of those parametric vectors will represent the texture image I. 7) and [START_REF] Snoussi | Information Geometry and Prior Selection[END_REF]. As observed in Fig. 1, the natural intra-class diversity of texture images is captured by the proposed normal law on a Riemannian manifold (intrinsic prior). In this paper, we propose an application to texture image recognition based on a multivariate modeling context. The wavelet coefficients located around the neighborhood p × q of the current spatial position are clustered in the random vector x. The realizations of vector x characterize the spatial dependency of wavelet coefficients. Spherically Invariant Random Vectors are a class of stochastic models which have shown promising results for modeling the spatial dependency of wavelet coefficients [START_REF] Lasmar | Multivariate statistical modeling for texture analysis using wavelet transforms[END_REF]. Hence, we propose to model those observations based on the SIRV representation. Let x be a pq-dimensional vector following a SIRV distribution, it yields that x admits the stochastic representation x = √ τ g where τ is a scalar random variable called multiplier (τ ∈ R + ) and g a real Gaussian vector with zero mean and covariance matrix Σ = E{gg T }. By exploiting the independence of the processes τ and g and by working on the joint vector y = (τ, g), the Jeffrey divergence of the joint model can be expressed as the sum of the Jeffrey divergence for the multivariate Gaussian process and the Jeffrey divergence for the multiplier part. Note that both terms admit a closed-form expression recalled in [START_REF] Schutz | Centroidbased texture classification using the SIRV representation[END_REF]. It yields that the centroid for a SIRV model y is composed by two centroids: one for the Gaussian part and one for the multiplier part. For more information dealing with the implementation of those centroids estimators, the interested reader is referred to [START_REF] Schutz | Centroidbased texture classification using the SIRV representation[END_REF].

Texture classification

Let T t be a test image. According to the classification rule presented in Section 3.2, this image is labeled to the class ĉ, corresponding to the class maximizing the likelihood p(T t |λ ĉ). Since the subbands of the wavelet decomposition are independent, one can consider the chain rule principle to obtain the pdf of p(T t ) as the product of the pdf computed for each subbands. After some computations, the decision rule is simply the sum of the decision rules [START_REF] Kass | Approximate bayesian inference in conditionally independent hierarchical models (parametric empirical bayes models)[END_REF]

(14)

When considering the multivariate SIRV model with Weibull distributed multiplier, the dimension d of the parameter space Θ is equal to (pq+1)pq 2 + 1 since one covariance matrix of dimension pq × pq and one shape parameter for the multiplier τ should be estimated. When a univariate model is consider to represent the wavelet coefficients such as the 2-parameters generalized Gaussian distribution (GGD), d is equal to 2.

Similarly, when the homoscedasticity assumption holds, the decision rule reduces to a nearest neighbor classifier according to the Jeffrey divergence.

Note that to the best of our knowledge, even if this last decision rule has been previously proposed in [START_REF] Choy | Supervised texture classification using characteristic generalized Gaussian density[END_REF][START_REF] Schutz | Kcentroids based supervised classification of texture images: handling the intra-class diversity[END_REF][START_REF] Shabbir | Multivariate texture discrimination based on geodesics to class centroids on a generalized Gaussian manifold[END_REF], no previous works had been fully formalized in the Bayesian framework.

Results and discussion

To evaluate the performance of the proposed supervised classification algorithm, the database is split into a training database and a disjoint testing database. From a practical point of view, N T r training samples are randomly selected for each texture class, the remaining samples are taken as testing samples. In the following, 100 Monte Carlo runs are used to evaluate the performance of the proposed classifiers. Performances are evaluated in terms of kappa index. The kappa index refers to the proportion of consistent classifications observed beyond that expected by chance alone [START_REF] Cohen | A coefficient of agreement for nominal scales[END_REF][START_REF] Gomez | Determining the accuracy in image supervised classification problems[END_REF]. This experiment is carried out on the MIT Vision texture (VisTex) [START_REF]MIT Vision and Modeling Group. Vision Texture[END_REF]. This database is composed of 40 classes and 64 images per class of size 64 × 64 pixels. All texture images are normalized in intensity to have zero mean and unit standard deviation. This normalization gives invariance to affine transformations in the illumination intensity. Here, the stationary wavelet decomposition with 2 scales and Daubechies' filter db4 have been used and a 3 × 3 neighborhood has been considered to model the spatial dependency of the wavelet coefficients.

Fig. 2 draws the evolution of the average kappa index as a function of the number of training samples on the VisTex database. Results for both univariate (GGD) and multivariate (SIRV) models are respectively displayed in red and blue. Moreover, experiments are carried out to evaluate the influence of the standard deviation σ in the decision rule. The solid and dashed lines correspond to the classification results when the standard deviation is considered [START_REF] Amari | Natural gradient works efficiently in learning[END_REF] and when the homoscedasticity assumption holds [START_REF] Arnold | Information-geometric optimization algorithms: A unifying picture via invariance principles[END_REF]. As observed, a gain of about 3 points is observed when a multivariate model (such as the SIRV) is used to take into account the spatial dependency compared to an univariate model (such as the univariate GGD). Note also that the proposed classifier (dash lines) has a significant gain of 5 points compared to the decision rule when the homoscedasticity assumption is considered (solid lines).

CONCLUSION

This paper has addressed the problem of classification based on an intrinsic prior. After introducing the proposed intrinsic prior distribution as the normal law on a Riemannian manifold when the Jeffrey divergence is considered, we have derived an intrinsic estimation and classification scheme. Next an application to supervised classification texture images has been proposed. Classification results on the VisTex database have shown a gain compared to other conventional approaches. Further works will deal with the extension of the proposed work to an intrinsic multi-barycentric classification algorithm in order to handle the intra-class diversity of natural texture images.

ACKNOWLEDGEMENT

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the Investments for the future Programme IdEx Bordeaux (ANR-10-IDEX-03-02).

c

  maximizing the likelihood p(χ t |λ c ), i.e. ĉ = arg max c p(χ t | λc ) (10) By following the same procedure as described in Section. 3.1 for the estimation process, one can rewrite the expression of the decision rule (10) as: (•| θt ), p(•| θc )) (11) Under the homoscedasticity assumption, ∀{i ∈ 1, . . . , c}, σi = σ, (11) reduces to the decision rule ĉ = arg min c J(p(•| θt ), p(•| θ)).

T 13 ) 4 . 2 .

 1342 I = {θ s,o |s = 1, . . . , N s , o = 1, . . . , N o } . (Intrinsic texture estimation Let (T c,1 , . . . , T c,N T r ) be N T r training samples from the same class c. From this collection of samples, the class c is represented by a collection of centroids θc,s,o and standard deviation σ c,s,o computed on each subband, since the subbands of the wavelet decomposition are assumed to be independent. Here, the subscripts c, s, o refer respectively to the texture class c, the scale s and the orientation o of the wavelet subband. Hence, for each wavelet subband, one centroid θc,s,o and standard deviation σ c,s,o are estimated according to the intrinsic estimation scheme developed in Section 3.1, see equations (

Fig. 1 .

 1 Fig. 1. Representation of the intra-class diversity modeled by the proposed intrinsic prior distribution.

  computed on each subband. It yields ĉ = arg min c s,o d ln σc,s,o + 1 2σ 2 c,s,o J(p(•| θt ), p(•| θc,s,o )).

  (•| θt ), p(•| θc,s,o )).

Fig. 2 .

 2 Fig. 2. Evolution of the average kappa index as a function of the number of training samples on the VisTex database for the univariate GGD and the SIRV models.