Maxime Lenormand
email: maxime.lenormand@irstea.fr

Sylvie Huet

Floriana Gargiulo

Guillaume Deffuant

A Universal Model of Commuting Networks

We show that a recently proposed model generates accurate commuting networks on 80 case studies from different regions of the world (Europe and United-States) at different scales (e.g. municipalities, counties, regions). The model takes as input the number of commuters coming in and out of each geographic unit and generates the matrix of commuting flows between the units. The single parameter of the model follows a universal law that depends only on the scale of the geographic units. We show that our model significantly outperforms two other approaches proposing a universal commuting model [1,2], particularly when the geographic units are small (e.g. municipalities).

Introduction

Billions of people move everyday from home to workplace and generate networks of socio-economic relationships that are the vector of social and economic dynamics such as epidemic outbreaks, information flows, city development and traffic [START_REF] Balcan | Multiscale mobility networks and the spatial spreading of infectious diseases[END_REF][START_REF] Ortúzar | Modeling Transport[END_REF]. Understanding the essential properties of these networks and reproducing them accurately is therefore a crucial issue for public health institutions, policy makers, urban development, infrastructure planners, etc. [START_REF] De Montis | The structure of interurban traffic: A weighted network analysis[END_REF][START_REF] De Montis | Modeling commuting systems through a complex network analysis: A study of the italian islands of sardinia and sicily[END_REF]. This challenge is the subject of an intensive scientific activity (see [START_REF] Barthélemy | Spatial networks[END_REF]7] for reviews), in which the analogy of the gravitational attraction inspires a majority of approaches [START_REF] Wilson | Land-use/transport interaction models: Past and future[END_REF][START_REF] Choukroun | A general framework for the development of gravity-type trip distribution models[END_REF]: the number of commuters between two geographic units (cities, counties, regions...) is supposed proportional to the product of the "masses" of each geographic unit (the population for example) and inversely proportional to a function of the distance between them. Unfortunately, numerous experiments showed that the optimum function and parameter values vary a lot with the case studies [START_REF] De Montis | The structure of interurban traffic: A weighted network analysis[END_REF][START_REF] De Montis | Modeling commuting systems through a complex network analysis: A study of the italian islands of sardinia and sicily[END_REF][START_REF] De Vries | Exponential or power distance-decay for commuting? an alternative specification[END_REF][START_REF] Fotheringham | Spatial structure and distance-decay parameters[END_REF]. This situation is not satisfactory because when one wants to generate a particular commuting network without having the total origin destination matrix of commuting, no practical heuristic is available for choosing the adequate type of function and parameter values. This paper addresses this problem.

We consider a recently proposed model [START_REF] Gargiulo | Commuting network model: getting to the essentials[END_REF][START_REF] Lenormand | Generating french virtual commuting network at municipality level[END_REF], differentiating itself from the usual gravity law models in two main features:

• It takes as input the total number of commuters in and out from each geographic unit. With this starting point, the model focuses directly on the influence of the distance between geographic units on the commuting probability. The model is data demanding, but these data are widely available.

• It builds the network progressively, allocating commuters one by one in the different flows, according to probabilities that increase with the number of commuters coming in the destination and decrease with the distance between the origin and destination. These probabilities are updated after each allocation.

Our model is close to the traditional doubly-constrained gravity model [START_REF] Wilson | Land-use/transport interaction models: Past and future[END_REF][START_REF] Choukroun | A general framework for the development of gravity-type trip distribution models[END_REF], but it is more flexible and less data demanding. Indeed, the doubly constrained model and the methods used to solve it require a closed network of commuters: they cannot take into account commuting links outside the considered geographical units. Our individual based stochastic approach overcomes this problem and can deal with the usually available data of total number of commuters in and out of geographic units.

We test this model on 80 case-studies with geographic units of different scales. For example in the same case-study the geographic unit can be either the municipality, the canton or the department, (see an example on Figure 1). More precisely, the case studies include: Czech Republic (municipality scale, 1 case-study), France (municipality scale, 34 case-studies), France (canton scale, 15 case-studies including whole France), France (département scale one case-study (whole France), Italy (municipality scale, 10 case-studies), Italy (province scale, 4 case-studies), USA (county level, 15 case-studies including whole USA). For a detailed description of the datasets see the Supplementary Information Datasets.

We show that the single parameter of our model follows a simple universal law that depends only on the average surface of the considered geographic units. This implies that, given the number of commuters in and out of each geographic unit and their average surface, we can derive the whole matrix of flows with a very good confidence.

Two other approaches [START_REF] Balcan | Multiscale mobility networks and the spatial spreading of infectious diseases[END_REF][START_REF] Simini | A universal model for mobility and migration patterns[END_REF] claim to catch universal properties of commuting networks. We show that our model yields significantly more accurate results, especially for case-studies with small geographic units (e.g. municipalities).

The model

We consider the basic double-constrained model setup, without adding any ingredient about the job market characteristics (professions, salary range, etc.). Instead of solving analytically the optimisation problem, we use an individual based procedure that allocates virtual individuals one by one in the different flows between geographic units, according to a probability that is updated after each allocation.

This individual based approach can deal with less constrained data than the doubly-constrained gravity model that requires the total number of commuters in to be equal to the total number of commuters out. In other words the doubly contrained model can only deal with the flows between the considered geographic units; it cannot take into account the commuting links with destinations outside the case study area. This is a problem when only the numbers of commuters in and out the geographic units are available (and not the complete matrix of the commuting flows), because the data do not distinguish between the flows inside and outside the case study area. It is therefore difficult to estimate the correct data to take as input to the doubly-constrained model in this case. Our approach is more flexible and overcomes this difficulty. It does not require that the total number of commuters in and out to be equal (for more details see [START_REF] Lenormand | Generating french virtual commuting network at municipality level[END_REF]), hence it can easily use directly the usually available data on the number of commuters in and out of each geographic unit.

Let s out i and s in j be respectively the global number of commuters starting from unit u i and the global number of commuters arriving in unit u j . These numbers are initialised from data and then they are progressively modified by the procedure. More precisely, at each step we select unit u i such that s out i > 0 at random, and we consider a virtual commuter starting from u i . We draw at random the working place u j * of this individual among all possible destinations u j according to probabilities P i→j :

P i→j = s in j e -βDij N k=1 s in k e -βD ik (1)
where D ij is the Euclidian distance in meter between units u i and u j (computable from the Lambert or GIS coordinates). Having drawn u j * , we decrement of one s out i and s in j * . Note that decrementing s in and s out at each step complicates significantly the derivation of an analytical expression of the model. We chose a probability decreasing exponentially with the distance, in accordance with the investigations carried out in [START_REF] Lenormand | Generating french virtual commuting network at municipality level[END_REF] and with the literature on commuting network models. The importance of the distance in the commuting choices is embedded in parameter β: for β → 0 the probability tends to be independent from the distance, while for high values of β, the probability tends to zero very rapidly when the distance increases, independently from the number of commuters arriving in the units.

To reduce the border effect (see [START_REF] Lenormand | Generating french virtual commuting network at municipality level[END_REF]), we consider the job-search basin in an extended (EXT) area, composed by the n residential units and m units surrounding the area. Thus, we have n units which are commuting origins and N = n + m units that are commuting destinations. The generated network is saved in matrix T ∈ M n×N (N) where each entry Tij represents the number of commuters between units u i and u j . The algorithm is summarized in Figure 2.

A universal law ruling parameter β

The model depends on a single parameter ruling the importance of the distance in commuting choice. We show that this parameter can be derived as a function of the scale of the problem, independently from the socio-geographical location of the case study area. This opens the possibility to reconstruct the commuting flows (origin-destination matrix) when they are not provided.

We calibrated parameter β by maximising the common part of commuters (CPC), based on the Sørensen index [START_REF] Sørensen | A method of establishing groups of equal amplitude in plant sociology based on similarity of species and its application to analyses of the vegetation on danish commons[END_REF].

CP C(T, T ) = 2N CC(T, T ) N C(T ) + N C( T ) (2)
with:

N CC(T, T ) = n i=1 n j=1 min(T ij , Tij ) N C(T ) = n i=1 n j=1 T ij (3)
where T is the observed origin-destination matrix and T is the simulated one. This is a similarity measure based on the Sørensen index in ecology computing which part of the commuting flows is correctly reproduced, on average, by the simulated network. It varies between 0, when no agreement is found, and 1, when the two networks are identical. We priviledged this indicator because of its direct interpretation. Indeed, when N C(T ) ≃ N C( T ) (it is the case for our model), the CPC represents the percentage of commuting connection correctly located (i.e. with the right pair origin -destination). Moreover, we tested on all case studies that the results obtained with the MAE, the RMSE or CPC 1 are equivalent (see the Supplementary Information Other indicators for more details). As an example on the F R1 case study, Figure 3 shows that the same β value maximizes the CPC and minimizes the MAE. In this Figure we can also note that the CPC is very sensitive to β and that its value does not vary much with the different replicas of the stochastic solving process. Moreover, in order to have an idea of the improvement of the model compared with complete randomness, we have computed the CPC of a random model where the probabilities presented in Eq. ( 1) are uniform (P i→j = 1 n , where n is the number of units). As shown on the Figure 4 we obtained an average CPC around 0.1. For our model, the CPC is always higher than 0.7 with an average around 0.8, which can be interpreted as 70 to 80 % of correctly predicted commuting connections.

Our goal is to derive the value of β from some easily available global characteristics of the case-study, giving the possibility to reconstruct the commuting flows when they are not available. Figure 5 gives strong evidence of such a universal relation.

The x-axis represents the average surface of the geographic units of the case-study ( S in logarithm scale) and the y-axis the optimal β value (in logarithm scale). The linear regression in the log-log plane shows a simple relation:

β = α S -ν (4)
with α = 3.15 • 10 -4 and ν = 0.177. α corresponds to the β value for the unitary surface 1 km 2 . The high value of the adjusted R 2 = 0.92 confirms the quality of the linear model. We observe that β decreases with the average surface of the units S , meaning that, when S is small (e.g. for municipalities in France) the distance is more important in the commuting choice than when S is large (e.g. for regions or counties).

We now evaluate the robustness of our estimation of α and ν using a common statistical procedure: the cross-validation. The cross-validation aims at evaluating the potential error of using the β value derived from the regression model intsead of deriving this value by optimisation for a new case study. This procedure repeats a large number of times the following steps: define a sub-sample of the total sample of case studies, derive a regression model of β from this sub-sample, for each case study that do not belong to the sub-sample, derive β from this regression model and compare the corresponding CPC with the value of β directly calibrated on the complete origin -destination data. The dataset (including 80 case-studies) is randomly cut into two sets, called the training set (comprising 53 case-studies) and the test set (composed of 27 case-studies). We build a regression model on the training set, providing α and ν, from which we derive estimates of β for each of the 27 case-studies of the testing set. We have 27 estimations of β using the relation 4 where α and ν are obtained from the random sub-sample of 53 case-studies. We repeat this process 10, 000 times obtaining 270,000 estimations of β (uniformly distributed over the 80 case-studies) corresponding to about 270,000 80 = 3, 375 estimations of β for each case study. Then we calculate the average, minimum and maximum CPC for each of these values of β, and we compare them with the CPC obtained with value of β directly calibrated on the data.

Figure 4 shows, for each case-study, the CPC associated with the calibrated β, the average CPC obtained with the β values estimated from the cross-validation and the confidence interval defined by the minimum and the maximum values (but it is too small to be seen in most cases). The CPC obtained with the calibrated β value (black triangle) is almost the same as the average CPC obtained with the estimated β in most cases (red square). Globally, we can conclude that the β estimated with the log-linear model and the calibrated β lead to very similar CPCs and also very similar MAE and the RMSE as shown in the appendix Other indicator. The method appears therefore fairly robust and this gives confidence for using it with the value of β derived from our loglog regression in new cases studies.

Comparaison with other universal derivations of commuting networks

Two other different approaches, [START_REF] Balcan | Multiscale mobility networks and the spatial spreading of infectious diseases[END_REF] and [START_REF] Simini | A universal model for mobility and migration patterns[END_REF], claim also to provide a universal derivation of commuting networks. The objective of [START_REF] Balcan | Multiscale mobility networks and the spatial spreading of infectious diseases[END_REF] is to generate a worldwide commuting network, and the model must deal with the wide variety of populations and surfaces of geographic units for which the data are available.

To solve this difficulty, the authors project these data on ad-hoc units defined with a Voronoi diagram. They define their basic unit as a cell approximately equivalent to a rectangle of 25 x 25 kilometers along the Equator. This allows them to calibrate their model because a unit is the same object whatever the country. This is an interesting solution for generating a world-wide commuting network but it leads to an average commuting distance of 250 km which is much larger than the average distance of daily commuting. For example for the USA case study the average distance of daily commuting is about 68 km for the observed network and about 64 km for the simulated network obtained with our algorithm.

For the Auvergne (France) case study at municipality scale the average distance of daily commuting is about 12 km for the observed network and about 11 km for the simulated one.

In the radiation model, proposed in [START_REF] Simini | A universal model for mobility and migration patterns[END_REF], the commuting flow between two geographic units is a function of the cumulated population in a circle at the distance between the two units. The model has an elegant analytical solution and the average flow T ij from unit u i to unit u j can be approximated by

T ij = m i P c P m i n j (m i + s ij ) (m i + n j + s ij ) ( 5
)
where m i and n j are respectively the population of units u i and u j , P c is the total number of commuters and P is the total population in the case-study region, and s ij the total population in the circle of radius r ij centred at u i (excluding the source and destination population). We implemented their analytical approximation and reproduced the graphs presented in their paper. Figure 6 shows the comparison between the radiation model and ours in the US for inter-county commuting and in the French Auvergne region for inter-municipality commuting. We observe that in both cases our approach yields significantly better results. Moreover, as shown on Figure 4, the average CPC for the radiation model on all the case studies is around 0.4, and lower for all case studies than the one obtained with our approach.

However, it should be reminded that our model uses more specific data (total number of commuters in and out of each geographic unit) than the radiation model, hence one could expect our results to be more accurate. Therefore, to be fair with the radiation model we implemented a modified version of this model using the number of out and in commuters of each units. This new approximation is presented in equation 6 where s ij the total number of in-commuters in the circle of radius r ij centred at u i (excluding the source and destination).

T ij = s out i s out i s in j (s out i + s ij ) s out i + s in j + s ij (6)
As shown on Figure 4, this new model reaches an average CPC around 0.5 which is higher than the original radiation model but still significantly lower than the results obtained with our model. Using the MAE and the RMSE leads to the same conclusions (see the appendix Other indicator for more details).

Discussion

The power law of our model's single parameter β with the average area of the case study geographic units, is surprising to us because of the high variety in our case studies in terms of scale, number of units, number of commuters and surface areas. For instance the Auvergne region in France is rural with a population density of about 50 hab./km 2 whereas the New York City region is very urban with a population density of about 6500 hab./km 2 . As far as we know, this is the first time that a single model is shown to fit such diverse group of datasets.

We show that our approach outperforms the radiation model and that the difference of input data plays a minor role in this superiority. This superiority is not due to our particular treatment of the border effects either. Indeed, we could check our approach outperforms the radiation model also on particular case studies (e.g. on islands such as Corsica) where this border effect does not play. We can conclude that the accuracy of our model comes from a proper use of the number of commuters in and out of each geographic unit and an adequate choice of the function of the distance.

The results of the cross validation procedure give a good confidence in the robustness of this law. However, we have to admit that, despite their diversity, our 80 case studies come all from western industrialised countries. Therefore it will be important to check the validity of our law on case studies coming from other continents and less industrialised countries. Moreover, we use a very rough approximation of the distance between the geographic units with the Euclidian distance between the unit centroids. More accurate approximations of this distance would certainly improve the results. Finally, we also intend to apply our approach to commuting networks inside urban areas because many cities of the world show an impressive growth and an increasing part of commuting takes place within them [START_REF] Roth | Structure of urban movements: Polycentric activity and entangled hierarchical flows[END_REF]. An important issue in our perspective is to check if our law holds at this scale.

∈ M n×N (R), s in ∈ N N , s out ∈ N n , β ∈ R + Output : T ∈ M n×N (N) Tij ← 0 while n k=1 s out k > 0 do Pick at random i ∈ |[1, n]|, such that s out i = 0 Pick at random j from |[1, N ]| with a probability P i→j Tij ← Tij + 1 s in j ← s in j -1 s out i ← s out i -1 end
CZ FR1 FR2 FR3 FR4 FR5 FR6 FR7 FR8 FR9 FR10 FR11 FR12 FR13 FR14 FR15 FR16 FR17 FR18 FR19 FR20 FR21 FR22 FR23 FR24 FR25 FR26 FR27 FR28 FR29 FR30 FR31 FR32 FR33 FR34 FRc0 FRc1 FRc2 FRc3 FRc4 FRc5 FRc6 FRc7 FRc8 FRc9 FRc10 FRc11 FRc12 FRc13 FRc14 FRd0 IT1 IT2 IT3 IT4 IT5 IT6 IT7 IT8 IT9 IT10 ITp0 ITp1 ITp2 ITp3 USA0 USA1 USA2 USA3 USA4 USA5 USA6 USA7 USA8 USA9 USA10 USA11 USA12 USA13 USA14

Figure 1 .

 1 Figure 1. Three scales of geographic units (Auvergne region, France)

Figure 3 .

 3 Figure 3. Plot of the average CPC (blue circle) and the average NMAE (red triangle) in term of β for 10 replications of the model for the Auvergne case study (FR1). The error bars represent the minimum value and maximum value obtain over the 10 replications.

Figure 4 .

 4 Figure 4. Common part of commuters (CPC) for the 80 case-studies. The red squares represent the CPC obtained with the value of β optimised from data on the case-study network. Black plain triangles represent the average CPC obtained with β values estimated with the rule linking β and the average surface of the units obtain with the cross-validation; Dark bars represent the minimum and the maximum CPC obtained with the estimated β but in most cases they are too close to the average to be seen. The green circles represent the CPC obtained with the random model. The blue triangles represent the CPC obtained with the radiation model. The purple crosses represent the CPC obtained with the modified version of the radiation model.

Figure 5 .Figure 6 .

 56 Figure 5. Log-log scatter plot of the calibrated β values in terms of average surface of the geographic units for 80 case-studies; the line represents the regression line predicting β. The surface is made non-dimensional by the unitary surface 1 km 2 .

Available at http://www.plosone.org

We have also shown in[START_REF] Gargiulo | Commuting network model: getting to the essentials[END_REF][START_REF] Lenormand | Generating french virtual commuting network at municipality level[END_REF] that the value of β yielding the maximum CPC also yields the maximum similarity between observed and simulated commuting distance distributionsAuthor-produced version of the paper Lenormand M, Huet S, Gargiulo F, Deffuant G (2012) A Universal Model of Commuting Networks. PLoS ONE 7(10): e45985. doi:10.1371/journal.pone.0045985 Available at http://www.plosone.org

Author-produced version of the paper Lenormand M, Huet S, Gargiulo F, Deffuant G (2012) A Universal Model of Commuting Networks. PLoS ONE 7(10): e45985. doi:10.1371/journal.pone.0045985 Available at http://www.plosone.org

Acknowledgments References


Dmitry Berdinsky

Min-Jae Oh

Tae-Wan Kim

Bernard Mourrain

On the problem of instability in the dimension of a spline space over a T-mesh

Keywords: Spline space, T-mesh, Dimension

In this paper, we discuss the problem of instability in the dimension of a spline space over a T-mesh. For bivariate spline spaces S (5, 5, 3, 3) and S (4, 4, 2, 2), the instability in the dimension is shown over certain types of T-meshes. This result could be considered as an attempt to answer the question of how large the polynomial degree (m, m ′ ) should be relative to the smoothness (r, r ′ ) to make the dimension of a spline space stable. We show in particular that the bound m ≥ 2r + 1 and m ′ ≥ 2r ′ + 1 are optimal.

Introduction

This paper is motivated by a recent work [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] devoted to the observation of instability in the dimension of a spline space S (m, m, m -1, m -1) over a T-mesh.

Splines over a T-mesh [START_REF] Deng | Dimension of spline spaces over T-meshes[END_REF] could be a useful tool in many areas such as surface modeling and finite element analysis. For instance, polynomial splines over hierarchical T-meshes (PHT-splines) [START_REF] Deng | Polynomial splines over hierarchical T-meshes[END_REF] appear to be a natural generalization of B-splines over hierarchical Tmeshes. A method has been proposed for stitching several surface patches and for fitting genus-zero meshes using PHT-splines [START_REF] Li | Surface modeling with polynomial splines over hierarchical T-meshes[END_REF]. Splines over T-meshes allow one to avoid the problem of superfluous control points, because inserting an extra node implies only local refinement of a control mesh. From this point of view, splines over T-mesh also have a good prospect for using in isogeometric analysis [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF]; for instance, see [START_REF] Nguyen-Thanh | Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids[END_REF], [START_REF] Vuong | A hierarchical approach to adaptive local refinement in isogeometric analysis[END_REF].

However, unfortunately, the problem of describing the basis functions of a spline space over a given Tmesh is generally non-trivial. For bicubic splines with continuity of order one, an efficient algorithm has been proposed to construct basis functions over a general Tmesh [START_REF] Li | Polynomial splines over general Tmeshes[END_REF]; these basis functions satisfy remarkable properties such as local support, nonnegativity, and partition of unity. However, for bicubic splines with continuity of order two, even the dimension is not easy to calculate, and it is exactly known only for some particular classes of T-meshes. Therefore, the dimension is one of the crucial characteristics of a spline space over a T-mesh that must be analyzed.

We briefly recall some basic notations and recent results concerning the dimension of a spline space over a T-mesh.

A T-mesh is basically a subdivision of some domain into rectangular cells that allow T-junctions (see Fig. 1). In this paper, only a rectangular domain is considered. Let S (m, m ′ , r, r ′ )(T ) denote the space of functions of variables x, y that are piecewise polynomials of degree m, m ′ and the class of smoothness C r,r ′ over the given Tmesh T . We also refer to S (m, m ′ , r, r ′ )(T ) as the spline space over T-mesh T .

In [START_REF] Deng | Dimension of spline spaces over T-meshes[END_REF], the exact formula for the dimension of the space S (m, m ′ , r, r ′ )(T ) is obtained when

m ≥ 2r + 1, m ′ ≥ 2r ′ + 1. In this case, dim S (m, m ′ , r, r ′ )(T ) = (m + 1)(m ′ + 1) f 2 - (m + 1)(r ′ + 1) f h 1 -(m ′ + 1)(r + 1) f v 1 + (r + 1)(r ′ + 1) f 0 , (1)
where f 2 is the number of cells of T , f h 1 and f v 1 respectively are the numbers of horizontal and vertical interior edges (these intersect the interior of the rectangular domain of T ), and f 0 is the number of interior vertices (these are in the interior of the rectangular domain). For instance, in Fig. 1, ν 1 , . . . , ν 7 are interior vertices; b 2 ν 2 , ν 1 ν 4 , vertical interior edges; and b 4 ν 1 , ν 6 ν 7 , horizontal interior edges.

The method used to calculate the dimension of S (m, m ′ , r, r ′ )(T ) in [START_REF] Deng | Dimension of spline spaces over T-meshes[END_REF] is based on Bézier nets. In other words, in [START_REF] Deng | Dimension of spline spaces over T-meshes[END_REF], the formula for the number of elements in a minimal determining set of Bézier ordinates is obtained.

Proofs of analogous results have been presented in [START_REF] Huang | New proof of dimension formula of spline spaces over T-meshes via smoothing cofactors[END_REF] and [START_REF] Li | Improvement on the dimension of spline spaces on T-mesh[END_REF] independently. However, the approach used in these studies is based on the smoothing cofactor method.

All these methods are based on the following simple observation. Let P 1 (x, y) and P 2 (x, y) be two polynomials of degree m, m ′ defined on rectangular domains Ω 1 and Ω 2 , respectively, and having the common interval (ν 1 , ν 2 ) (see Fig. 2). Then, the piecewise polynomial defined on Ω 1 Ω 2 has r-th order of smoothness for variable x, if and only if

P 1 (x, y) -P 2 (x, y) = λ(x, y)(x -x 0 ) r+1 ,
for some polynomial λ(x, y) of degree mr -1, m ′ for variables x, y, respectively.

Here, we do not provide any detailed explanations of the approach used in these studies, but only mention the fact that the smoothing cofactor method appears to be the basic method for all further results concerning the dimension of a spline space over a T-mesh.

Thus, when m ≥ 2r + 1, m ′ ≥ 2r ′ + 1, one can express the dimension formula in simple topological quantities such as the number of cells, edges, and vertices. This fact implies that any small variation, saving the topology of T-mesh T , does not change the dimension of the spline space S (m, m ′ , r, r ′ )(T ). We will show hereafter that the conditions m ≥ 2r+1, m ′ ≥ 2r ′ + 1 are optimal by giving explicit cases where m = 2 r, m ′ = 2 r ′ or m = 2 r -1, m ′ = 2 r ′ -1 for which the dimension of the spline space changes with the position of the vertices.

In [START_REF] Li | Improvement on the dimension of spline spaces on T-mesh[END_REF], Li et al. made the following improvement to the dimension formula. They found the constraint on the T-mesh that guarantees the dimension of the corresponding spline space to be calculated using formula [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF].

It appears that the last progress in the analysis of the dimension of a spline space, in a general case, was made in [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF]. In this study, the results presented in [START_REF] Huang | New proof of dimension formula of spline spaces over T-meshes via smoothing cofactors[END_REF] and [START_REF] Li | Improvement on the dimension of spline spaces on T-mesh[END_REF] were reproved. In addition to the smoothing cofactor method, [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] exploited homology techniques to provide additional insight on the analysis of the dimension of the spline spaces. Finally, the dimension of a space S (m, m ′ , r, r ′ )(T ) is represented in this study in the following manner:

dim S (m, m ′ , r, r ′ )(T ) = (m + 1)(m ′ + 1) f 2 - (m + 1)(r ′ + 1) f h 1 -(m ′ + 1)(r + 1) f v 1 + (r + 1)(r ′ + 1) f 0 + h r,r ′ m,m ′ (T ), (2)
where the meaning of the term h r,r ′ m,m ′ (T ) is described in the next section.

Recall that the hierarchical T-mesh is either the initial square or is obtained from a hierarchical T-mesh by splitting a cell along a vertical or horizontal line. A valuable observation was made in [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] about the upper estimation of the summand h r,r ′ m,m ′ (T ) for a hierarchical T-mesh T . This estimation allows one to establish the constraint on a hierarchical T-mesh T for the summand h r,r ′ m,m ′ (T ) to be zero, and this constraint, in some manner, generalizes the one obtained in [START_REF] Li | Improvement on the dimension of spline spaces on T-mesh[END_REF]. Therefore, if a hierarchical T-mesh T satisfies the constraint proposed in [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF], then the dimension of a spline space S (m, m ′ , r, r ′ )(T ) is given by [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF].

In a general case, the formula for the dimension of a space S (m, m ′ , r, r ′ )(T ) cannot be expressed using only the simple topological characteristics of T , such as f 2 , f h 1 , f v 1 , and f 0 . Recently, Li and Chen [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] demonstrated that the dimension of a space S (m, m, m -1, m -1)(T ) may also depend on the geometry of a T-mesh T .

For instance, the dimension of a space S (3, 3, 2, 2)(T ) for T , shown in Fig. 3, depends on the values of the coordinates s i , t i , i = 1 . . . 6. There are certain values of the coordinates s i , t i , i = 1 . . . 6 for which the dimension jumps with a small variation in these values; this indicates the instability in the dimension of space S (3, 3, 2, 2)(T ).

The interesting question stated in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] is how large the polynomial degree should be relative to the smoothness for the dimension of a spline space to be stable. In this work, we treat the new cases, S (4, 4, 2, 2)(T ) and S (5, 5, 3, 3)(T ), and provide certain examples of Tmeshes T with instability in the dimension of the spline space. This shows that the constraints m ≥ 2 r + 1, m ′ ≥ 2 r ′ + 1 known to yield a combinatorial dimension formula are optimal.

We note that the key approach to construct examples of T-meshes with instability in the dimension was in using the abovementioned constraint on a T-mesh. The key was to violate this constraint to guarantee that the dimension is not given by [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF]. Unfortunately, the summand h r,r ′ m,m ′ (T ) for some T-mesh T should be calculated as a rank of the corresponding matrix. However, it is difficult to control even for small degrees m, m ′ . Therefore, we specified coordinates of some nodes to simplify the calculations. As shown in the next section, the number of rows of this matrix linearly depends on the number of maximal interior segments. Therefore, to minimize the matrix size and make its rank unstable, we treated T-meshes with four maximal interior segments.

We investigated the cases S (4, 4, 2, 2)(T ) and S (5, 5, 3, 3)(T ) because after the result of [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF], they were the first natural candidates to test the stability in the dimension. The other motivation was that splines with small polynomial degree could be of interest in many applications.

We also assume that further generalizations of our results require more sophisticated techniques. We believe that investigating the instability in the dimension of a space S (m, m ′ , r, r ′ )(T ) is an important issue because it provides a better understanding of the main problem, namely, how the basis of a space S (m, m ′ , r, r ′ )(T ) can be described.

The remainder of this work is organized as follows.

In section 2, we explain the dimension formula (2) in detail and interpret it in a form pertinent for showing the instability in the dimension. In section 3, we explain how the dimension formula (2) can be used to prove the instability in the dimension of a spline space S (3, 3, 2, 2)(T ) for certain types of T-meshes T described in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF]. In section 4, we describe certain types of T-meshes T and prove the instability in the dimension of spline spaces S (5, 5, 3, 3)(T ) and S (4, 4, 2, 2)(T ). In section 5, we present the conclusions of this work and discuss future works.

Dimension of a spline space S(m, m ′ , r, r ′ )(T) over T-mesh T: General formula

In this section, we explain the dimension formula (2) in detail. We also transform it into a form pertinent for showing the instability in the dimension. For reference, one can also find similar types of formulas for the dimension of a spline space in previous works [START_REF] Huang | New proof of dimension formula of spline spaces over T-meshes via smoothing cofactors[END_REF], [START_REF] Li | Improvement on the dimension of spline spaces on T-mesh[END_REF].

The general formula (2) for the dimension of the space S (m, m ′ , r, r ′ )(T ) for some T-mesh T of the rectangular domain is found by using the homology technique [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF].

The summand h r,r ′ m,m ′ (T ) in (2) may be considered as the dimension of the factor space M K [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF]. The space M is the direct sum

M = ρ∈MIS (T ) [ρ]R (m,m ′ )-δ(ρ)
over the entire set of maximal interior segments MIS (T ).

The maximal interior segment is the maximal union of connected edges of the same direction, such that there is no intersection with the border of the domain (for example, the segments ν 1 ν 6 and ν 4 ν 5 are all maximal interior segments of the T-mesh shown in Fig. 1).

Here, R (m,m ′ ) denotes the space of polynomials of degree m, m ′ of variables x, y, respectively, and δ(ρ) = (0, r ′ + 1) if ρ is the horizontal segment and δ(ρ) = (r + 1, 0) if ρ is the vertical segment.

One can consider the space M as simply a vector space formed by the formal sums Σ ρ∈MIS (T ) β ρ [ρ], where β ρ ∈ R (m,m ′ -r ′ -1) if ρ is the horizontal maximal interior segment, and

β ρ ∈ R (m-r-1,m ′ ) if ρ is the vertical maximal interior segment.
The subspace K ⊂ M has the following form:

K = Σ γ∈T 0 0 (∆ r,r ′ h [ρ v (γ)] -∆ r,r ′ v [ρ h (γ)])R (m-r-1,m ′ -r ′ -1) ,
where T 0 0 is the set of vertices of maximal interior segments (for instance ν 1 , ν 4 , ν 5 , and ν 6 are all points of T 0 0 for the T-mesh shown in Fig. 1). In addition,

∆ r,r ′ h = (y -y γ ) r ′ +1 , ∆ r,r ′ v = (x -x γ ) r+1
, where x γ , y γ denotes the x, y coordinates of the vertex γ. Furthermore, ρ h (γ) and ρ v (γ) are the corresponding horizontal and vertical maximal interior segments containing the vertex γ (for example, in Fig. 1, ρ h (ν 4 ) = ν 4 ν 5 , ρ v (ν 4 ) = ν 1 ν 6 ). For more details, see [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF].

To investigate the question of instability in the dimension of a space S (m, m ′ , r, r ′ )(T ) over some T-mesh T , it is only necessary to consider the summand h r,r ′ m,m ′ (T ). Let n be the number of vertices of T 0 0 , and l h , l v be the numbers of horizontal and vertical maximal interior segments, respectively. One can interpret h r,r ′ m,m ′ (T ) for some T-mesh T in terms of the rank of the corresponding matrix. Indeed,

h r,r ′ m,m ′ = dimM -dimK = (m + 1)(m ′ -r ′ )l h + (m ′ + 1)(m -r)l v -(m -r)(m ′ -r ′ )n+ dimN.
where N is the subset

N ⊂ γ∈T 0 0 R (m-r-1,m ′ -r ′ -1) ,
so that each element of N is a set of polynomials

α γ ∈ R (m-r-1,m ′ -r ′ -1)
, γ ∈ T 0 0 , such that it satisfies the following conditions:

Σ γ∈ρ h α γ (x -x γ ) r+1 = 0 (3)
for each horizontal maximal interior segment ρ h , and

Σ γ∈ρ v α γ (y -y γ ) r ′ +1 = 0 (4)
for each vertical maximal interior segment ρ v . However, the space N can be represented as the space of solutions of the system of linear equations Lv = 0, where v is the column vector of coefficients of polynomials α γ , γ ∈ T 0 0 . Here, we do not describe the structure of the matrix L in the general case because it is simpler to understand it from the examples below.

We simply note that the number of columns in matrix L is (mr)(m ′ -r ′ )n and the number of rows is (m

+ 1)(m ′ -r ′ )l h + (m ′ + 1)(m -r)l v . One can compute dim N as dim N = (m -r)(m ′ -r ′ )n -rank L.
Thus, for some T-mesh T,

h r,r ′ m,m ′ (T ) = (m + 1)(m ′ -r ′ )l h + (m ′ + 1)(m -r)l v -rank L. ( 5
)

Example of T-mesh with instability in the dimension of a spline space S(3, 3, 2, 2)(T)

In order to establish the instability in the dimension of a spline space, it is necessary to find a T-mesh in which the dimension jumps for some small change in the mesh. The proof of dimension jumping is based on calculating the rank of the corresponding matrix L. For m = m ′ , r = r ′ = m -1, examples of such T-meshes are found in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF].

We describe in detail the example given in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] for the case m = m ′ = 3, r = r ′ = 2. For the T-mesh shown in Fig. 3, every element of the space N (see (3),( 4)) is simply the 16-th dimensional column vector v of the real numbers α j , j = 1 . . . 16 (each α j corresponds to the vertex γ j , j = 1 . . . 16), such that the following equations hold for arbitrary x and y:

Σ 5 i=1 α i (x -x γ i ) 3 = 0, Σ 9 i=5 α i (y -y γ i ) 3 = 0, Σ 13 i=9 α i (x -x γ i ) 3 = 0, Σ 16
i=13 α i (yy γ i ) 3 + α 1 (yy γ 1 ) 3 = 0. Note that these equations correspond to the maximal interior segments γ 2 γ 4 , γ 6 γ 8 , γ 10 γ 12 , and γ 14 γ 16 , respectively.

Because the coefficients of the polynomials in these equations are zero, one can obtain the following system of linear equations:

α 2 s k 1 + α 3 s k 2 + α 5 s k 3 + α 1 s k 4 + α 4 s k 5 = 0, 4 α 6 t k 1 + α 7 t k 2 + α 5 t k 3 + α 9 t k 4 + α 8 t k 5 = 0, α 10 s k 2 + α 9 s k 3 + α 13 s k 4 + α 11 s k 5 + α 12 s k 6 = 0, α 14 t k 2 + α 1 t k 3 + α 13 t k 4 + α 15 t k 5 + α 16 t k 6 = 0, for k = 0 . . . 3.
By writing these linear equations in the form Lv = 0, we obtain the following form for the 16 × 16 matrix L:

      1 1 1 1 1 0 0 0 s 4 s 1 s 2 s 5 s 3 0 0 0 s 2 4 s 2 1 s 2 2 s 2 5 s 2 3 0 0 0 s 3 4 s 3 1 s 3 2 s 3 5 s 3 3 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 t 3 t 1 t 2 t 5 0 0 0 0 t 2 3 t 2 1 t 2 2 t 2 5 0 0 0 0 t 3 3 t 3 1 t 3 2 t 3
5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 t 3 0 0 0 0 0 0 0 t 2 3 0 0 0 0 0 0 0 t 3 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 t 4 0 0 0 0 0 0 0 t 2 4 0 0 0 0 0 0 0 t 3 4 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 s 3 s 2 s 5 s 6 s 4 0 0 0 s 2 It is shown in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] that the determinant of the matrix L is zero if and only if

(s 3 -s 1 )(s 6 -s 4 ) (t 3 -t 1 )(t 6 -t 4 ) = (s 4 -s 1 )(s 6 -s 3 ) (t 4 -t 1 )(t 6 -t 3 ) .
This equation holds, for instance, when s i = t i , i = 1 . . . 6. However, a small change in coordinates s i , t i , i = 1 . . . 6 makes the matrix L non-degenerate, which implies jumping of the rank(L), and therefore, because of ( 2) and ( 5) jumping of the dimension of the spline space. The bound of [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] yields in this case 0 ≤ h 2,2 3,3 (T ) ≤ 1.

Note that the structure of the T-meshes in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF] is not complicated, and these T-meshes only have four maximal interior segments.

The observation of instability in the dimension of S (m, m, m -1, m -1)(T ), given in [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF], inspired us to construct examples of T-meshes with instability in the dimension for the cases S (4, 4, 2, 2)(T ) and S (5, 5, 3, 3)(T ).

Examples of T-meshes with instability in the dimension of spline spaces S(4, 4, 2, 2)(T) and S(5, 5, 3, 3)(T)

In this section, we describe certain types of T-meshes T and prove the instability in the dimension of spline spaces S (4, 4, 2, 2)(T ) and S (5, 5, 3, 3)(T ).

We start with the case S (5, 5, 3, 3)(T ), because for the example of the T-mesh given for S (2, 2, 1, 1)(T ) [START_REF] Li | On the instability in the dimension of splines spaces over T-meshes[END_REF], instability in the dimension of the space S (5, 5, 3, 3)(T ) has been demonstrated. This T-mesh is shown in Fig. 4.

Case: m

= m ′ = 5, r = r ′ = 3
= m ′ = 2, r = r ′ = 1.
The polynomials α γ i ∈ R (1,1) , γ i ∈ T 0 0 , i = 1 . . . 12 can be represented as

α γ i = k 1i (x -x γ i )(y -y γ i ) + k 2i (x -x γ i ) +k 3i (y -y γ i ) + k 4i .
Let v denote the column vector of the length 48, such that for j = 4 l + r, r = j mod 4, the j-th coordinate of v equals k (r+1)(l+1) . Then, the set of polynomials α γ i , i = 1 . . . 12 is in N (see [START_REF] Deng | Polynomial splines over hierarchical T-meshes[END_REF] and ( 4)) if and only if the following polynomials are zero:

Σ γ i ∈ρ k 1i (x -x γ i ) 5 + k 3i (x -x γ i ) 4 = 0, Σ γ i ∈ρ k 2i (x -x γ i ) 5 + k 4i (x -x γ i ) 4 = 0 ( 6
)
if ρ is the maximal horizontal interior segment, and

Σ γ i ∈ρ k 1i (y -y γ i ) 5 + k 2i (y -y γ i ) 4 = 0, Σ γ i ∈ρ k 3i (y -y γ i ) 5 + k 4i (y -y γ i ) 4 = 0 ( 7
)
if ρ is the maximal vertical interior segment. Let V t denote, for some given t, the rectangular 12×4 matrix of the form

V t =                                                          1 0 0 0 t -1 5 0 0 t 2 -2 5 t 0 0 t 3 -3 5 t 2 0 0 t 4 -4 5 t 3 0 0 t 5 -t 4 0 0 0 0 1 0 0 0 t -1 5 0 0 t 2 -2 5 t 0 0 t 3 -3 5 t 2 0 0 t 4 -4 5 t 3 0 0 t 5 -t 4       1 0 0 0 s 0 -1 5 0 s 2 0 -2 5 s 0 s 3 0 -3 5 s 2 0 s 4 0 -4 5 s 3 0 s 5 0 -s 4 0 0 1 0 0 0 s 0 -1 5 0 s 2 0 -2 5 s 0 s 3 0 -3 5 s 2 0 s 4 0 -4 5 s 3 0 s 5 0 -s 4               k 1i k 2i k 3i k 4i               = 0 ( 8
)
if ρ is the maximal horizontal interior segment, and

Σ γ i ∈ρ V y γ i               k 1i k 2i k 3i k 4i               = 0 ( 9
)
if ρ is the maximal vertical interior segment. We consider the system of linear equations ( 8) and ( 9) for coefficients k ji , j = 1 . . . 4, and i = 1 . . . 12 in terms of the linear system Lv = 0. Then, the matrix L has 4 × 12 = 48 columns and 12 × 4 = 48 rows. By using equations ( 8) and ( 9) and sequentially taking the maximal interior segments γ 2 γ 3 , γ 5 γ 6 , γ 7 γ 9 , and γ 1 γ 11 (see Fig. 4), one can represent the matrix L in terms of 12 × 4 blocks H x γ i , V y γ i , i = 1 . . . 12 and zero blocks 0 filled with nulls in the following manner:

              V t 3 V t 2 V t 5 V t 4 0 0 0 0 0 H s 4 H s 5 H s 2 0 0 0 0 0 0 H s 4 0 0 0 0 0 0 0 0 0 0 0 H s 3 0 0 0 0 0 V t 4 V t 2 V t 1 V t 3 0 0 0 0 0 H s 3 H s 1 H s 2              
Thus, according to [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], for the given T-mesh, h

equals 48 -rank (L).

In this work, we do not consider the question of calculating the rank of L for an arbitrary set of values t 1 , . . . , t 5 and s 1 , . . . , s 5 to describe all possible positions with instability in the dimension.

We are simply interested in the topology of the T-mesh and in providing some concrete values for t 1 , . . . , t 5 and s 1 , . . . , s 5 such that dimension jumping occurs for some small change in these values.

For instance, let us consider the following values:

t 1 = s 1 = -2, t 2 = s 2 = -1, t 3 = s 3 = 0, s 4 = t 4 = 1, s 5 = t 5 = 2.
Then, rank L = 46. In order to be both rigorous and simple, we fix all values except t 1 and calculate the determinant of the matrix L. One can check that the determinant is nonzero when t 1 -2, implying that the rank of L is equal to 48. Therefore, the dimension of S (5, 5, 3, 3) is increased by 2 when t 1 = -2, and therefore, we have dimension jumping.

Notice that the bound of [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] yields 0 ≤ h 3,3 5,5 (T ) ≤ 4.

Case: m

= m ′ = 4, r = r ′ = 2.
For S (4, 4, 2, 2)(T ), the example is a bit more curious. It is shown in Fig. 5.

The approach used in this subsection is the same as that used in the previous one. Again, the polynomials α γ i ∈ R (1,1) , γ i ∈ T 0 0 , i = 1 . . . 10 can be represented as

α γ i = k 1i (x -x γ i )(y -y γ i ) + k 2i (x -x γ i ) + k 3i (y -y γ i ) + k 4i .
Let v denote the column vector of length 40, such that for j = 4 l + r, r = j mod 4, the j-th coordinate of v equals k (r+1)(l+1) . Then, the set of polynomials α γ i , i = 1 . . . 10 is in space N (see (3) and ( 4)) if and only if the following polynomials are zero:

Σ γ i ∈ρ k 1i (x -x γ i ) 4 + k 3i (x -x γ i ) 3 = 0, Σ γ i ∈ρ k 2i (x -x γ i ) 4 + k 4i (x -x γ i ) 3 = 0 ( 10
)
if ρ is the maximal horizontal interior segment, and

Σ γ i ∈ρ k 1i (y -y γ i ) 4 + k 2i (y -y γ i ) 3 = 0, Σ γ i ∈ρ k 3i (y -y γ i ) 4 + k 4i (y -y γ i ) 3 = 0 ( 11
)
if ρ is the maximal vertical interior segment. We define the 10 × 4 matrices V t and H s as follows:

     V t 2 V t 1 V t 4 V t 3 0 0 0 0 H s 5 H s 1 0 0 0 0 0 H s 5 0 0 0 0 0 0 0 0 0 H s 3 H s 2 0 0 0 0 V t 3 V t 1 V t 2 0 0 0 0 H s 2 H s 4              
We take the following values for t 1 , t 2 , t 3 , t 4 and s 1 , s 2 , s 5 :

t 1 = 0, t 2 = 1, t 3 = 3, t 4 = 4, s 1 = 0, s 2 = 1, s 5 = 3.
If we substitute s 3 = s 4 , then the rank of L is 38, otherwise it is 39. This means that the dimension of S (4, 4, 2, 2) for the given T-mesh increases by 1 when s 3 = s 4 . Again, the bound of [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] yields 0 ≤ h 2,2 4,4 (T ) ≤ 4.

Note that both examples have only four maximal interior segments and their structures are quite simple.

Conclusions and future work

In this work, we proved the instability in the dimension of a spline space S (m, m ′ , r, r ′ To prove the instability in the dimension, we specified the values of the node coordinates of the T-mesh for the following two reasons. First, our main focus is on understanding the topology of a T-mesh with instability in the dimension, and second, calculating the rank of the matrix L is simplified considerably. This work does not accurately analyze the algebraic conditions for the node coordinates for the dimension to be unstable. However, such an analysis could be carried out by implementing a row reduction procedure to calculate the rank of matrix L.

In the future, the following problems would be interesting to explore. The first is the investigation of the instability in the dimension for the following, more general cases: m = m ′ , r = r ′ and 2r + 1 > m. The problem is that the analysis of the rank of matrix L is nontrivial in more general cases; therefore, we might need to use another interpretation of the dimension formula. The second is the description of classes of T-meshes for the dimension of a space S (m, m ′ , r, r ′ )(T ) to be stable and controlled by some topological characteristics. The work of [START_REF] Mourrain | On the dimension of spline spaces on planar T-subdivisions[END_REF] made progress in this regard, but we believe that it could be considerably strengthened by a more careful analysis of h r,r ′ m,m ′ . We intend to investigate these questions in our future work.

10 bFigure 1 :

 101 Figure 1: Example of a T-mesh.

Figure 2 :

 2 Figure 2: Two adjacent cells.

Figure 3 :

 3 Figure 3: Unstable T-mesh for m = m ′ = 3, r = r ′ = 2.

Figure 4 :

 4 Figure 4: Unstable T-mesh for cases m = m ′ = 5, r = r ′ = 3 and m = m ′ = 2, r = r ′ = 1.

Figure 5 :

 5 Figure 5: Unstable T-mesh for m = m ′ = 4, r = r ′ = 2

  )(T ) over certain types of T-meshes T for the new cases m = m ′ = 4, r = r ′ = 2 and m = m ′ = 5, r = r ′ = 3. The constructed examples have only four maximal interior segments and appear to be the simplest ones in some sense. One can construct more complicated examples of T-meshes with instability in the dimension by some combinations of the T-meshes shown in Figs. 4 and 5 for the cases m = m ′ = 5, r = r ′ = 3 and m = m ′ = 4, r = r ′ = 2, respectively.

Acknowledgement

This research was supported by the National Research Foundation of Korea(NRF) grant funded by the Korea government(MEST) (No. 2011-0018023) and partially supported by the Industrial Strategic technology development program, 10035474, Development of inspection platform technology based on 3-dimensional X-ray images funded by the Ministry of Knowledge Economy(MKE, Korea).


Matthieu Gautier
email: matthieu.gautier@irisa.fr1

Vincent Berg
email: vincent.berg@cea.fr

Dominique Noguet
email: dominique.noguet@cea.fr

Wideband frequency domain detection using Teager-Kaiser energy operator

This paper addresses wireless microphone sensing in the TV white space and efficient detection of narrowband FM modulation signals. To this end, a wideband frequency domain analysis is proposed. The required Fast Fourier Transform for this operation may be shared between sensing analysis and modulation functions. A particular decision metric is then studied for the analysis of wireless microphone signals based on the Teager-Kaiser energy operator. Simulation results show that 6 dB of detection gain could be achieved when using a frequency domain analysis compared to time domain methods. The Teager-Kaiser detection leads to further improvement of 1.5 dB. This performance could be reached at no extra cost in term of complexity.

I. INTRODUCTION 1 T

The performance of sensing algorithms is fundamental to establish the opportunistic communication of a cognitive radio (CR) system [START_REF] Mitola | Cognitive radio: making software radios more personal[END_REF]. The UHF band is a under-used candidate band where the CR system could operate [START_REF]Official announcement of FCC[END_REF]. In this band, the primary users are TV transmitters and wireless microphones. If the detection of TV signals has been addressed in the literature, solutions should be proposed for the detection of wireless microphones [START_REF] Chen | Spectrum sensing for wireless microphone signals[END_REF]. Wireless microphone signals are narrowband in comparison to their overall possible band of operation. This makes wireless microphone sensing difficult. Indeed, the European TV band is composed of 48 channels of 8 MHz bandwidth. Each band has to be analyzed and TV and microphone signals must be detected. A statistic test should be provided to the CR system where the final decision on occupancy assessment is made. The analyzed TV channel is extremely wide compared to the frequency band occupied by the microphone signal (i.e.: around 200kHz). Thus, the detection of wireless microphones may be considered a narrowband signal detection problem. Most of the references available in the literature use a blind detection for the specific case of wireless microphones [START_REF] Chen | Spectrum sensing for wireless microphone signals[END_REF]. These studies are based on eigenvalue decomposition [START_REF] Xu | SVD based Sensing of a Wireless Microphone Signal in Cognitive Radio Networks[END_REF], spectral correlation [START_REF] Han | Spectral correlation based signal detection method for spectrum sensing in IEEE 802.22 WRAN systems[END_REF] and energy detection [START_REF] Ghosh | Spectrum Sensing Prototype for Sensing ATSC and Wireless Microphone Signals[END_REF]. One common property of these algorithms is that they assume the detection of wideband signals, i.e.: the analyzed bandwidth is of the same order as the signal bandwidth. In this paper, the proposed architecture considers the narrowband property of the wireless microphone signal in the TVWS and proposes sharing a wideband Fast Fourier Transform 1 Matthieu Gautier is currently an associate professor at the CAIRN team of the IRISA laboratory. The study and results presented in this paper have been carried out when he was a research engineer at CEA-LETI.

(FFT) operator that can be an element used by the cognitive device to simultaneously detect primary users and make the communication link. Indeed, the concept of "spectrum pooling" has been proposed as a candidate modulation scheme for CR system [START_REF] Weiss | Spectrum pooling: an innovative strategy for the enhancement of spectrum efficiency[END_REF]. It proposes to use a wideband FFT to access an important bandwidth and allocate frequency bins to CR users according to white spaces opportunities. The same FFT operation may also be used to perform the sensing step [START_REF] Noguet | Advances in Opportunistic Radio Technologies for TVWS[END_REF]. The outputs of the FFT are analyzed using specific frequency domain metrics. This leads to a wideband sensing operation and a low complexity architecture with shared FFT operation. Furthermore, a frequency domain metric is specificaly proposed to detect wireless microphone signals. The authors proposed in [START_REF] Gautier | Teager-Kaiser energy detector for narrowband wireless microphone spectrum sensing[END_REF] the Teager-Kaiser detector that allows a better estimation of the energy of the wireless microphone signal. Introduced in the time domain [START_REF] Gautier | Teager-Kaiser energy detector for narrowband wireless microphone spectrum sensing[END_REF], a frequency domain counterpart of this operator is introduced and analyzed in this paper. This paper consists of 5 parts. Section II gives the models of the wireless microphone signals and of the CR system. In Section III, the wideband frequency domain detector is introduced and its application to Teager-Kaiser energy operator is detailed. Simulation results and complexity evaluation are given in Section IV. Finally, conclusions are drawn and outlook is provided.

II. WIRELESS MICROPHONE SENSING

A. The wireless microphone signal model

Wireless microphones operate in the UHF band. Most of them use an analog FM modulation [START_REF] Chen | Spectrum sensing for wireless microphone signals[END_REF]. The signal has a spectral bandwidth B x of 200 kHz. But, most of the signal energy is concentrated in a bandwidth of 100 kHz. The transmit power is a few tens of mW. The coverage area is therefore relatively low, about 500 meters for the most powerful microphones. The microphone signal, x(t), may be modeled as follows:

x(t) = A cos 2πf 0 t + κ f s m τ s(τ )dτ , (1)
where f 0 is the carrier frequency, κ f the frequency deviation of the FM modulation, and s(t) the modulating signal having an amplitude s m . The signal x(t) has a power σ 2 x equals to A 2 /2. Let y(t) be the microphone signal received by the opportunist receiver and n(t) an Additive White Gaussian Noise (AWGN) with a zero mean and a variance σ 2 n : y(t) = x(t) + n(t).

(2)

x(t) and n(t) being independant, the Signal to Noise Ratio (SN R) received by the opportunistic user is:

SN R = σ 2 x σ 2 n = A 2 2σ 2 n . (3)

B. The cognitive radio system description

The CR detector has to detect signals in the presence of noise. The problem can be described by the following hypothesis:

H 0 : y(t) = n(t) H 1 : y(t) = x(t) + n(t) , (4)
where H 0 is the null hypothesis for the event "free band" and H 1 is the alternative hypothesis for the event "occupied band". By choosing one of the two assumptions H 0 and H 1 , two kinds of errors can occur: -a false alarm error: it corresponds to the case when the hypothesis H 1 is chosen while the band is free. The false alarm probability is denoted by p F A , -a non-detection error: it occurs when the hypothesis H 0 is chosen while the band is occupied. The non-detection probability is denoted by p N D . The performance of the detectors is evaluated by computing the detection probability p D (p D = 1 -p N D ) for a given p F A . The p F A imposes the value of the absolute detection threshold. The probability, p D , is given for different powers of the received signal in order to determine the minimum SN R that can be detected. A challenge is to achieve high sensitivity at a low cost of complexity. The digital complexity is calculated in number of operations needed to provide a decision. This paper aims at proposing a detector with the best tradeoff between its detection sensitivity and its complexity. In this study, the detection band is a TV UHF channel, its bandwidth B c is equal to 8 MHz in Europe and 6 MHz in the US.

C. State of the art of wireless microphone detection

Unlike the detection of digital TV signals which can use the characteristics of the Orthogonal Frequency Division Multiplexing modulation [START_REF] Bouzegzi | A second order statistics based algorithm for blind recognition of OFDM based systems[END_REF], wireless microphone sensing is difficult due to the few characteristics of its signal. Most of the literature references use a blind detection for the specific case of wireless microphones [3][4][5] [START_REF] Ghosh | Spectrum Sensing Prototype for Sensing ATSC and Wireless Microphone Signals[END_REF]. Energy and autocorrelation detections are detailed as comparison techniques for our study. The energy detector computes a variable which is proportional to the energy of the received signal [START_REF] Urkowitz | Energy detection of unknown deterministic signals[END_REF] [START_REF] Sonnenschein | Radiometric detection of spreadspectrum signals in noise of uncertain power[END_REF]. The test statistic T of this detector is given by:

T = 1 N s Ns-1 k=0 |y[k]| 2 , (5)
where N s the number of samples of the analyzed signal.

The autocorrelation based detection tests the stationarity of the signal by computing the samples autocorrelation function. First, the receiver estimates the autocorrelation function C y [τ ] of the received signal:

C y [τ ] = Ns-1 k=0 y[k]y * [k -τ ], (6)
where * represents complex conjugation. Then, the following test statistic is calculated [START_REF] Shellhammer | Spectrum sensing in IEEE[END_REF]:

T = Ns-1 τ =0 |C y [τ ]| 2 |C y [0]| 2 . ( 7
)
When no signal is present, the two terms should be roughly equal, since the non-central values (τ =0) should be approximately zero. When the wireless microphone signal is present, the signal is not white and the statistical test should increase.

D. Teager-Kaiser detector

Introduced in [START_REF] Gautier | Teager-Kaiser energy detector for narrowband wireless microphone spectrum sensing[END_REF], this method is suitable to detect wireless microphones. Instead of using the conventional energy detector, we propose to use the Teager-Kaiser energy operator to measure the energy activity of a sample. This operator better reflects the energy of FM signal. In 1990 [START_REF] Kaiser | On simple algorithm to calculate the energy of signal[END_REF], Kaiser used the results of Teager and Teager (especially the energy curve needed to produce speech) which showed the non-linear model of the speech. These changes of the characteristics of the speech signal can be modeled as a linear combination of AM-FM signals.

Based on this model, Kaiser has proposed a very simple and fast algorithm [START_REF] Kaiser | On simple algorithm to calculate the energy of signal[END_REF] to estimate the energy, called the Teager-Kaiser energy operator, whenever the restriction related to the bandwidth of the signal (narrowband signal) is respected. One of the first applications of this operator is the detection of FM modulations. The Teager-Kaiser energy operator Ψ extracts directly the energy from the instantaneous signal and is expressed by:

Ψ {x[k]} = (x[k]) 2 -x[k + 1]x * [k -1]. (8)
In the previous part, we have seen that the wireless microphone is a FM modulation of a speech signal, so the Teager-Kaiser detector should be adapted to this kind of signal. From ( 8), the semi-blind detection could be performed by computing the test statistic:

T =E Ψ {y[k]} , (9) =E Ψ {x[k]} + σ 2 n , ( 10
)
where E is the mathematical expectation operator. The main result in [START_REF] Gautier | Teager-Kaiser energy detector for narrowband wireless microphone spectrum sensing[END_REF] is the sensitivity of the Teager-Kaiser detector. Taking into account the non-linear model of the FM signal energy, this operator leads to a 2 dB gain compared to classical energy detection without an important increase of the complexity. Indeed, the number of real multiplications needed to process N s samples is N s +2 for the Teager-Kaiser detector instead of N s for the energy detector.

III. WIDEBAND FREQUENCY DOMAIN DETECTION

The main purpose of this study is to analyse the band in the frequency domain by using a wideband FFT operator. First, the general architecture of the detector is introduced. Its application to Teager-Kaiser detection is then detailed.

A. General architecture

The proposed scheme follows the idea of spectrum pooling introduced in [START_REF] Weiss | Spectrum pooling: an innovative strategy for the enhancement of spectrum efficiency[END_REF]. Proposed as a transmission technique used by an opportunist user, the idea is to match the size of the FFT used in the digital baseband of the transceiver to the possible bandwith of the application. The bandwidth of one subband should then be an integer multiple of the carrier spacing. Two key advantages can be identified for the transceiver architecture [START_REF] Noguet | Advances in Opportunistic Radio Technologies for TVWS[END_REF]. First, the input of the FFT can be fed with zeros at frequencies where incumbent users are present, while vacant subbands are used for transmission by the opportunistic radio. Second, the FFT operation required at the receiver to invert the modulation may also be used to monitor spectral activity of the other users. The proposed architecture is based on this second idea. Fig. 1 shows the general architecture of the proposed detector. It is composed of three main blocs: an FFT operation, a "Frequency Bin Selection" bloc and several "Frequency Analysis blocs". In order to compute a large FFT, the N s samples of the signal are split into N t slots of size N F F T which is the size of the FFT. The outputs of the FFT operation are:

Y (ν) = N F F T -1 k=0 y[k]e -j2π ν N F F T k . (11)
The "Frequency Bin Selection" is a way of scanning the whole band, its parameters are the number of subbands M and the overlapping ratio OR between two adjacent bands. In that case, the number of frequency bins per subbands (denoted by N ) is:

N = N F F T M (1 -OR) . ( 12
)
Theses parameters are set according to the signal. The frequency domain energy detection could be directly computed from the frequency bins Y (ν). It is defined by the following metric:

T i = 1 N Bi+N -1 ν=Bi |Y (ν)| 2 , (13)
B i =N (1-OR)i being the first sample of the subband i (i=0, . . . M -1). Contrary to energy detection, the autocorrelation detection must be computed in the time domain by using an IFFT of size N . From the frequency bins Y (ν), the autocorrelation function per subband C i (i=0, . . . M -1) is computed by:

C i [τ ] = Bi+N -1 ν=Bi Y [ν]Y * [ν]e j2π ν N τ . (14)
After the frequency analysis bloc, the decision is made by averaging N t FFT frames.

B. Application to Teager-Kaiser detection

In the particular case of detection of wireless microphones as described in Section II, we propose an architecture of the Frequency Analysis bloc based on the Teager-Kaiser energy operator. Using the autocorrelation function described in ( 6), the Teager-Kaiser detector could be written in the time domain by:

T = E Ψ {y[k]} = E y 2 [k] -E y[k + 1]y * [k -1] , ( 15
) = C y [0] -C y [2]. (16)
C y [0] is the energy of the signal and C y [START_REF]Official announcement of FCC[END_REF] is the tap 2 of the autocorrelation function. So to compute the Teager-Kaiser energy operator from the FFT of the signal, the two following terms should be calculated:

C y [0] = N -1 ν=0 Y [ν]Y * [ν], (17)
and

C y [2] = N -1 ν=0 Y [ν]Y * [ν]e j2π ν N 2 . ( 18
)
Fig. 2 describes the structure of the frequency analysis bloc in the case of the Teager-Kaiser detection.

C. Architecture settings

In the following, the wireless microphone study case is used to analyze the proposed architecture. Based on the European UHF band specifications, the bandwidth is 8 MHz. We consider that a large FFT of 4096 bins is used as a fixed parameter of the system. The baseband microphone signal is generated following the model of ( 1) with these parameters: the frequency bandwidth of s(t) is 20 kHz, the FM deviation κ f is 3, the carrier frequency f 0 is equal to 100 kHz and the noise power σ 2 n is set to 0 dBm. Under these conditions, the transmitted microphone signal x(t) has a frequency bandwidth of around 100 kHz.

The first characteristic that must be specified is the number of input samples that can be analyzed by this architecture. Indeed, according to the test statistic, a bias could appear between the time and the frequency domain computations of the metric. For example, a FFT size of N F F T can be used to compute the frequency domain energy of N F F T input samples without any bias. However, for the autocorrelation and Teager-Kaiser metrics, a bias may occur depending the size of the inputs samples. Fig. 3 gives the total introduced bias for both metrics as a function of the ratio between the input sample size and the FFT size. Results show that, for the autocorrelation function, a ratio of more than 0.5 induces a significant bias. The bias increases then more rapidely as the ratio is increased. For the Teager-Kaiser energy operator, the bias starts to become significant only for a ratio closed to 1. This bias is lower than the autocorrelation for the same ratio of input size over FFT size. This analysis outlines that the comparison of the different "Frequency Analysis bloc" is not obvious. Indeed, using a fixed FFT size, the same number of samples could not be analyzed without any bias. Thus, in the following, only the energy detection and the Teager-Kaiser detection will be computed in the frequency domain. For the Teager-Kaiser detection, 2 samples (insignificant compared to the FFT size of 4096) won't be used in order to have no bias.

In a second step, the number of subbands M has to be set. Obtained by simulation, the sensitivity of the two frequency domain detectors is plotted versus the number of subbands. The sensitivity is the minimum SN R that the detector could sense with p D = 95% and p F A = 1%. The results are shown in Fig. 4 for an OR of 20% and N t =1. Simulation results show that the sensitivity increases with the number of subchannel. For large number of subchannel, the SN R decreases with a lower step and seem to reach a SN R floor. This occurs when the subband bandwidth becomes smaller that the signal bandwidth. In the following, M =16 will be used.

IV. PERFORMANCE ANALYSIS

A. Detection sensitivity

The main feature of a detector is its sensitivity level expressed in SN R. Fig. 5 shows the performance of the two frequency domain algorithms. Performance is compared with the time domain algorithms. The detection probability p D is computed versus the SN R of the wireless microphone received signal. The architecture settings are M =16, OR=20%, N t =1 and N F F T =4096 and the p F A is 1%. For a target p D of 95%, the frequency domain detectors enable the detection of signals with SN R of -14 dB and -15.5 dB for the energy and the Teager-Kaiser metrics respectively. In the time domain, the sensitivity levels are -7.5 dB for the energy detector, -9.5 dB for the Teager-Kaiser detector and -11.5 dB for the autocorrelation detector. The simulation results show that, to narrowband the analyzed bandwidth leads to a 6 dB improvement of the time domain energy detection. A further 1.5 dB gain is achieved using the Teager-Kaiser energy operator. Both frequency domain detectors outperform the autocorrelation detector with a gain of 4 dB for the Teager-Kaiser detector and a gain of 2.5 dB for the energy detector.

B. Complexity evaluation

After the sensitivity evaluation, the complexity of the frequency domain algorithms should be evaluated and compared to their time domain counterpart. A simple way to estimate complexity is to determine the total number of real multiplications for each algorithm. Assuming split radix-2 implementation of the FFT, the number of real multiplications to process N F F T samples are:

• for the frequency domain energy detector: • for the frequency domain Teager-Kaiser detector:

N F F T * log 2 (N F F T ) + N F F T , (19)
(N F F T + 2) * log 2 (N F F T + 2) + 2 * (N F F T + 2). ( 20
)
The complexity is compared with the one of a filterbank based architecture as introduced in [START_REF] Gautier | Teager-Kaiser energy detector for narrowband wireless microphone spectrum sensing[END_REF]. In order to split the frequency band into subbands, an alternative to FFT is the use of a filterbank in order to channelized the wideband. To understand the requirements on the filter used by the filterbank, we derive a rejection requirement A (in dB) on the filter and then extrapolate the number of taps required for a finite impulse response implementation of the filter. The filter length L f ilter is estimated using an equiripple lowpass implementation. Using a low complexity dyadic implementation of the filterbank [START_REF] Akansu | Multiresolution Signal Decomposition: Transforms, Subbands and Wavelets[END_REF], the complexity of this "channelized structure" is:

L f ilter * M * log 2 (M ) * N F F T /M, (21)
to process N F F T samples. Then, N F F T and 2 * N F F T real multiplications must be added to compute the channelized energy detection and the channelized Teager-Kaiser detection respectively. Fig. 6 shows the complexity of each option as a function of the FFT size. For the channelized structure, two rejection requirements are used: A=10 dB and A=30 dB. It is obvious that the time domain algorithms have lower complexity when no channelization is done. However, the results underline that the proposed FFT based detector has a lower complexity compared with the filterbank based architecture. Furthermore, this study doesn't take into account the fact that the FFT operation could be performed with the same FFT operation as the demodulation scheme. In that case, the FFT comes for free and the frequency domain detectors have nearly the same complexity as their time domain counterpart.

V. CONCLUSIONS

In this paper, the detection of narrowband wireless microphone signals has been addressed. The motivation of this study was to identify a wireless microphone detector that could be implemented with an acceptable complexity. This paper focused on the complexity-performance tradeoffs of the algorithms and two main conlusions can be drawn. First, wideband signal may be monitored using FFT operation. Since the FFT implementation may be shared between communication and sensing schemes, the approach leads to a low complexity detector. In the wireless microphone study case, the proposed architecture brings a gain of 6 dB gain compared to time domain detectors. Second, the computation of the Teager-Kaiser energy operator may favorably be implemented in the frequency domain and brings a gain of 1.5 dB compared to frequency domain energy detection. Future work will be dedicated to the study of a detection threshold calibration. Then, a hardware implementation of the proposed algorithms and experimental measurements should be realized.

Fig. 1 .

 1 Fig. 1. Architecture of the wideband frequency domain detector.

  For example, if a 40 MHz band is analyzed and a 8 MHz channel is sensed, M =5 subbands and no overlapping, OR=0, are needed. If wireless microphone signals are sensed in a 8 MHz band, overlapping is needed because of the random frequency position of the signal and the bandwidth of each subband must be set to 200 kHz. The "Frequency Analysis bloc" provides a test statistic from the frequency bins selected by the previous bloc. One bloc uses N samples to compute the metric. Several test statistics can be computed in the frequency domain. Starting with the detectors introduced in section II-C as the state of the art for wireless microphone detection, two test statistics are first introduced in the frequency domain.

Fig. 2 .

 2 Fig. 2. Frequency analysis bloc for Teager-Kaiser detection.

Fig. 3 .

 3 Fig. 3. Bias versus ratio between Input size and FFT size.

Fig. 4 .

 4 Fig. 4. Sensitivity in SN R versus number of subchannel.

Fig. 5 .

 5 Fig. 5. Detection probability versus SN R.

Fig. 6 .

 6 Fig. 6. Complexity in number of real multiplications versus FFT size.

ACKNOWLEDGMENT

The research leading to these results was derived from the European Community's Seventh Framework Program (FP7) under Grant Agreement number 248454 (QoSMOS).


Manuel Bouillon
email: manuel.bouillon@irisa.fr

Eric Anquetil
email: eric.anquetil@irisa.fr

Abdullah Almaksour
email: abdullah.almaksour@irisa.fr

Decremental Learning of Evolving Fuzzy Inference Systems Using a Sliding Window

Keywords: Incremental Learning, Decremental Learning, Evolving Fuzzy Inference System, Recursive Least Squares, Concept Drifts, On-line Gesture Recognition

This paper tackles the problem of decremental learning of an evolving classification system. We study the use of decremental learning to improve performance of evolving recognizers in non-stationary scenarios. Our on-line recognizer is based on an evolving fuzzy inference system. In this paper, we propose a new strategy to introduce decremental learning, with the use of a sliding window, in the optimization of fuzzy rules conclusions. This approach is based on a downdating technique of least squares solutions for unlearning old data. This technique is evaluated on handwritten gesture recognition tasks. In particular, it is shown that this downdating techniques allow to adapt to concept drifts and that we face a precision reactiveness trade-off. It is also demonstrated that decremental learning is necessary to maintain the system learning capacity over time, making decremental learning essential for the lifetime use of an evolving classification system.

I. INTRODUCTION

Evolving classification systems have appeared in the last decade to meet the need for recognizers that work in changing environments. They use incremental learning to adapt to the data flow and to cope with class adding (or removal) at runtime. This paper focuses on the decremental learning of such an evolving recognizer, based on a fuzzy inference system. The aim of decremental learning is twofold. First, to maintain the system learning capacity, and second, to forget obsolete data and focus the learning process on current environment.

The target application of this work is the use of on-line handwriting gesture classifiers to facilitate interactions with computers using pen-based interfaces like tablet computers, smartphones, whiteboards, etc. Gestures can be drawn differently from one user to another, and users may want to add or remove gestures, as long as they use the application. Moreover, users would often change progressively the manner by which they draw gestures. Novice users start drawing carefully and slowly their gestures, while they do them in a more fluid and rapid manner as they become expert. To cope with these requirements, forgetting capacity must be used to increase system reactivity and performance in such dynamic environment.

We extend in this paper our evolving classification system Evolve [START_REF] Almaksour | Improving premise structure in evolving takagi-sugeno neuro-fuzzy classifiers[END_REF] by integrating a decremental learning technique. A new decremental strategy, relying on a sliding window of data samples, is proposed. This window is used to unlearn old data -so-called downdating of least squares solutions.

We present briefly the architecture of Evolve and its incremental learning algorithm in Section II. Our proposed decremental learning method is detailed in Section III. This new approach is then evaluated on some handwritten gesture recognition tasks in section IV. Section V concludes and discusses future work.

II. SYSTEM ARCHITECTURE

We focus here on Fuzzy Inference Systems (FIS) [START_REF] Lughofer | Evolving fuzzy models: incremental learning, interpretability, and stability issues, applications[END_REF], which have demonstrated their good performance for classification for quite some time. Moreover, they can easily be trained incrementally (in real time) and have a good behavior when new classes are added. [START_REF] Angelov | Evolving fuzzy-rule-based classifiers from data streams[END_REF] and [START_REF] Almaksour | Improving premise structure in evolving takagi-sugeno neuro-fuzzy classifiers[END_REF] are recent examples of evolving FIS used for on-line classification.

Our evolving system Evolve [START_REF] Almaksour | Improving premise structure in evolving takagi-sugeno neuro-fuzzy classifiers[END_REF] is a first-order Takagi-Sugeno (TS) fuzzy inference system. It consists of a set of fuzzy inference rules. Rule premises are the fuzzy membership to the rule prototype, which are clusters in the input space. Rule conclusions are fuzzy membership to all classes, as a linear function of the input, that are combined to produce the system output.

In Evolve, both the rule premises and conclusions are incrementally adapted. Rule prototypes are statistically updated to model the runtime data. Rule conclusions parameters are optimized on the data flow, using Recursive Least Squares (RLS) algorithm [START_REF] Haykin | Adaptive Filter Theory[END_REF]. New rules, with their associated prototypes and conclusions, are created by an incremental clustering method when needed.

As rule conclusions induce decision boundaries, we decided to introduce decremental learning in the conclusions optimization process, namely the RLS algorithm.

III. DECREMENTAL LEARNING OF EVOLVING FUZZY

INFERENCE SYSTEMS This paper focuses on decremental learning for the lifelong use and learning of an evolving classification system optimized with the Recursive Least Squares (RLS) algorithm.

The interest of integrating forgetting in a learning system is twofold. First, decremental learning is necessary to limit the weight of past data, and thus to maintain the learning capacity of the system. If every data sample has the same relative weight, the system will tend to become set with time. Second, decremental learning allows the system to follow any change -concept drift [START_REF] Widmer | Learning in the presence of concept drift and hidden contexts[END_REF] -of the data flow by unlearning obsolete data samples.

Several decremental systems already exist in different contexts. The use of a forgetting factor in RLS algorithm [START_REF] Kulhavy | On a general concept of forgetting[END_REF] is well known in the literature dedicated to control of complex systems. It comes to using the RLS algorithm on an (exponentially weighted) sliding window. However, this algorithm behaves poorly when systems are not uniformly exited -it is known as the covariance "wind-up" problem [START_REF] Liavas | On the numerical stability and accuracy of the conventional recursive least squares algorithm[END_REF] -which is the case in classification problems.

We present here a new decremental learning strategy applied to FIS optimization with the RLS algorithm. We propose a new decremental learning approach based on downdating the least squares solutions with a sliding window. We use a real time algorithm to unlearn old data that leave the window: the de-recursive least squares algorithm.

A. Principle

The principle of this approach is simple, we maintain a sliding window over the latest data, and we optimize the rules conclusions only on this window of data.

As we can't afford to rebuild rules conclusions at the arrival of each new data, they are updated using the recursive least squares (RLS) algorithm [START_REF] Almaksour | Improving premise structure in evolving takagi-sugeno neuro-fuzzy classifiers[END_REF]. In the same way, we need to downdate the least squares solutions at departure of each old data from the window, without complete rebuilding. To do so, we propose to use the de-recursive least squares algorithm in order to recursively unlearn old data.

B. De-Recursive Least Squares (DRLS) Algorithm

Let θ (i) s→t ∈ R n×c be the i th rule conclusion, optimized on data samples from (x s , y s ) to (x t , y t ) (with x k ∈ R n the input vectors, and y k ∈ R c the binary output vectors). We propose to downdate least squares solutions with the following DRLS formulas to unlearn "old" data sample (x s , y s ).

θ (i) (s+1)→t = θ (i) s→t -P (i) (s+1)→t x s β (i) s (y s -x T s θ (i) s→t ) (1)
Where the covariance matrices P (i) are updated as follows.

P (i) (s+1)→t = P (i) s→t + P (i) s→t x s x T s P (i) s→t -1 β (i) s + x T s P (i) s→t x s (2)
These DRLS formulas can easily be proven in a similar way to the standard RLS formulas.

C. Discussion of the Window Size

The sensitive issue of this approach is the size of the sliding window. Indeed, performance is directly linked to the window length. In steady environment, the longer is the window, the lower is the error rate (until a minimum rate). However, a too long window reduces system reactiveness and thus deteriorates performance in changing environment.

We focus on having a large enough window to avoid performance reduction in stationary environment. Such a large window is generally sufficient to limit old data weight and to adapt to concept drifts in non-stationary environment. Even if a fixed length sliding window is not optimal, it provides a good behavior as it will be demonstrated experimentally.

The window minimum length is a problem dependent variable that must be empirically determined. We recommend to dynamically link the window size to the number of classes (noted c) to avoid performance scale-down when new classes are added to the system.

IV. EXPERIMENTAL RESULTS

Starting from the state-of-the-art recognizer Evolve [START_REF] Almaksour | Improving premise structure in evolving takagi-sugeno neuro-fuzzy classifiers[END_REF], we implemented our new approach: Evolve D -with Downdating. All systems use the Heterogeneous Baseline Feature set (HBF49) [START_REF] Delaye | Hbf49 feature set: A first unified baseline for online symbol recognition[END_REF].

In this section, we first present the incremental evaluation protocol, and the datasets, used for testing. Then, we measure system inertia to novelty to asses the need for decremental learning. Next, we study the impact of the window length in both stationary and non-stationary environment.

A. Incremental Evaluation Protocol

To evaluate our systems in a realistic way, we used an incremental evaluation protocol called predictive sequential -or prequential [START_REF] Gama | Issues in evaluation of stream learning algorithms[END_REF] -with a sliding window to converge to the holdout error.

As an incremental system first tries to recognize a data sample, and then learn from it once it has the true label, we evaluated our systems in a similar way. Each data sample is first used as test sample, and then as learning sample. Error rates are then computed between every test points.

B. Datasets

As this work is applied to on-line handwritten gesture recognition, we evaluated our new approach on handwritten gestures databases: ILGDB 1 [START_REF] Renau-Ferrer | The IL-GDB database of realistic pen-based gestural commands[END_REF] and IRONOFF-digit [START_REF] Viard-Gaudin | The ireste on/off (ironoff) dual handwriting database[END_REF].

1) ILGDB: This database contains handwritten gestures that have been collected in an immersive environment. It is composed of 6629 mono-stroke gestures, belonging to 21 classes, which were written by 38 writers. This database is very interesting for several reasons.

First, gestures are ordered chronologically in their drawing order which allows us to see changes in writer style 1 Freely available at http://www.irisa.fr/intuidoc/ILGDB.html The only drawback of this database is the low number of gesture per writer (less than 180) and per class for each writer (from 5 to 17 gestures). To make up for this inconvenience, we also used another database, containing more data, presented below.

2) IRONOFF-digits: The interest of this database is that it offers more gestures (in writer independent mode), but like most classic benchmark databases, IRONOFF-digit is not ordered (no chronological evolution of the data with time).

C. System Inertia to Novelty

We test the reactiveness of the systems, and the evolution of the inertia of their model with time. We train the different systems with seven classes during a varying period time and then introduce three other classes. We measure the time needed by the different systems to learn those new classes. Results averaged over 20 different data orders are shown in figure 1.

If all systems behave similarly after 70 training samples, results are different after 700 or 2100 samples. Looking at the results of our reference system without forgetting Evolve, we can clearly see that its reactiveness decreases with time, and that the model tends to become set. After 2100 training samples, the system needs 400 more samples to resume to an error rate under 10%, which makes it partially unusable for quite some time. On the contrary, systems using decremental learning need the same time to learn and adapt their model to the novelty, whenever it is introduced. Their reactiveness is independent of their age.

This test scenario shows the necessity of integrating forgetting into an incremental learning system operating in a non-stationary environment.

D. System Performance in Stationary Environment

We test our different systems on a scenario simulating soft concept drifts. For this purpose, we used the ILGDB 11 writers of group 3 (whose gestures for each class are identical) in a row. We computed the error rate for each writer on the two last phases, approximately the 56 last samples of the 173 of each writer, to measure system performance after the concept drift. Mean results over 20 writer orders are plotted figure in 3. Writers, and their drawing style, changes but the base gestures of each class stay unchanged. This test scenario is thus nearly stationary and doesn't require the use of decremental learning. Our reference system Evolve achieves quite good results here with an average error rate of 5.25%.

Evolve D performance is limited by the restrained number of data it is learning from. With a sliding window of 20 and 10 times the number of classes, Evolve D reaches an error rate of 5.38% and 5.61% respectively, which is as good as Evolve. A smaller window length reduces performance.

E. System Performance in Non-Stationary Environment

We also test our three systems on a scenario simulating abrupt concept drifts. For this purpose, we used the ILGDB 21 writers of group 1 (whose gestures for each class are different) in a row. We computed the error rate for each writer on the two last phases, approximately the 56 last samples of the 173 of each writer, to measure performance after the concept. Mean results over 20 writer orders are plotted in figure 4. The results on this scenario are sharply contrasted. Without forgetting, Evolve is not able to adapt its model to the changes in the data and ends up with an error rate of 50%. The use of decremental learning allows Evolve D to follow the abrupt concept drifts. We can note that a too long window size (20 times the number of classes) deteriorates performance in this changing environment.

F. Results Discussion

Results on those non-stationary scenarios are quite clearcut. They show the necessity of using decremental learning to maintain the reactiveness of the system over time. Without downdating, the system model tends to become set, and take ages to learn some novelty. Decremental learning is necessary to maintain system learning capacity.

In the same way, we showed that decremental learning is essential to face (abrupt) concept drifts. Without downdating, the system model becomes more and more complex and performance slowly but surely collapses. Decremental learning allows to discard obsolete data and thus enables the system to focus on current system environment.

We also highlighted the trade-off in the choice of the window length. A short window insures high reactiveness, but deteriorates performance in steady environment. Whereas a long window enables good performance in stationary environment but deteriorates performance in changing environment. For ILGDB, a window length of 10 times the number of classes (210 samples) is a fine compromise that gives good performance in both cases.

V. CONCLUSION

In this paper, we investigated the decremental learning of an evolving fuzzy inference system used for classification. We proposed a new approach to integrate downdating in the optimization of rules conclusion. Our approach is based on a sliding window and uses de-recursive least squares formulas to decrementally downdate least squares solutions and unlearn "old" data samples from rule conclusions.

This work is applied on handwritten gesture recognition. In particular, we demonstrated that decremental learning is essential to maintain the system reactiveness over time and to adapt to concept drifts. We showed that our approach performs well in changing environment, without deteriorating performance in stationary environment.

Future Work: These approaches could be improved to lessen the accuracy reactiveness trade-off. The sliding window length could be adapted on-line to fit the data flow. The window length could be increased to improve performance during when the data flow is stationary, and reduced when concept drifts are detected.

Figure 1 :

 1 Figure 1: Scenario "inertia and reactiveness" -IRONOFF-digit database -Introducing novelty after 70, 700 and 2100 training samples

Figure 2 :

 2 Figure 2: Gesture samples from ILGDB group 1 (free gestures)

Figure 3 :

 3 Figure 3: Scenario "soft concept drifts" -Using ILGDB 11 writers of group 3 (fixed gestures) in a row

Figure 4 :

 4 Figure 4: Scenario "abrupt concept drifts" -Using ILGDB 21 writers of group 1 (free gestures) in a row