Olivier Cappé Ltci

Telecom

Online EM Algorithm for Hidden Markov Models

Keywords: Hidden Markov Models, Expectation-Maximization Algorithm, Online Estimation, Recursive Estimation, Stochastic Approximation, Smoothing

Online (also called "recursive" or "adaptive") estimation of fixed model parameters in hidden Markov models is a topic of much interest in times series modelling. In this work, we propose an online parameter estimation algorithm that combines two key ideas. The first one, which is deeply rooted in the Expectation-Maximization (EM) methodology consists in reparameterizing the problem using complete-data sufficient statistics. The second ingredient consists in exploiting a purely recursive form of smoothing in HMMs based on an auxiliary recursion. Although the proposed online EM algorithm resembles a classical stochastic approximation (or Robbins-Monro) algorithm, it is sufficiently different to resist conventional analysis of convergence. We thus provide limited results which identify the potential limiting points of the recursion as well as the large-sample behavior of the quantities involved in the algorithm. The performance of the proposed algorithm is numerically evaluated through simulations in the case of a noisily observed Markov chain. In this case, the algorithm reaches estimation results that are comparable to that of the maximum likelihood estimator for large sample sizes.

Introduction

Hidden Markov modelling is a key concept of statistical time series analysis, which has had a wide-ranging practical impact over the latest forty years. Hidden Markov models (HMMs) in their classical form (i.e., when the state variable is finite-valued) are sufficiently simple to give rise to efficient inference procedures while allowing for useful modelling of various practical situations. Ever since the pioneering contributions of [START_REF] Baum | An inequality with applications to statistical estimation for probalistic functions of Markov processes and to a model for ecology[END_REF], [START_REF] Baum | A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains[END_REF], the EM (Expectation-Maximization) algorithm has been the method of choice for parameter inference in HMMs. The EM algorithm is a dedicated numerical optimization routine which aims at maximizing the (log) likelihood of a batch of observations. It tends to be preferred to its alternatives due to its robustness and ease of implementation.

This contribution is devoted to online parameter estimation for HMMs, in which the available observations are only scanned once and never stored, allowing for a continuous adaptation of the parameters along a potentially infinite data stream. In the case of HMMs, online parameter estimation is a challenging task due to the non-trivial dependence between the observations. The EM-inspired methods proposed so far have been either based on finite-memory approximations of the required smoothing computations [START_REF] Krishnamurthy | On-line estimation of hidden Markov model parameters based on the kullback-leibler information measure[END_REF] or on finite-memory approximations of the data log-likelihood itself [START_REF] Rydén | On recursive estimation for hidden Markov models[END_REF]). An alternative consists in using gradient-based methods [START_REF] Le Gland | Recursive estimation in HMMs[END_REF] which do not directly follow the principles of 2 Online EM Algorithm for HMMs

Model and Notations

It is assumed that the state and observation sequences, (X t , Y t ) t∈Z are generated under a stationary Hidden Markov model with unknown parameter θ ⋆ , where (X t ) takes its values in some finite set X . The notations ℓ θ⋆ , p θ⋆ , P θ⋆ and E θ⋆ refer to, respectively, the likelihood, the joint density of the states and observations, the probability, and, the expectation under the model parameterized by θ ⋆ . In practice θ ⋆ is unknown and one only has access to the observation sub-sequence (Y t ) t≥0 , where, by convention, the initial observation time is taken to be 0. In this context, the initial pdf ν of X 0 is arbitrary and ℓ ν,θ , p ν,θ , P ν,θ and E ν,θ are used to denote, respectively, the likelihood, the joint density of the states and observations, the probability, and, the expectation under the non-stationary model that has initial pdf ν and θ as parameter. Hence, P θ refers to the probability under θ of the stationary HMM process (X t , Y t ) t∈Z while P ν,θ denotes the probability of the nonstationary process (X t , Y t ) t∈N which is started with initial distribution ν for X 0 . Note that ν itself is not considered as a model parameter as it cannot be estimated consistently from a single trajectory (see Chapters 10 and 12 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF] for further discussion of this issue). The state transition matrix and state conditional pdf (probability density function of Y t given X t ) that characterize the HMM are denoted, respectively, by q θ (x, x ′ ) and g θ (x, y).

To make the EM recursion explicit, a key requirement is that the model belongs to an exponential family. In the following, we will thus make the following assumptions.

Assumption 1.

(i) Exponential Family p θ (x t , y t |x t-1 ) = h(x t , y t ) exp ( ψ(θ), s(x t-1 , x t , y t ) -A(θ)) ,

(1)

where • denotes the scalar product, s : (x ′ , x, y) ∈ X 2 × Y → s(x ′ , x, y) ∈ S is the vector of complete-data sufficient statistics, ψ(•) is the (non-necessarily invertible) function that maps θ to the natural parameterization and A(•) is the log-partition function.

(ii) Explicit M-Step For all S ∈ S, the complete-data maximum likelihood equation

∇ θ ψ(θ)S -∇ θ A(θ) = 0 ,
where ∇ θ denotes the gradient, has a unique solution denoted by θ(S).

Assumption 1-(ii) states that the function θ : S ∈ S → θ(S) ∈ Θ that returns the completedata maximum likelihood estimator corresponding to any feasible value of the sufficient statistics is available in closed-form. Note that the form used in (1) is in fact slightly more general than the HMM case. Indeed, if both q θ and g θ belong to exponential families:

q θ (x ′ , x) = h q (x ′ , x) exp ψ q (θ), s q (x ′ , x) -A q (θ) , g θ (x, y) = h g (x, y) exp ( ψ g (θ), s g (x, y) -A g (θ)) , (2)
we then have A(θ) = A q (θ) + A g (θ),

ψ(θ) = ψ q (θ) ψ g (θ)
and s(x ′ , x, y) = s q (x ′ , x) s g (x, y) .

We will not make use of this specific structure when describing the algorithm and we thus stick to the concise representation of (1). Sections and 3.2 and 5 provide more details on the nature of the function s in specific examples of HMMs. Under Assumption 1, the k-th iteration of the usual EM algorithm applied to observations Y 0:n = (Y 0 , . . . , Y n ) takes the following familiar form:

E-Step Compute S k+1 = 1 n E ν,θ k n t=1 s(X t-1 , X t , Y t ) Y 0:n . (3)

M-

Step Update the parameter estimate to θ k+1 = θ(S k+1 ).

To avoid unnecessary notational complexity, we have omitted in (3) the initial term 1 n log p ν,θ (x 0 , y 0 ) from the normalized complete-data log-likelihood. The influence of this term is vanishing with n and it is not necessary to take it into account for online estimation.

Recursive Form of Smoothing

We briefly recall here a key ingredient of the proposed algorithm which makes it possible to compute recursively the normalized sum

1 n E ν,θ [ n t=1 s(X t-1 , X t , Y t )| Y 0:n ]
by use of an auxiliary recursion. Curiously, this idea which dates back to, at least, [START_REF] Zeitouni | Exact filters for the estimation of the number of transitions of finite-state continuous-time Markov processes[END_REF] and has been extensively studied by [START_REF] Elliott | Hidden Markov Models: Estimation and Control[END_REF] remains largely under-exploited (see discussion in Chapter 4 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF].

In addition to the usual filter,

φ n,ν,θ (x) = P ν,θ ( X n = x| Y 0:n ) , (4)
define the following intermediate quantity

ρ n,ν,θ (x) = 1 n E ν,θ n t=1 s(X t-1 , X t , Y t ) Y 0:n , X n = x . (5)
Obviously, these two quantities allow the computation of the sum of interest as

x∈X φ n,ν,θ (x)ρ n,ν,θ (x) = 1 n E ν,θ n t=1 s(X t-1 , X t , Y t ) Y 0:n .
The appeal of this decomposition is that φ n,ν,θ and ρ n,ν,θ can be updated recursively according to the following proposition.

Proposition 1. Initialization For x ∈ X , set

φ 0,ν,θ (x) = ν(x)g θ (x, Y 0 ) x ′ ∈X ν(x ′ )g θ (x ′ , Y 0 ) , ρ 0,ν,θ (x) = 0 .
Recursion For n ≥ 0 and x ∈ X , it holds that

φ n+1,ν,θ (x) = x ′ ∈X φ n,ν,θ (x ′ )q θ (x ′ , x)g θ (x, Y n+1 ) x ′ ,x ′′ ∈X 2 φ n,ν,θ (x ′ )q θ (x ′ , x ′′ )g θ (x ′′ , Y n+1 ) (6) ρ n+1,ν,θ (x) = x ′ ∈X 1 n + 1 s(x ′ , x, Y n+1 ) + 1 - 1 n + 1 ρ n,ν,θ (x ′ ) φ n,ν,θ (x ′ )q θ (x ′ , x) x ′′ ∈X φ n,ν,θ (x ′′ )q θ (x ′′ , x) . (7)
In Proposition 1 above, the rightmost term in (7),

r n+1,ν,θ (x ′ |x) = φ n,ν,θ (x ′ )q θ (x ′ , x) x ′′ ∈X φ n,ν,θ (x ′′ )q θ (x ′′ , x) , corresponds to the backward retrospective probability P ν,θ (X n = x ′ |X n+1 = x, Y 0:n ),
P ν,θ (X t = x t , X t+1 = x t+1 |X n+1 = x n+1 , Y 0:n+1 ) = xn∈X P ν,θ (X t = x t , X t+1 = x t+1 |X n = x n , Y 0:n )P ν,θ (X n = x n |X n+1 = x n+1 , Y 0:n ) ,
for all indices 0 ≤ t ≤ n -1 which implies the claimed result by summation on t.

Online EM Algorithm

Proposition 1 constitutes a recursive rewriting of the computation required to carry out the Estep in the batch EM algorithm. By analogy with the case of independent observations studied in [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF], the proposed online EM algorithm for HMMs takes the following form.

Algorithm 1. Chose a decreasing sequence (γ n ) n≥1 of step-sizes, which satisfy the usual stochastic approximation requirement that n≥1 γ n = ∞ and n≥1 γ 2 n < ∞. Also select a parameter initialization θ0 and a minimal number of observations n min required before performing the first parameter update.

Initialization Compute, for x ∈ X , φ0 (x) = ν(x)g θ0 (x, Y 0 ) x ′ ∈X ν(x ′ )g θ0 (x ′ , Y 0 ) , ρ0 (x) = 0 . Recursion For n ≥ 0, Compute, for x ∈ X , φn+1 (x) = x ′ ∈X φn (x ′ )q θn (x ′ , x)g θn (x, Y n+1 ) x ′ ,x ′′ ∈X 2 φn (x ′ )q θn (x ′ , x ′′ )g θn (x ′′ , Y n+1 ) , (8) ρn+1
(x) = x ′ ∈X γ n+1 s(x ′ , x, Y n+1 ) + (1 -γ n+1 )ρ n (x ′ ) φn (x ′ )q θn (x ′ , x)
x ′′ ∈X φn (x ′′ )q θn (x ′′ , x)

.

If n ≥ n min , update the parameter according to

θn+1 = θ x∈X ρn+1 (x) φn+1 (x) , (10)
otherwise, set θn+1 = θn .

In Algorithm 1, the role of n min is only to guarantee that the M-step update is numerically wellbehaved as it is well-known that in most models the maximum likelihood estimation equation is degenerate for very small numbers of observations and, hence, the function θ may not be properly defined. For this purpose, a very small value of n min is usually sufficient (for instance, n min = 20 is used in the simulations of Section 5.3).

Discussion

We first discuss connections of the proposed algorithm with earlier works before discussing its numerical complexity in more details.

Quite obviously, Algorithm 1 is intended to generalize the online EM algorithm of [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF] to the case of dependent observations. In the case of independent observations, (8) and (9) reduces to a simpler recursion of the form:

Ŝn+1 = γ n+1 E θ( Ŝn) [ s(X n+1 , Y n+1 )| Y n+1 ] + (1 -γ n+1 ) Ŝn , (11)
which can be analyzed using the arguments developed for the analysis of stochastic approximation (or Robbins-Monro) schemes under Markovian random perturbations. The HMM framework however implies several key differences. First, in HMMs it is necessary to maintain an approximate filter φn through (8) which, hopefully, becomes an acceptable proxy for φ n,ν, θn when approaching convergence. More importantly, as it is no more possible to compute the conditional expectation of the complete-data sufficient statistics online, we require the auxiliary recursion of (9). Although directly inspired by the exact recursive smoothing formula of (7), this auxiliary update cannot be put in the usual stochastic approximation form. And thus, in contrast to the online EM algorithm of [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF], Algorithm 1 cannot be analyzed using off-the-shelf mathematical arguments. Both the need to maintain an estimate of the filter and the presence of the backward retrospective probability in (9) constitute significant departures from the usual stochastic approximation framework. Thus a complete analysis of the convergence of Algorithm 1 is a challenging mathematical problem and we provide in Section 4 below important -though limited-arguments in that direction. Algorithm 1 is related to the EM-based approach suggested in Eq. ( 51) of [START_REF] Kantas | An overview of sequential Monte Carlo methods for parameter estimation in general state-space models[END_REF] which however lacks the idea of a recursive implementation of smoothing (in this work, the authors propose to use a particle filtering method to approximate the smoothing functional). Chapters 4 and 5 of [START_REF] Ford | Adaptive hidden Markov model estimation and applications[END_REF] and [START_REF] Elliott | On-line almost-sure parameter estimation for partially observed discrete-time linear systems with known noise characteristics[END_REF] consider the use of Proposition 1 for online estimation of the parameter of a Markov chain observed in additive noise -an example that we will discuss in more detail in Section 5-and for a linear state-space model with known noise characteristics. Ford (1998) however does not acknowledge the generality of the approach and the role of complete-data sufficient statistics and is constrained to the choice of γ n = n -1 , which in simulations gives poor performances compared to alternative choices of step-sizes (see Section 5 below). To the best of our knowledge, the first instance of a particular case of Algorithm 1 appears in the work of [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] who considered the specific case of finite-valued observations. [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] however entirely rediscovered Proposition 1 under an equivalent form discussed below and thus failed to identify some of the general principles underpinning Algorithm 1. [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] This section is devoted to a more detailed analysis of the difference between Algorithm 1 applied to the case of finite-valued observations and the actual proposal of [START_REF] Mongillo | Online learning with hidden Markov models[END_REF]. This difference is not significant from a practical point of view but is important for the understanding of the behavior of the algorithm. [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] considered the particular case of finite valued HMMs in which the observations (Y t ) t≥1 also take their values in a finite set Y. In such a situation, it is easily checked that for any parameterization of the model, the complete-data sufficient statistics may be chosen as s

Comparison with the Algorithm of

(X t-1 , X t , Y t ) = (½{X t-1 = i, X t = j, Y t = k}) (i,j,k)∈X 2 ×Y .
The recursion derived by [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] for this case is based on recursively updating the product

τ n,ν,θ (x) = φ n,ν,θ (x)ρ n,ν,θ (x) rather than ρ n,ν,θ (x). The probabilistic interpretation of the new term τ n,ν,θ (x) is E ν,θ [ ( n t=1 s(X t-1 , X t , Y t )) ½{X n = x}| Y 0:n ].
By multiplying (9) by φn+1 (x) and using (8), one obtains the following online update

τn+1 (x) = γ n+1 x ′ ∈X s(x ′ , x, Y n+1 ) φn (x ′ )q θn (x ′ , x)g θn (x ′ , Y n+1 ) x ′ ,x ′′ ∈X 2 φn (x ′ )q θn (x ′ , x ′′ )g θn (x ′′ , Y n+1 ) + (1 -γ n+1 ) x ′ ∈X τn (x ′ ) q θn (x ′ , x)g θn (x ′ , Y n+1 ) x ′ ,x ′′ ∈X 2 φn (x ′ )q θn (x ′ , x ′′ )g θn (x ′′ , Y n+1 ) , (12)
which coincides with Eqs. ( 15)-( 16) of [START_REF] Mongillo | Online learning with hidden Markov models[END_REF], for the particular choice of complete-data sufficient statistics discussed above. Of course, using either ( 8) and ( 9) or ( 12) is practically equivalent. But the results of Section 4 indicate that φn and ρn have very different limiting behaviors: the auxiliary recursion ( 9) should converge to a fixed deterministic limit while the same is not true for the approximate filtering recursion (8). Hence the use of decreasing step-sizes for (9) only is justifiable while its use in ( 12) is less natural (and/or potentially misleading) as τn should not be expected to converge to a deterministic limit.

Regarding the choice of the step-size, [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] consider the cases where, either, the step-size γ n is small but non-decreasing, which may be useful for tracking potential changes but is not sufficient to guarantee the consistency of the approach. The other option mentioned by [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] is to use γ n = n -1 by analogy with the work of [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] and the case of the batch EM algorithm. The range of step-sizes mentioned in Algorithm 1 is chosen in reference to the theory of stochastic approximation and by analogy with [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF]. As will be illustrated below, in the numerical simulations of Section 5, the choice of γ n = n -1 should definitely be avoided for HMMs and we instead recommend using step-sizes of the form γ n = n -α with α in the interval (0.5,0.8), possibly combined with Polyak-Ruppert averaging (see Section 5.3 below).

Implementation and Numerical Complexity

Regarding the numerical complexity of Algorithm 1, observe that in the case considered by [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] where

s(X t-1 , X t , Y t ) = (½{X t-1 = i, X t = j, Y t = k}) (i,j,k)∈X 2 ×Y , s(X t-1 , X t , Y t ) is a vector of dimension |X | 2 × |Y| (where |•|
denotes the cardinal of the set). Thus, the numerical complexity of ( 12) is of order |X | 4 × |Y| per observation. For this case, it is indeed possible to bring down the numerical complexity to the order of |X | 4 + |X | 3 × |Y| operations by updating separately the terms corresponding to the two statistics (½{X t-1 = i, X t = j}) (i,j)∈X 2 and (½{X t = j, Y t = k}) (j,k)∈X ×Y (see the example considered in Section 5 for more details). Interestingly, the numerical complexity of the batch EM algorithm for this model, when implemented using traditional forward-backward smoothing [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF], is of the order of (|X | 2 + |X |×|Y|) per observation and per iteration of the EM algorithm. The comparison is not directly meaningful as the batch EM algorithm does necessitate several iterations to converge (see numerical illustrations in Section 5.3). The scaling of the numerical complexity of the online-EM algorithm with respect to |X | can constitute an hindrance in models with a large number of states. This being said, the complexity of online gradient-based approaches, is equivalent as the main burden comes from the necessity of updating, via a recursion related to (9), one coordinate of the gradient for each of the pairs (x, x ′ ) ∈ X 2 of state values (see, e.g., Le [START_REF] Le Gland | Recursive estimation in HMMs[END_REF]. When the transition matrix is structured -i.e., parametered by a low dimensional parameter rather than by all its individual entries-, the numerical cost of implementing the online EM approach is reduced to an order of the number of parameters times |X | 2 .

Some Results on Convergence

In this section, we provide two important elements regarding Algorithm 1. The first is a generalization of the argument of [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF] that makes it possible to identify a deterministic limit for the EM update as the number of observations tends to infinity. Interestingly, the form of this limit is non trivial and substantially different from the case of independent observations. The second result pertains to the limiting behavior of the auxiliary quantity ρn (x) of ( 10) that is instrumental in Algorithm 1. As discussed in the previous section, a complete analysis of the convergence of Algorithm 1 is still lacking but it is possible to show that when the parameter is frozen (i.e. when the M-step update of (10) is inhibited), ρn converges to a deterministic quantity. This limit does not depend on x (or, in other words, |ρ n (x)ρn (x ′ )| → 0) and is related to the limiting EM mapping. This result, although limited, is very important to understand the nature of the asymptotic attractors of Algorithm 1. In this section, we will work under the assumptions of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF] which guarantee the asymptotic normality of the MLE, adapted to the (simplest) case of a finite state space X .

We start by a closer inspection of the limiting behavior of the normalized score function (gradient of the log-likelihood). Under suitable assumptions (see below), the normalized HMM log-likelihood 1 n log ℓ ν,θ (Y 0 , . . . , Y n ) converges, P θ⋆ almost surely and in L 1 , to the limiting contrast

c θ⋆ (θ) = E θ⋆ [log ℓ θ (Y 0 |Y -∞:-1 )] . ( 13
)
The same is true for the normalized score 1 n ∇ θ log ℓ ν,θ (Y 0 , . . . , Y n ) which converges to ∇ θ c θ⋆ (θ). Such consistency results have been established, under various assumptions, by (among others) [START_REF] Baum | Statistical inference for probabilistic functions of finite state Markov chains[END_REF], [START_REF] Bickel | Asymptotic normality of the maximum likelihood estimator for general hidden Markov models[END_REF], [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF]. Now, thanks to Fisher identity, for all n,

1 n ∇ θ log ℓ ν,θ (Y 0 , . . . , Y n ) = 1 n E ν,θ n t=1 ∇ θ log p θ (X t , Y t |X t-1 ) Y 0:n + 1 n E ν,θ [ ∇ θ log p ν,θ (X 0 , Y 0 )| Y 0:n ] . ( 14
)
As already discussed, the last term on the r.h.s., whose influence is vanishing with increasing values of n, can be ignored. Hence, the consistency result for the score function combined with ( 14) implies that

1 n E ν,θ [ n t=1 ∇ θ log p θ (X t , Y t |X t-1 )| Y 0:n ] converges P θ⋆ almost surely to ∇ θ c θ⋆ (θ)
, the gradient of the limiting contrast. Now using the exponential family representation in (1), the non-vanishing term in the r.h.s. of ( 14) may rewritten as

1 n E ν,θ n t=1 ∇ θ log p θ (X t , Y t |X t-1 ) Y 0:n = ∇ θ ψ(θ) 1 n E ν,θ n t=1 s(X t-1 , X t , Y t ) Y 0:n -∇ θ A(θ) . ( 15
)
The following theorem defines the limiting behavior of the r.h.s. of the above equation, and thus, the limiting behavior of the EM update for HMMs (see Appendix A for the corresponding proofs).

Theorem 1. In addition to Assumption 1, assume that (i) X is a finite set; (ii) the parameter space Θ is compact and θ ⋆ lies in the interior of Θ; (iii) the transition matrix is such that

q θ (x, x ′ ) ≥ ǫ > 0 for all θ ∈ Θ; (iv) sup θ sup y ḡθ (y) < ∞ and E θ⋆ [|log inf θ ḡθ (Y 0 )|] < ∞,
where ḡθ (y) =

x g θ (x, y); and, (v) ψ q , A q , ψ g , A g in (2) are continuously differentiable functions on the interior of Θ. Then the following properties hold.

(i)

1 n E ν,θ n t=1 s(X t-1 , X t , Y t ) Y 0:n -→ E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ]) , P θ⋆ a.s. ( 16
)
(ii) The fixed points of the limiting EM algorithm

θ k+1 = θ {E θ⋆ (E θ k [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ])} (17)
are the stationary points of the limiting likelihood contrast c θ⋆ (θ).

Theorem 1 provides an interpretation of the limiting form of the classical EM algorithm when used on very long sequences of observations. In the case of HMMs, this form is quite complicated and it is interesting to compare it to the case of independent observations investigated by [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF]. If the observations are assumed to be independent, the law of large number implies the convergence of the normalized log-likelihood to c θ⋆ (θ) = E θ⋆ [log ℓ θ (Y 0 )]. In this case, it is obvious that maximizing c θ⋆ (θ) is equivalent to minimizing the Kullback-Leibler divergence D(ℓ θ⋆ |ℓ θ ) = log ℓ θ⋆ (y) ℓ θ (y) ℓ θ⋆ (y)dy. In the HMM case, D(ℓ θ⋆ |ℓ θ ) needs to be replaced by the expression of c θ⋆ (θ) given in ( 13), which is a consequence of the tower property of conditional expectation and of the forgetting property of the filter [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF]. However, in contrast to the case of independent observations, it is no more straightforward to provide an explicit expression for the gradient of c θ⋆ (θ). The key idea here is the use of Fisher identity in ( 14) which yields the limiting term found in the r.h.s. of ( 16). The use of Fisher identity also explains the particular form of conditioning found in (16), which involves both the infinite future and past of the trajectory. This expression also suggests that being able to compute or approximate smoothed functionals of the state recursively is indeed a key requirement for online estimation in HMMs.

The next result shows that under parameter freeze, the online EM update equation (8) converges to a deterministic constant function, equal to the r.h.s. of ( 16).

Corollary 1. Under the assumptions of Theorem 1, Algorithm 1 used without (10) -that is, with θn equal to a fixed value θsatisfies

ρn (x) -→ E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ]) , P θ⋆ a.s., for all x ∈ X .
Corollary 1 shows that under parameter freeze the auxiliary quantity ρn (x) converges to a constant limit that does not depend on x and that is equal to the limit obtained in ( 16). Note that since φn (x) is a probability on X , this implies that θ

( x∈X ρn+1 (x) φn+1 (x)) tends to θ {E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ])}.
This result together with Theorem 1-(ii) suggests that Algorithm 1 can only be stable at the stationary points of the limiting contrast c θ⋆ (θ). The argument is only heuristic at this stage as Corollary 1 is obtained under the artificial assumption that the evolution of the parameter is frozen. Corollary 1 however highlights an essential characteristic of Algorithm 1: in contrast to φn (x) which does approximate the infinite past filter and hence varies with each observation, ρn (x) converges to a deterministic limit that is independent of x. This property, which can be verified in simulations, justifies the form of Algorithm 1 and in particular the use of a stochastic approximation type of update for ρn (x) only.

Application to Gaussian HMMs

HMM with Product Parameterization

For the sake of concreteness, we consider in following the case where the state variables (X t ) take their values in the set {1, . . . , m}. In addition, assume that, as is often the case in practise, the parameter θ may be split into two sub-components that correspond, respectively, to the state transition matrix q θ and to the state-conditional densities {g θ (i, •)} 1≤i≤m . In the fully discrete case considered in [START_REF] Mongillo | Online learning with hidden Markov models[END_REF] for instance, the parameter θ consist of the transition matrices q θ and g θ parametered by their respective entries, with the constraint that each line of a transition matrix must sum to one. In the case of Gaussian HMMs used in speech processing as well as many in other applications, the parameters are the state transition matrix q θ and the mean vector and covariance matrix associated with each of the m state-conditional densities {g θ (i, •)} 1≤i≤m [START_REF] Rabiner | A tutorial on hidden Markov models and selected applications in speech recognition[END_REF].

In such a model, there are two distinct types of EM complete-data sufficient statistics which give rise to two separate forms of the auxiliary function ρ n,ν,θ :

ρ q n,ν,θ (i, j, k; θ) = 1 n E ν,θ n t=1 ½{X t-1 = i, X t = j} Y 0:n , X n = k , (18)
ρ g n,ν,θ (i, k; θ) = 1 n E ν,θ n t=0 ½{X t = i}s(Y t ) Y 0:n , X n = k , (19)
where the form of s itself depend on the nature of the state-conditional distribution g θ (x, •) -see Gaussian example below. There's a slight difference between ( 19) and ( 5), which is that (19) also incorporates the initial (t = 0) conditional likelihood term, i.e., the contribution corresponding to the rightmost term on the r.h.s. of ( 14). As noted earlier, this difference is minor and does not modify the long-term behavior of the algorithm. With these notations, Eq. ( 9) in Algorithm 1 is implemented as

ρq n+1 (i, j, k) = γ n+1 δ(j -k)r n+1 (i|j) + (1 -γ n+1 ) m k ′ =1 ρq n (i, j, k ′ )r n+1 (k ′ |k) , (20)
ρg n+1 (i, k) = γ n+1 δ(i -k)s(Y n+1 ) + (1 -γ n+1 ) m k ′ =1 ρg n (i, k ′ )r n+1 (k ′ |k) , (21)
where δ denotes the Kronecker delta (i.e., δ(i) = 0 iff i = 0) and the notation rn+1 (i|j) refers to the approximate retrospective conditional probability :

rn+1 (i|j) = φn (i)q θn (i, j) m i ′ =1 φn (i ′ )q θn (i ′ , j) . ( 22
)
A complete iteration of the online algorithm involves the approximate filter update (8) and the statistics updates (20) and ( 21) followed by an application of the M-step function θ to Ŝq

n+1 (i, j) = m k=1 ρq n+1 (i, j, k) φn+1 (k) and Ŝg n+1 (i) = m k=1 ρg n+1 (i, k) φn+1 (k).
The form of the M-step depends on the exact nature of q θ and g θ . If the transition matrix q θ is parametered simply by its entries, the update is generic and is given by

q θn (i, j) = θ Ŝq n = Ŝq n (i, j) m j=1 Ŝq n (i, j) . ( 23
)
For the update of the state-dependent parameters, one needs to be more specific and the form of the equations depend on the choice of the state conditional density g θ (x, •). In the multivariate Gaussian case, the function s has to be chosen such that s(y) consists of the three components {1, y, yy t }. The corresponding components of the approximated EM extended statistics are denoted, respectively, by Ŝg n,0 , Ŝg n,1 , Ŝg n,2 . If the state conditional Gaussian densities are parametered by their mean vectors, µ θ (i), and covariances matrices, Σ θ (i), the M-step update is defined as

µ θn (i) = θ Ŝg n,0 , Ŝg n,1 = Ŝg n,1 (i) Ŝg n,0 (i) , (24)
Σ θn (i) = θ Ŝg n,0 , Ŝg n,1 , Ŝg n,2 = Ŝg n,2 (i) Ŝg n,0 (i) -µ θn (i)µ t θn (i) . ( 25
)
The derivation of ( 23) and ( 24)-( 25) is straightforward but some more details are provided in the next section for a particular case of Gaussian HMM.

Markov Chain Observed in Gaussian Noise

In the numerical experiments described below, we consider the simple scalar model

Y t = X t + V t ,
where (V t ) is a scalar additive Gaussian noise of variance υ and (X t ) is a Markov chain with transition matrix q, which takes its values in the set {µ(1), . . . , µ(m)}. Although simple, this model is already statistically challenging and is of some importance in several applications, in particular, as a basic model for ion channels data [START_REF] Chung | Characterization of single channel currrents using digital signal processing techniques based on hidden Markov models[END_REF]) -see also, e.g., [START_REF] Roberts | An EM algorithm for ion-channel current estimation[END_REF] for discussion of the continuous time version of the model as well [START_REF] De Gunst | Statistical analysis of ion channel data using hidden Markov models with correlated state-dependent noise and filtering[END_REF], for an up-to-date account of models for ion channels. The parameter θ is comprised of the transition matrix q, the vector of means µ and the variance υ.

In this case, the intermediate quantity of the batch EM algorithm may be written -with constant terms omitted, as m i=1 m j=1 S q n (i, j) log q(i, j) -

1 2υ m i=1 S g n,2 (i; θ) -2µ(i)S g n,1 (i) + µ 2 (i)S g n,0 (i) , (26)
where

S q n (i, j) = 1 n E ν,θ n t=1 ½{X t-1 = i, X t = j} Y 0:n , S g n,0 (i) = 1 n E ν,θ n t=0 ½{X t = i} Y 0:n , S g n,1 (i) = 1 n E ν,θ n t=0 ½{X t = i}Y t Y 0:n , S g n,2 (i) = 1 n E ν,θ n t=0 ½{X t = i}Y 2 t Y 0:n .
Maximization of (26) with respect to q θ , µ θ and v θ directly yields (23) as well as

µ(i) = θ S g n,0 , S g n,1 = S g n,1 (i) S g n,0 (i) , (27)
υ = θ S g n,0 , S g n,1 , S g n,2 = m i=1 S g n,2 (i) -µ 2 (i)S g n,0 (i) m i=1 S g n,0 (i) . ( 28
)
It is easily checked that, as usual, the M-step equations ( 23) and ( 28) satisfy the constraints that q be a stochastic matrix and υ be non-negative. Note that for this particular model, the use of the statistic S g n,2 could be avoided as it is only needed in the M-step under the form m i=1 S g n,2 (i), which is equal to

1 n n t=0 Y 2 t .
Algorithm 2 below recaps the complete online EM algorithm pertaining to this example.

Algorithm 2 (Online EM algorithm for noisily observed m-state Markov chain).

Initialization Select θ0 and compute, for all 1 ≤ i, j, k, ≤ m and 0

≤ d ≤ 2, φ0 (k) = ν(k)g θ0 (k, Y 0 ) m k ′ =1 g θ0 (k, Y 0 ) , ρq 0 (i, j, k) = 0 , ρg 0,d (i, k) = δ(i -k)Y d 0 .
Recursion For n ≥ 0, and

1 ≤ i, j, k, ≤ m, 0 ≤ d ≤ 2, Approx. Filter Update φn+1 (k) = m k ′ =1 φn (k ′ )q n (k ′ , k)g θn (k, Y n+1 ) m k ′ ,k ′′ =1 φn (k ′ )q n (k ′ , k ′′ )g θn (k ′′ , Y n+1 )
, where g θn (k, y) = exp -(yμn (k)) 2 /2υ n .

Stochastic Approximation E-step

ρq n+1 (i, j, k) = γ n+1 δ(j -k)r n+1 (i|j) + (1 -γ n+1 ) m k ′ =1 ρq n (i, j, k ′ )r n+1 (k ′ |k) , ρg n+1,d (i, k) = γ n+1 δ(i -k)Y d n+1 + (1 -γ n+1 ) m k ′ =1 ρg n,d (i, k ′ )r n+1 (k ′ |k) ,
where rn+1 (i|j) = φn (i)q n (i, j)

m i ′ =1 φn (i ′ )q n (i ′ , j). M-step If n ≥ n min , Ŝq n+1 (i, j) = m k ′ =1 ρq n+1 (i, j, k ′ ) φn+1 (k ′ ) , qn+1 (i, j) = Ŝq n+1 (i, j) m j ′ =1 Ŝq n+1 (i, j ′ ) , Ŝg n+1,d (i) = m k ′ =1 ρg n+1,d (i, k ′ ) φn+1 (k ′ ) , μn+1 (i) = S g n+1,1 (i) S g n+1,0 (i) , υn+1 = m i ′ =1 S g n+1,2 (i ′ ) -μ2 n+1 (i ′ )S g n+1,0 (i ′ ) m i ′ =1 S g n+1,0 (i ′ )
.

Numerical Experiments

Algorithm 2 is considered in the case of a two-state (m = 2) model estimated from trajectories simulated from the model with parameters

q ⋆ (1, 1) = 0.95 , µ ⋆ (1) = 0 , q ⋆ (2, 2) = 0.7 , µ ⋆ (2) = 1 , υ ⋆ = 0.5 .
With these parameters, state identification is a difficult task as the separation of the means corresponding to the two states is only 1.4 times the noise standard deviation. The optimal filter associated with the actual parameter does for instance misclassify the state (using Bayes' rule) in about 10.3% of the cases. As will be seen below, this is reflected in slow convergence of the EM algorithm.

All estimation algorithms are systematically started from the initial values

q ⋆ (1, 1) = 0.7 , µ ⋆ (1) = -0.5 , q ⋆ (2, 2) = 0.5 , µ ⋆ (2) = 0.5 , υ ⋆ = 2 ,
and run on 100 independent trajectories simulated from the model. Figure 1 illustrates the consequences of Theorem 1 by plotting the estimates of the parameters q(1, 1) and µ(1) obtained by the batch EM algorithm, as a function of the number of EM iterations, for two different sample sizes: n = 500 (dotted lines) and n = 8000 iterations. To give an idea of the variability of the estimates, Figure 1 feature the median estimate (bold line) as well as the lower and upper quartiles (lighter curves) for both sample sizes. The first striking observation is the slow convergence of EM in this case, which requires about 50 iterations or so to reach decent estimates of the parameters. When comparing the curves corresponding to the two samples sizes, it is also obvious that while the variability is greatly reduced for n = 8000 compared to n = 500, the median learning curve is very similar in both cases. Furthermore, the plots corresponding to n = 8000 provide a very clear picture of the deterministic limiting EM trajectory, whose existence is guaranteed by Theorem 1. Indeed, the large sample behavior of the batch EM algorithm is rather disappointing as using a fixed number of iteration of EM does involve a computational cost that grows proportionally to n but will converge, as n grows, to a deterministic limit which only depends on the parameter initialization. This is all the more regrettable that from a statistical perspective, it is expected that the true maximum likelihood estimator converges, at rate n -1/2 towards the actual value θ ⋆ of the parameter. This behavior of the batch EM algorithm is illustrated on Figure 2 which displays, from left to right, the estimation results for the parameters associated with the first component (q(1, 1), µ(1)), as in Figure 1, together with the noise variance υ (rightmost display). The results obtained on the 100 independent runs are here summarized as box and whiskers plots. Figure 2, which should be compared with Figure 3 below, shows that when using a potentially large (here, 50) but fixed number of iterations the variability of the batch EM estimates does decrease but the accuracy does not improve as n grows. Clearly, statistical consistency could only be achieved by using more batch EM iterations as the number n of observations grows.

In contrast, Figure 3 which corresponds to the online EM algorithm (Algorithm 2) used with γ n = n -0.6 does suggest that online EM estimation is consistent. For the smallest sample sizes (n = 500 or n = 2000), the estimation variance is still quite large and the online estimates are not as good as those obtained using 50 batch EM iterations. But for sample sizes of n = 8000 and larger, the online EM estimates are preferable despite their somewhat larger variance. In this application, the choice of a slowly decreasing step-size appears to be of utmost importance. In particular, the choice γ n = n -1 , despite its strong analogy with the batch EM case (see Section 2 as well as [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF] provides estimates that are not robust enough with respect to the choice of the initial parameter guess. Simple step-size schemes of the form γ n = n -α perform very poorly in this example when α is set to values higher than 0.8. Figure 4 features a more complex choice of step-sizes which was hand-tuned based on pilot runs. Even with this preliminary phase of step-size tuning -which would hardly be feasible in real-life applications of the method- Figure 4: Estimation results when using the online EM algorithm with γ n = 0.01 for n ≤ n 0 and γ n = 0.5(n -n 0 ) -1 for values of n > n 0 , with n 0 = 10000. From left to right, estimated values of q(1, 1), µ(1) and υ for values of n ranging from 0.5 to 128 thousands of samples. Box and whiskers plot based on 100 independent runs.

it is observed on Figure 4 that the results are not significantly improved compared to the simple choice of γ n = n -0.6 as in Figure 3. The comparatively higher variability observed for the largest sample sizes (32 and 128 thousands) on Figure 4 is caused by a very slow reduction of the bias of the parameter estimates when γ n = n -1 . In contrast, parameter trajectories corresponding to the choice of γ n = n -0.6 look much less smooth but provide more reliable estimates, without requiring model-specific tuning of the step-sizes. This observation -which is also true, but to a lesser extent, in the case of independent observations-is certainly a consequence of the temporal dependence between the observations and, correlatively, of the time needed for the filtering and smoothing relations to forget their initial state. The comparison between Figures 2 and3 is not meant to be fair in terms of computing time, as shown by Table 1 1 . The 50 batch EM iterations used to produce Figure 2 take about 10 to 40, depending on the implementation of batch EM, times longer than for the corresponding online estimates of Figure 3. Being fair in this respect would have mean using just five batch EM iterations which, as can be guessed from Figure 1, is not competitive with online EM, even for the smallest sample sizes. Note that a different option would have been to also consider running the online EM algorithm several times on the same batch data, following [START_REF] Neal | A view of the EM algorithm that justifies incremental, sparse, and other variants[END_REF]. In the case of hidden Markov models however, this way of using the online EM algorithm for fixed-sample maximum likelihood estimation of the parameters appears to be less straightforward than in the case of i.i.d. data and has not been considered.

Method

The two batch EM implementations featured in Table 1 correspond, respectively, to the use of the recursive form of smoothing based on Proposition 1 and to the usual forward-backward form of smoothing. The former implementation is obviously related to the online EM algorithm, which explains that both of them lead to rather similar running times. As discussed in Section 3, due to the fact that the whole m × m transition matrix q is here used as a parameter, the numerical complexity of the online EM algorithm and of the recursive implementation of batch EM scale as m 4 , compared to m 2 only for the batch EM algorithm when implemented with forward-backward smoothing. Hence, it is to be expected that the forward-backward implementation of batch EM would be even more advisable for models with more than m = 2 states. On the other hand, when m is large it is usually not reasonable to parameterize the transition matrix q by its individual entries.

In order to provide a more detailed idea of the asymptotic performance of the algorithm, Figure 5 displays results similar to those of Figure 3 but centered and scaled as follows. Each parameter estimates, say θn is represented as √ n( θn -θ ⋆ ) and further scaled by the asymptotic standard deviation of θ deduced from the inverse of the Fisher information matrix. The Fisher information matrix has been estimated numerically by applying Fisher's identity to (26) so as to obtain

1 n ∇ q(i,j) log ℓ θ (Y 0:n ) = S q n (i, j) q(i, j) - S q n (i, m) q(i, m) (for 1 ≤ j < m) , 1 n ∇ µ(i) log ℓ θ (Y 0:n ) = S g n,1 -µ(i)S g n,0
υ ,

1 n ∇ v log ℓ θ (Y 0:n ) = m i=1
1 2v 2 S g n,2 (i; θ) -2µ(i)S g n,1 (i) + µ 2 (i)S g n,0 (i) .

The information matrix has then been estimated by averaging the gradient computed in θ ⋆ for 100 independent sequences of length one million simulated under θ ⋆ . Additionally, Figure 5 also displays results that have been post-processed using Polyak-Ruppert averaging [START_REF] Polyak | A new method of stochastic approximation type[END_REF][START_REF] Ruppert | Efficient estimation from a slowly convergent Robbins-Monro process[END_REF]. In Figure 5, Polyak-Ruppert averaging is used starting from n avg = 8000. That is, for n > 8000, θn is replaced by 1/(n-8000) n i=8001 θn . For time indices n smaller than 8000, averaging is not performed and the estimates are thus as in Figure 3, except for the centering and the scaling. Under relatively mild assumptions, averaging has been shown to improve the asymptotic rate of convergence of stochastic approximation algorithm making it possible to recover the optimal rate of convergence of n -1/2 . At least for µ(1) and υ, Figure 5 suggests that in this example the proposed algorithm does reach asymptotic efficiency, i.e., becomes asymptotically equivalent to the maximum likelihood estimator. For q(1, 1) the picture is less clear as the recentered and scaled estimates present a negative bias which disappears quite slowly. This effect is however typical of the practical trade-off involved in the choice of the index n avg where averaging is started. To allow for a significant variance reduction, n avg should not be too large. On the other hand, if averaging is started too early, forgetting of the initial guess of the parameters occurs quite slowly. In the present case, the negative bias visible on the left panel of Figure 5 is due to n avg being too small (see corresponding panel in Figure 3). Although, this could be corrected Figure 5: Estimation results when using the online EM algorithm with γ n = n -0.6 with Polyak-Ruppert averaging started after n = 8000. From left to right, estimated values of q(1, 1), µ(1) and υ for values of n ranging from 0.5 to 128 thousands of samples. The estimated values are centered and scaled so as to be comparable with a unitary asymptotic standard deviation. Box and whiskers plot based on 100 independent runs. here by setting n avg to twenty thousands or more, it is important to underline that optimally setting n avg is usually not feasible in practice.

Conclusions

The algorithm proposed in this paper for online estimation of HMM parameters is based on two ideas. The first, which is inspired by [START_REF] Sato | Convergence of on-line EM algorithm[END_REF] and [START_REF] Cappé | On-line expectation-maximization algorithm for latent data models[END_REF], consists in reparameterizing the model in the space of sufficient statistics and approximating the limiting EM recursion by a procedure resembling stochastic approximation. Theorem 1 provides a first argument demonstrating that this idea can also be fruitful in the case of HMMs. The second element is more specific to HMMs and relies on the recursive implementation of smoothing computations for sum functionals of the hidden state which is provided by Proposition 1. As discussed in Section 2, this possibility requires that the auxiliary quantity ρ n,ν,θ defined in (5) be approximated during the course of the algorithm.

Although the performance reported in Section 5 is encouraging, there are several questions raised by this approach. The first is of course the theoretical analysis of the convergence of Algorithm 1, which is still missing. Although originally inspired by stochastic approximation ideas, it seems that Algorithm 1 would be difficult to analyze using currently available stochastic approximation results due to the backward kernel operator rn+1 (x ′ |x) involved in (9). As discussed in Sections 3 and 5.3, the proposed algorithm may become less attractive, from a computational point of view, when used in models with many distinct state values. In such cases and, more generally, in cases where the state-space X of the hidden chain is no longer finite, a promising approach consists in using some form of Monte Carlo computation to approximate (8) and (9). [START_REF] Cappé | Online sequential Monte Carlo EM algorithm[END_REF][START_REF] Del Moral | Forward smoothing using sequential Monte Carlo[END_REF] report encouraging first results in that direction.

A Proofs

Theorem 1 mainly relies on the use of a two-sided forgetting result which is first proved in Corollary 2 below. This result generalizes the one-sided forgetting bounds of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF] and allows conditioning with respect to both past and future observations, which is required for studying the asymptotic behavior of ( 14) and related quantities. The proof of Theorem 1 then mostly relies on the results of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF].

Lemma 1. Given q a transition matrix on the finite set X such that q(x, x ′ ) ≥ ǫ > 0 and α and β probabilities on X , define

J α,q,β (x, x ′ ) = α(x)q(x, x ′ )β(x ′ ) x,x ′ ∈X 2 α(x)q(x, x ′ )β(x ′ ) . Then J α 1 ,q,β 1 -J α 2 ,q,β 2 1 ≤ 1 ǫ ( α 1 -α 2 1 + β 1 -β 2 1 ) , (29)
where µ 1 = x |µ(x)| denotes the L 1 or total variation norm.

Proof. Lemma 1 is obviously related to the application of Bayes' formula. Hence, one may apply Lemma 3.6 of [START_REF] Künsch | State space and hidden Markov models[END_REF]

to obtain J α 1 ,q,β 1 -J α 2 ,q,β 2 1 ≤ 1 ǫ α 1 ⊗ β 1 -α 2 ⊗ β 2 1 . The r.h.s. of (29) is obtained by noting that |α 1 (x)β 1 (x ′ ) -α 2 (x)β 2 (x ′ )| ≤ |α 1 (x) -α 2 (x)|β 1 (x ′ ) + |β 1 (x ′ ) - β 2 (x ′ )|α 2 (x).
Corollary 2. Under the assumptions of Theorem 1, for any function f such that 0 < f < f ∞ and probabilities µ 1 and µ 2 on X 2 , and any index

1 ≤ t ≤ n, x,x ′ ∈X 2 E θ f (X t-1 , X t )| Y 0:n , X 0 = x, X n = x ′ µ 1 (x, x ′ ) - x,x ′ ∈X 2 E θ [ f (X t-1 , X t )| Y 0:n , X 0 = x, X n = x] µ 2 (x, x ′ ) ≤ f ∞ ǫ ρ t-1 + ρ n-t ,
where ρ = (1 -ǫ).

Proof. First apply Lemma 1 to the familiar forward-backward decomposition α i (x) = P θ (X t-1 = x|Y 0:t-1 , X 0 = x 0,i ) ,

β i (x ′ ) ∝ P θ (Y t:n , X n = x n,i |X t = x ′ ) ,
for i = 1, 2 (where the normalization factor in the second equation is determined by the constraint

x∈X β i (x) = 1) to obtain E θ [ f (X t-1 , X t )| Y 0:n , X 0 = x 0,1 , X n = x n,1 ] -E θ [ f (X t-1 , X t )| Y 0:n , X 0 = x 0,2 , X n = x n,2 ] ≤ f ∞ ǫ ( α 1 -α 2 1 + β 1 -β 2 1 ) ,
observing that P θ (X t-1 = x, X t = x ′ |Y 0:t-1 , X 0 = x 0,i , X n = x n,i ) = J α i ,q θ ,β i . Next use, the onesided forgetting bounds of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF] (Corollary 1 and Eq. ( 20)) to obtain α 1 -α 2 1 ≤ ρ t-1 and β 1 -β 2 1 ≤ ρ n-t . The result of Corollary 2 follow by the general inequality |µ 1 (g) -µ 2 (g)| ≤ 1 2 µ 1 -µ 2 1 sup z 1 ,z 2 ∈Z 2 |g(z 1 ) -g(z 2 )|.

Note that the fact that the backward function β i (x ′ ) = P θ (Y t:n , X n = x n,i |X t = x ′ ) may be normalized to a pseudo-probability, has been used in the above proof. This is generally not the case outside of the context where X is finite [START_REF] Cappé | Inference in Hidden Markov Models[END_REF][START_REF] Briers | Smoothing algorithms for state-space models[END_REF] but it is easily checked that Corollary 2 holds in greater generality under the "strong mixing conditions" discussed in Section 4.3 of [START_REF] Cappé | Inference in Hidden Markov Models[END_REF].

Proof of Theorem 1. Corollary 2 implies that

|E θ [ s(X -1 , X 0 , Y 0 )| Y -n:n ] -E θ [ s(X -1 , X 0 , Y 0 )| Y -m:m ]| ≤ 2 ǫ ρ n M (Y 0 ),
for m ≥ n, where M (y) ≥ sup x,x ′ |s(x, x ′ , y)|. M (y) may be chosen as M (y) = x 1 + |s g (x, y)| due to the fact that s q (x, x ′ ) is a vector of indicator functions when X is finite. As, θ ⋆ lies in the interior of Θ, standard results on exponential family imply that M (Y 0 ) has finite first and second order moments under P θ⋆ . Hence, the a.s. limit of E θ [ s(X -1 , X 0 , Y 0 )| Y -m:m ] as m → ∞, which is denoted by E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ], exists and has finite expectation under P θ⋆ . Similarly,

1 n n t=1 (E ν,θ [ s(X t-1 , X t , Y t )| Y 0:n ] -E θ [ s(X t-1 , X t , Y t )| Y -∞:∞ ]) ≤ 1 nǫ n t=1 ρ t-1 + ρ n-t M (Y t ) .
As E θ⋆ [M (Y 0 )2 ] < ∞, standard applications of Markov inequality and Borel-Cantelli Lemma imply that the r.h.s. of the above expression tends P θ⋆ -a.s. to zero. Hence, the quantities For the second statement, one can check that the assumptions of Theorem 1 imply (A1)-(A3) of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF] as well as a form of (A6)-(A8) 2 . Hence, proceeding as in proof of Theorem 3 of [START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF], shows that (14) converge a.s. to the gradient ∇ θ c θ⋆ (θ) of the limiting contrast defined in (13). Eq. ( 15) combined with the previous result then shows that parameter values θ for which ∇ θ c θ⋆ (θ) vanishes are also such that ∇ θ ψ(θ)

{E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ])} -∇ θ A(θ) = 0, that is, θ {E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ])} = θ.
In order to prove Corollary 1, one first needs the following simple lemma. (ii) r n and 1 n n k=1 e k have the same limit.

Lemma 2-(ii) can be proved using Abel transformation noticing that n k=1 ω n k is by definition equal to one and that for the range of step-sizes considered here, (ω n k ) 1≤k≤n is a non-decreasing sequence (it is equal to n -1 when γ k = k -1 and strictly increasing for sequences of step-sizes that decrease more slowly). Lemma 2 shows that we may concentrate on the case where γ k = k -1 to prove Corollary 1. Applying Proposition 1, ρn is then equal under parameter freeze to ρ n,ν,θ defined in (5). In light of Corollary 2, ρ n,ν,θ has the same behavior as the term 1 n n t=1 E ν,θ [ s(X t-1 , X t , Y t )| Y 0:n ] that was analyzed above when proving Theorem 1.

Figure 1 :

 1 Figure1: Estimated values of q(1, 1) (top) and µ(1) (bottom) as a function of the number of batch EM iterations for n = 500 (dotted lines) and n = 8000 (solid lines) observations. The plot is based on 100 independent runs summarized by the median (bold central line) and the upper and lower quartiles (lighter lines).

Figure 2 :

 2 Figure2: Estimation results when using 50 batch EM iterations. From left to right, estimated values of q(1, 1), µ(1) and υ for values of n ranging from 0.5 to 128 thousands of samples. Box and whiskers plot based on 100 independent runs.

Figure 3 :

 3 Figure3: Estimation results when using the online EM algorithm with γ n = n -0.6 . From left to right, estimated values of q(1, 1), µ(1) and υ for values of n ranging from 0.5 to 128 thousands of samples. Box and whiskers plot based on 100 independent runs.

1n

  E ν,θ [ n t=1 s(X t-1 , X t , Y t )| Y 0:n ] and 1 n n t=1 E θ [ n t=1 s(X -1 , X 0 , Y 0 )| Y -∞:∞] have the same limit, where the latter expression converges to E θ⋆ (E θ [ s(X -1 , X 0 , Y 0 )| Y -∞:∞ ]) by the ergodic theorem. This proves the first assertion of Theorem 1.

Lemma 2 .

 2 Let (e k ) k≥0 denote a deterministic sequence and(γ k ) k≥1 a sequence of step-sizes such that γ k ∈ [0, 1], k≥1 γ k = ∞, and, k≥1 γ 2 k < ∞; define r k = (1 -γ k )r k-1 + γ k e k for k ≥ 1. Then (i) r n = n k=1 ω n k e k + ω n 0 r 0 , where ω n k = γ k n j=1(1 -γ j ).

  which does not depend on the newly available observation Y n+1 . The main argument in proving Proposition 1 is to check that

Table 1 :

 1 Computing times in seconds for a record of length n = 10000 observations (Intel Core 2 E6600 2.4 GHz processor).

	MATLAB 7.7 OCTAVE 3.0

The MATLAB/OCTAVE code used for the simulations is very simple and is available as supplementary material. It is mostly vectorized, except for a loop through the observations, and hence it is expected that the differences in running times are indeed representative, despite the use of an interpreted programming language.

Theorem

[START_REF] Douc | Asymptotic properties of the maximum likelihood estimator in autoregressive models with Markov regime[END_REF] deals with the Hessian of the normalized log-likelihood. As we are only concerned with the gradient here, one can drop the second order conditions in (A6)-(A7). Furthermore, as the assumption of Theorem 1 are supposed to hold uniformly on Θ, the set G can be dropped in (A6)-(A7), which provides a law of large number for the score that holds for all values of θ ∈ Θ.


Mohamed Rabie Oularbi

Francois-Xavier Socheleau

Abdeldjalil Aïssa-El-Bey

Sebastien Houcke

Blind Estimation of The Time-Frequency Activity Rate of OFDMA Signals

Keywords: Vertical handover, time-frequency activity rate, OFDMA

The emerging trend to provide users with ubiquitous seamless wireless access leads to the development of multi-mode terminals able to smartly switch between heterogeneous wireless networks. This switching process known as vertical handover requires the terminal to first measure various network metrics. A relevant metric is the allocation rate of network physical channels. In the context of OFDMA and based on the physical signal properties, we propose a method to estimate this rate. Detailed numerical examples demonstrate the efficiency of the proposed estimator. Comparisons with a classical energy detector are also returned.

I. INTRODUCTION

In parallel to the deployment of new mobile services, there is an emerging trend to provide users with ubiquitous wireless access. This can be made possible by taking advantage of the coexistence of complementary heterogeneous networks such as 3G(LTE), WiMAX, Wifi etc. In such environment, a challenge is to develop multi-mode terminals able to smartly switch from one wireless interface to another while maintaining IP or voice connectivity and required QoS. This switching process is known as vertical handover.

Unlike the horizontal handover between two base stations using the same technology, the signal strength metrics are not suited to trigger a vertical handover. In such a case, metrics based on the downlink effective bandwidth or throughput are preferred. As detailed in [START_REF] Dai | Vertical handover criteria and algorithm in IEEE 802.11 and 820.16 hybrid networks[END_REF], this effective bandwidth can be estimated from the data and the packet error rate that can be deduced from the Signal-to-Interference plus Noise Ratio and from the allocation rate of the cell's physical channels. Orthogonal Frequency-Division Multiple Access (OFDMA) is a multi-user version of the popular Orthogonal frequencydivision multiplexing (OFDM) digital modulation scheme. Multiple access is achieved in OFDMA by assigning subsets of subcarriers to individual users in a given time-slot, in view to support differentiated Quality of Service (QoS), i.e. to control the data rate and error probability individually for each user. When the number of active subcarriers is small, the data traffic should also be. Therefore, providing a satisfying downlink signal strength, it is better for a multimode-terminal to connect on such a base station rather than on one where the data traffic is high (high number of active subcarrier).

In this paper, we focus on the passive estimation of the allocation rate of OFDMA physical channels time-frequency slots. The allocation rate is defined as the number of active slots (allocated symbols) divided by the total number of slots per frame. In the context of vertical handover, only the passive estimation is relevant since the terminal seeks to know a priori if a network satisfies its bandwidth needs and not after wasting time and power to get connected to this network.

In some networks such as WiMAX, the physical channels allocation rate is regularly broadcasted by the base station so that it can be known by any terminal. However, this requires a multi-mode terminal that listens to the surrounding networks to intercept every frame preamble. If the multi-mode terminal has to decode every intercepted preamble to get this information, the vertical handover can be a very time and power consuming process.

An alternative approach developed in this contribution is to get the OFDMA physical channels allocation rate by blindly estimating the time-frequency activity rate of OFDMA physical signals. Such approach focuses on the signal properties and therefore does not require any message decoding (assuming this message is made available by the base station which may not be the case in all OFDMA networks). To the best of our knowledge, there is no algorithm published to date that addresses the blind estimation of the time-frequency activity rate of OFDMA signals.

The paper is organized as follows: Section II presents the OFDMA signal model. The proposed estimation algorithm is detailed in section III. The method is based on the estimation of the first and second order moments of the received signal. The performance of this estimator is assessed through simulations in Section IV. Finally, conclusions are presented in Section V.

II. SIGNAL MODEL

Assuming that an OFDMA symbol consists of up to 𝑁 active subcarriers, the discrete-time baseband equivalent transmitted signal is given by

𝑥(𝑚) = √ 𝐸 𝑠 𝑁 ∑ 𝑘∈ℤ 𝑁 -1 ∑ 𝑛=0 𝜀 𝑘,𝑛 𝑎 𝑘,𝑛 𝑒 2𝑖𝜋 𝑛 𝑁 (𝑚-𝐷-𝑘(𝑁 +𝐷)) (1)
.𝑔(𝑚 -𝑘(𝑁 + 𝐷)).

𝐸 𝑠 is the signal power, 𝑎 𝑘,𝑛 are the transmitted data symbols at the 𝑛-th subcarrier of the 𝑘-th ODFM block. These data symbols are assumed to be independent identically distributed (i.i.d), 𝐷 is the cyclic prefix length; 𝑚 → 𝑔(𝑚) is the pulse shaping filter. 𝜀 𝑘,𝑛 is a set of i.i.d random variable valued in {0, 1}, expressing the absence or presence of signal activity in the (𝑘, 𝑛) time frequency slot. Let ℎ(𝑙) 𝑙=0,...,𝐿-1 be the baseband equivalent discrete-time channel impulse response of length 𝐿, such that 𝐿 + 1 < 𝐷. The received samples of the OFDMA signals are then expressed as

𝑦(𝑚) = 𝐿 ∑ 𝑙=0 ℎ(𝑙)𝑥(𝑚 -𝑙) + 𝑤(𝑚) (2)
where 𝑤(𝑚) is a zero mean circularly-symmetric complexvalued white Gaussian noise of variance 𝜎 2 𝑤 . We assume that the receiver is synchronized with the transmitter in time and in frequency. This synchronization can be realized thanks to the frame preamble or thanks to blind techniques presented in [START_REF] Van De Beek | ML estimation of time and frequency offset in OFDM systems[END_REF] and [START_REF] Park | A blind OFDM synchronization algorithm based on cyclic correlation[END_REF]. We also assume that the noise variance is known or at least estimated thanks to blind method such as those detailed in [START_REF] Socheleau | Blind noise variance estimation for OFDMA signals[END_REF], [START_REF] Socheleau | Non Data-Aided SNR Estimation of OFDM Signals[END_REF].

III. ESTIMATION ALGORITHM

The estimation of the time-frequency activity rate 𝛼 is equivalent to detect the active slots from the non active ones.

𝛼 = ∑ 𝑘,𝑛 𝐼(𝜀 𝑘,𝑛 = 1) 𝑀𝑁 ( 3
)
where 𝐼(𝐴) is the indicator function of any event 𝐴 and 𝑀 is the number of observed OFDM symbols. Intuitively, considering that 𝜎 2 𝑤 is known, a classic detector structure could be used so that

α = ∑ 𝑘,𝑛 𝐼(|𝑌 𝑘,𝑛 | > 𝑠(𝜎 𝑤 )) 𝑀𝑁 ( 4
)
where 𝑠(𝜎 𝑤 ) is a threshold function and 𝑌 𝑘,𝑛 is the signal observation on the slot (𝑘, 𝑛).

𝑌 𝑘,𝑛 Δ = 1 √ 𝑁 𝑁 -1 ∑ 𝑚=0 𝑦[𝑘(𝑁 + 𝐷) + 𝐷 + 𝑚]𝑒 -2𝑖𝜋 𝑛𝑚 𝑁 (5) = 𝜀 𝑘,𝑛 𝑎 𝑘,𝑛 𝐻 𝑘,𝑛 √ 𝐸 𝑠 + 𝑊 𝑘,𝑛 (6)
where 𝐻 𝑘,𝑛 and 𝑊 𝑘,𝑛 are respectively the channel frequency response at subcarrier 𝑛 and the noise at subcarrier 𝑛 of the 𝑘-th received symbol.

The limitation of such approach is that the performance are strongly impacted by the choice of a threshold. In order to avoid this constraint, we hereafter propose an alternative method.

The proposed technique relies on the first and second order moments estimate of the observed samples absolute value. This moments are indeed dependent of the activity rate 𝛼.

For all (𝑘, 𝑛) such that 𝜀 𝑘,𝑛 = 0, the observations are made of noise-only slots such that they satisfy 𝑌 𝑘,𝑛 ∼ 𝒞𝒩 (0, 𝜎 2 𝑤 ).

Therefore, in this case the absolute value |𝑌 𝑘,𝑛 | has a Rayleigh distribution of expectation

𝔼[|𝑌 𝑘,𝑛 |/𝜀 𝑘,𝑛 = 0] = √ 𝜋 2 𝜎 𝑤 . ( 7
)
where 𝔼[./.] defines the conditional expectation.

When the observations are made of signal plus noise samples (i.e 𝜀 𝑘,𝑛 = 1), their distribution is quite difficult to know in a context of passive estimation. The various 𝑎 𝑘,𝑛 can indeed have different distributions that depend on the symbols constellation (usually BPSK, QPSK, 16-QAM or 64-QAM [START_REF]Part 16: air interface for broadband wireless access systems, Amendment 2 : Physical and Medium Access Control layers for Combined Fixed and Mobile Operation in License Bands and Corrigendum 1[END_REF]). Actual systems employing adaptive modulation, this constellation can be different from a slot to another. Therefore, according to the principle of maximum entropy [START_REF] Jaynes | Probability Theory: The Logic of Science[END_REF], the state of ignorance on the constellation distribution is here modeled by a uniform law. Hence, without prior information, we assume that the probability to get each constellation equals 1/4 (note that the impact of this assumption is discussed in section 4). Consequently the expectation of |𝑌 𝑘,𝑛 | when 𝜀 𝑘,𝑛 = 1 can be written as

𝔼[|𝑌 𝑘,𝑛 | / 𝜀 𝑘,𝑛 = 1] = 𝔼[|𝑎 𝑘,
∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 ] = 𝜑( ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 ) ( 9
)
where 𝜑 is a function that associate to each 𝐸 𝑠 the expectation

𝔼[|𝑌 𝑘,𝑛 |/𝜀 𝑘,𝑛 = 1], when 𝜎 2
𝑤 is assumed to be known. Since 𝛼% of the slots are active and (1 -𝛼)% are not, the expectation of the overall observed samples module is expressed as

𝔼[|𝑌 𝑘,𝑛 |] =𝛼𝜑( ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 ) + (1 -𝛼) √ 𝜋 2 𝜎 𝑤 . ( 10
)
Moreover, the second order moment

𝔼[|𝑌 𝑘,𝑛 | 2 ] satisfies 𝔼[|𝑌 𝑘,𝑛 | 2 ] = 𝜎 2 𝑤 + 𝛼 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 , ∀𝜀 𝑘,𝑛 . ( 11
)
It follows that

∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 is deduced as ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 = 𝔼[|𝑌 𝑘,𝑛 | 2 ] -𝜎 2 𝑤 𝛼 , ( 12
)
If we denote by

𝑚 1 = 𝔼[|𝑌 𝑘,𝑛 |] and 𝑚 2 = 𝔼[|𝑌 𝑘,𝑛 | 2 ] then m1 = 1 𝑀𝑁 𝑀 -1 ∑ 𝑘=0 𝑁 -1 ∑ 𝑛=0 |𝑌 𝑘,𝑛 | (13) m2 = 1 𝑀𝑁 -1 𝑀 -1 ∑ 𝑘=0 𝑁 -1 ∑ 𝑛=0 |𝑌 𝑘,𝑛 | 2 (14)
substituting this value in equation ( 10), an estimate of the channel occupancy rate α can be obtained solving the following equation

α𝜑( m2 -𝜎 2 𝑤 α ) + (1 -α) √ 𝜋 2 𝜎 𝑤 -m1 = 0 (15)
This equation has no analytical solution, we propose to solve it by a binary search algorithm. The whole corresponding technique is presented in algorithm 1. Notice that no channel estimation is required.

Algorithm 1 Moments method

Observe 𝑀 OFDM symbols Estimate 𝜎 2

𝑤

Compute 𝑌 𝑘,𝑛 Compute m1 and m2 thanks to (13) and ( 14) Deduce α solving (15) thanks to the binary search algorithm.

IV. SIMULATIONS

In this section, OFDMA signals with 512 subcarriers are considered. 𝐷 is set to 128, 𝑀 = 24 and the time-frequency slots allocation is supposed i.i.d. Each 𝑎 𝑘,𝑛 is randomly chosen within BPSK, QPSK, 16-QAM and 64-QAM constellations according to a uniform law. These constellations are the main constellations used by the WiMAX Adaptive Modulation and Coding (AMC) [START_REF]Part 16: air interface for broadband wireless access systems, Amendment 2 : Physical and Medium Access Control layers for Combined Fixed and Mobile Operation in License Bands and Corrigendum 1[END_REF].

The estimator performance is assessed in WMAN environment where the channel is highly selective [START_REF] Erceg | Channel Models for Fixed Wireless Applications[END_REF]. Figure 1

∑ 𝐿 𝑙=0 𝔼[|ℎ 𝑘 (𝑙)| 2 ] = 1.
The channel is assumed to be time variant with a Doppler frequency equal to 100 Hz. Figure 1 compares the performance of the estimator in the two cases where 𝜎 𝑤 is first assumed to be perfectly known and when it is estimated thanks to the method presented in [START_REF] Socheleau | Non Data-Aided SNR Estimation of OFDM Signals[END_REF] and [START_REF] Socheleau | Blind noise variance estimation for OFDMA signals[END_REF]. We can observe, that the estimator's performance deteriorates when 𝜎 𝑤 is estimated, but still offers satisfying performance for the targeted application. A NMSE of -15 dB can indeed be considered as sufficiently accurate to decide whether to trigger a handover or not.

In figure 2, the performance of the proposed estimator is compared to the Constant false alarm rate (CFAR) [START_REF] Scharf | Statistical Signal Processing: Detection, Estimation, and Time Series Analysis[END_REF] energy detector per slot. This classical approach, relies on a threshold associated to a false alarm 𝑃 𝑓𝑎 . It assumes that 𝜀

𝑘,𝑛 = 1 if |𝑌 𝑘,𝑛 | 2 > -𝜎 2 𝑤 ln(𝑃 𝑓𝑎 )
The results show the problems induced by the choice of the threshold. We can clearly observe that for a given SNR = 10dB, the choice of the threshold greatly impacts the performance of the CFAR method. The proposed approach offers better and more stable results even when 𝜎 𝑤 is estimated.

Figure 3, compares the performances of the proposed algorithm in the case where the 𝑎 𝑘,𝑛 are not uniformly chosen from one of the possible constellation present in WiMAX. Indeed, we assess the performance when the probability that a chosen symbol belongs to any constellation is P(𝑎 𝑘,𝑛 ∈ BPSK) = 1 10 , P(𝑎 𝑘,𝑛 ∈ QPSK) = 1 10 , P(𝑎 𝑘,𝑛 ∈ 16-QAM) = 4 10 , P(𝑎 𝑘,𝑛 ∈ 64-QAM) = 4 10 whereas the receiver expect a uniform distribution. We can clearly note that the proposed algorithm is robust to a non-equirepartition of the constellation due to the AMC scheme.

V. CONCLUSION

In this paper, a new blind estimation method of OFDMA time-frequency activity rate has been presented. The method is computationally cheap and exhibits good performance for vertical handover. This approach does not rely on a choice of a threshold and show good results compared to the classical CFAR energy detector even when 𝜎 𝑤 is estimated.

Fig. 1 .Fig. 2 .

 12 Fig. 1. NMSE of the proposed method versus the activity rate, 𝜎𝑤 known and estimated.

Fig. 3 .

 3 Fig.3. NMSE of the proposed method in the case of uniform and nonuniform repartition of the constellations.

  Rayleigh fading channel and a known 𝑎 𝑘,𝑛 , the distribution of the observed slots is Gaussian:𝑌 𝑘,𝑛 /𝑎 𝑘,𝑛 , 𝜀 𝑘,𝑛 = 1 ∼ 𝒞𝒩 (0, 𝐸 𝑠 |𝑎 𝑘,𝑛 | 2 + 𝜎 2 𝑤 ). It then follows that the absolute value |𝑌 𝑘,𝑛 /𝑎 𝑘,𝑛 , 𝜀 𝑘,𝑛 = 1| has a Rayleigh distribution. After performing integration over all the possible values of 𝑎 𝑘,𝑛 in each 𝐶 𝑀𝑗 constellation noted 𝑐 𝑝 we find that

	and consequently,
	𝔼[|𝑌 𝑘,𝑛 |/𝜀 𝑘,𝑛 = 1] = √ 𝜋 8 ⎡ ⎣ 5 2 √ ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 + + 1 4 √ 9 5 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 + 1 4 √ ∑ 𝑙 1 16 √ ∑ 𝜎 2 ℎ(𝑙) 𝜎 2 𝐸 𝑠 5 ℎ(𝑙) + 𝜎 2 𝑤 𝐸 𝑠 + 𝜎 2 𝑤 21 𝑙 + 1 8 √ 5 21 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 + 1 16 √ 3 ∑ 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 7 𝑙 + 1 8 √ 13 21 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 + 1 8 √ 17 ∑ 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 21 𝑙 + 3 16 √ 25 21 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 + 1 8 √ 29 ∑ 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 𝑤 21 𝑙 + 13 8 √ 37 21 ∑ 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 + 𝜎 2 13 √ 7 𝑤 + 16 3
	=	1 4
	∑ 𝐿-1 0 ℎ(𝑙) 𝔼[|𝑌 𝑘,𝑛 |/𝜀 𝑘,𝑛 = 1] = 𝜎 2
	√ 𝜋	1	4 ∑	1	𝑀𝑗 ∑	√ ∑	(8)
	2	4	𝑗=1	𝑀 𝑗	𝑝=1	𝑤

𝑛 𝐻 𝑘,𝑛 √ 𝐸 𝑠 + 𝑊 𝑘,𝑛 |] 4 ∑ 𝑗=1 𝔼 [ |𝑎 𝑘,𝑛 𝐻 𝑘,𝑛 √ 𝐸 𝑠 + 𝑊 𝑘,𝑛 |/𝑎 𝑘,𝑛 ∈ 𝒞 𝑀𝑗 ] where the 𝐶 𝑀𝑗 constellations are 𝑀 𝑗 -QAM such that for 𝑗 = 1, . . . , 4 𝑀 𝑗 is respectively equal to 2,4,16,64. Assuming a Gaussian noise, a 𝑙 𝜎 2 ℎ(𝑙) 𝐸 𝑠 |𝑐 𝑝 | 2 + 𝜎 2

  shows the NMSE (Normalized Mean Square Error) of the proposed estimators for different SNR versus the activity rate.

						The SNR
	is defined as	SNR =	𝑆 𝑤 𝜎 2			(16)
	with	𝑆 = 𝐸 𝑠 𝔼[|𝜀 𝑘,𝑛 𝑎 𝑘,𝑛 | 2 ]	𝐿-1 ∑	𝜎 2 ℎ(𝑙)	(17)
				𝑙=0
	The propagation channel {ℎ 𝑘 (𝑙)} 𝑙=0,⋅⋅⋅ ,𝐿 has an exponential decay profile for its non-null component (i.e., 𝔼[|ℎ 𝑘 (𝑙)| 2 ] = 𝐺𝑒 -𝑙/𝜇 for 𝑙 = 0, ⋅ ⋅ ⋅ , 𝐿) with 𝐿 = 𝐷, 𝜇 = 32 and 𝐺
	is chosen such that


Bogdan Ušćumlić

Annie Gravey

Philippe Gravey

Michel Morvan

Scheduling in Optical Packet Rings

Keywords: optical packet ring, scheduling, stability, fairness, MaxWeight

In the metropolitan area, traditional SONET/SDH circuit switched rings are likely to be replaced with optical packet/burst switching technologies. In this paper we consider a slotted WDM optical packet ring operating without resource reservation mechanisms. In such rings, optical packets in transit have priority over traffic to be inserted by the node. Packets to be inserted are thus queued according to their destination, in order to avoid head-of-line blocking. We focus on scheduling policies and compare several MaxWeight scheduling policies, including Oldest Packet First (OPF) which emulates FIFO queueing while avoiding head-of-line blocking. We show that there is a trade-off between implementation complexity and fairness, and identify the Largest Virtual Waiting Time First (LVWTF) scheduling policy as presenting both a low complexity and a good fairness performance.

I. INTRODUCTION

In this paper we consider an optical packet switching ring, designed for metro area, based on WDM technology and using both tunable lasers and Optical Add/Drop Multiplexers (OADM). This network is developed within the ECOFRAME project. 1 .

In the past years, several projects have studied optical packet rings ( [START_REF] White | A summary of the HORNET project: a next-generation metropolitan area network[END_REF], [START_REF] Carena | RingO: an experimental WDM optical packet network for metro applications[END_REF], [START_REF] Dey | FLAMINGO: A Packet-switched IP-over-WDM Alloptical MAN[END_REF], [START_REF] Stavdas | IST-DAVID: concept presentation and physical layer modeling of the metropolitan area network[END_REF], [START_REF] Cai | The multitoken interarrival time (MTIT) access protocol for supporting variable size packets over WDM ring network[END_REF], [START_REF] Hu | Performance of MAC layer and fairness protocol for the dual bus optical ring network (DBORN)[END_REF], [START_REF] Herzog | RINGOSTAR: an evolutionary performanceenhancing WDM upgrade of IEEE 802.17 resilient packet ring[END_REF], [START_REF] Lazaro | Remotely amplified SARDANA: Single-fibre-tree advanced ring-based dense access network architecture[END_REF]). They differ in terms of insertion and extraction methods, in terms of framing issues or in terms of architecture scenarios. In [START_REF] Uscumlic | Impact of Peerto-Peer Traffic on the Efficiency of Optical Packet Rings[END_REF] it is shown that allowing any-to-any traffic in a single-wavelength, unidirectional optical packet ring can improve its capacity compared to a classical concentration/distribution scenario. The positive impact of WDM dimension on performance of optical packet switched rings has also been studied in [START_REF] Uscumlic | Traffic Grooming in WDM Optical Packet Rings[END_REF]. These studies, as the present one, assume that packets are inserted in an opportunistic manner, i.e. as soon as an appropriate slot is identified that can carry the packet to the destination. In other words, we assume that the network has been correctly dimensioned for the offered traffic and that no reservation mechanism is used.

In the present paper, we analyze the performance offered by several scheduling policies, in terms of stability, of implementation complexity and of fairness.

The reminder of the paper is organized as follows. In Section II we give the basic characteristics of the network in study. In Section III we address the choice of suitable Dr Bogdan Ušćumlić is with the Institute Mihailo Pupin, University of Belgrade, Serbia (e-mail: bogdan.uscumlic@institutepupin.com). Dr Ušćumlić is a former PhD student of the Institut TELECOM -TELECOM Bretagne, Brest, France.

Prof Annie Gravey, Prof Philippe Gravey and Michel Morvan are with the Institut TELECOM -TELECOM Bretagne, Brest, France (e-mails: {annie.gravey, philippe.gravey, michel.morvan}@telecom-bretagne.eu) 1 This collaborative project project has been funded by the French National Research Agency between 2007 and 2009.

scheduling policy for optical packet ring. The numerical results comparing different scheduling policies are given in Section IV. Finally, last section concludes the document.

II. THE ECOFRAME OPTICAL PACKET RING

The ECOFRAME Optical Packet Ring is a unidirectional network, where nodes are connected to the ring via packet based OADMs, and to the client layers via an adaptation interface, responsible for multiplexing and demultiplexing client frames into optical packets.

All the operations in the network are synchronized and occur at discrete time periods. Data is transported by using optical containers, so-called DATA packets, while control information is transported on a separate channel in control packets. Control packets report the occupancy status of the current time slot. A time slot on a wavelength can be either busy or free, depending on whether it carries a DATA packet or not.

Each node is equipped with a fully-tunable transmitter, able to dynamically change the transmission wavelength, from slot to slot. An important feature of this architecture is that, although multiple wavelengths are used to carry data packets, a node is only allowed to insert at most one packet per time slot. In the general case, a node can receive packets on a fixed set of wavelengths. Received packets are queued in reception node before being delivered to the client layers.

The ECOFRAME network is transparent to transit traffic. Packets to be inserted at a given node are queued till a free time slot is found on an appropriate wavelength.

In the present paper, we assume that for each pair of nodes, a single wavelength is used to carry all the traffic between the source and the destination nodes. However, the system is still WDM since many wavelengths can be used on the ring. It has been shown in [START_REF] Uscumlic | The Impact of Network Design on Packet Scheduling in Slotted WDM Packet Rings[END_REF] that this policy, although suboptimal in terms of network dimensioning, allows to consider simple per-destination FIFO queues accessing each a single wavelength instead of a complex set of queues accessing several wavelengths.

III. SCHEDULING IN OPTICAL PACKET RINGS

Scheduling can be defined as the process of selecting a packet and a wavelength for insertion, by a ring node, per time slot basis. In the present case, scheduling consists in selecting, per time slot, a non-empty queue and a wavelength on which this queue should send its packets.

Contrarily to opaque systems where transit traffic is electronically handled in each node, and scheduling involves both new traffic and transit traffic, scheduling in ECOFRAME rings only affects new traffic since transit traffic is served transparently. A drawback of transparency is that in the general case, schedulers are "non-work-conserving": if in a given time slot, queued packets cannot be sent on free wavelengths, no new traffic is served.

Transparency offered to transit traffic also impacts the stability of the ring. Indeed, only one packet can be served per time slot, and this packet can only be transmitted on a single wavelength. This implies that the scheduling policy has a significant impact on stability. This is a first point to be addressed.

A scheduling policy is deemed fair if all insertion FIFO queues within a given node present the same delay performance. This does not mean that end-to-end delays are identical since in particular transit delays depend on the respective positions of the nodes on the ring. Fairness is the second point to be addressed when selecting a scheduling policy.

Lastly, it is also important to assess the complexity associated with a given scheduling policy. Complexity is expressed here as the number of states and counters to maintain in each node. Stateless scheduling can be implemented when there is a single queue per node. Another extreme is observed when several counters per packet to be inserted have to be maintained.

A. Impact of scheduling on stability

We have chosen to mostly focus on MaxWeight scheduling disciplines [START_REF] Stolyar | MaxWeight Scheduling in a Generalized Switch: State Space Collapse and Workload Minimization in Heavy Traffic[END_REF]. MaxWeight scheduling guarantees the stability of different instances of generalized switch model, as first shown by Tassiulas & Ephremides in [START_REF] Tassiulas | Dynamic server allocation to parallel queues with randomly varying connectivity[END_REF] and generalized by Stolyar in [START_REF] Stolyar | MaxWeight Scheduling in a Generalized Switch: State Space Collapse and Workload Minimization in Heavy Traffic[END_REF]. Dimensioning methods for ECOFRAME rings relying on MaxWeight scheduling policies are further studied in [START_REF] Uscumlic | Stable dimensioning of optical packet rings[END_REF].

Simple and popular policies may not be MaxWeight, which leads to sub-optimal dimensioning or to instability. We adapt here an example from [START_REF] Tassiulas | Dynamic server allocation to parallel queues with randomly varying connectivity[END_REF] to show that Priority Queueing is non MaxWeight, and presents a smaller stable domain than MaxWeight policies.

Consider a ring with 2 wavelengths. Assume that a given node sends traffic from Queue1 on wavelength Λ 1 and from Queue2 on wavelength Λ 2 ; let λ 1 =0.5 and λ 2 be the arrival rates to Queue1 and Queue2. We further assume that wavelength Λ 1 is always available, which we note with p(Λ 1 ) = 1, whereas a slot is available on wavelength Λ 2 with probability p(Λ 2 ) = 0.5.

According to [START_REF] Tassiulas | Dynamic server allocation to parallel queues with randomly varying connectivity[END_REF], a necessary and sufficient stability condition for this system with two wavelengths, employing a "Longuest Queue First" (LQF) scheduling policy, is given by the following set of inequalities: λ 1 < p(Λ 1 ),λ 2 < p(Λ 2 ), and

λ 1 + λ 2 < p(Λ 1 ) + p(Λ 2 ) -p(Λ 1 )p(Λ 2 ).
With the LQF scheduling policy, the above system is stable as long as λ 2 < 0.5.

Consider now a policy that implements priority to Queue1 (PQ). At each time slot, the node first attempts to transmit a packet from Queue1, and if there are no such packets, the node attempts to transmit a packet from Queue2. Obviously, under the PQ policy, the probability that the server attempts to serve Queue2 is only (1 -λ 1 )p(Λ 2 ) = 0.25, which implies that the system is only stable if λ 2 is smaller than 0.25, which is significantly smaller that the rate accepted by LQF.

This justifies our present focus on MaxWeight policies, for which stability conditions can straightforwardly be specified. Due to lack of space, we do not address here another popular scheduling policy, Round Robin, which is not either MaxWeight, although it is frequently considered.

B. Considered Scheduling Policies

Each packet is characterized by its destination, and by the wavelength it should be transmitted to reach this destination. Therefore, we can classify packets according to:

1) either the destination (destination address queue), 2) or the destination wavelength (wavelength queue). In the present case, we consider queues per destination. Each such queue is FIFO, i.e. we do not consider that insertion is class based within a queue. Note that such approach differs from the one taken in [START_REF] Eido | Packet Filling Optimization in Multiservice Slotted Optical Packet Switching MAN Networks[END_REF], for instance, where the authors consider separate queues for each (class of service, destination) pair. In our study, we try to reduce the number of queues, and thus to simplify the scheduling.

All scheduling mechanisms have to start by identifying, for each time slot, a set of eligible pairs (queue, wavelength). In a given slot, a pair is eligible if the wavelength is free, the queue is not empty and has to send packets on this particular wavelength. If the pair set is not empty, one eligible packet will be selected according to the rules that characterize the scheduling policy.

Four following scheduling policies are now considered: 1) Random (RAND): a queue to be served is determined by uniformly selecting one pair within the set of eligible pairs; 2) Oldest First Packet (OPF): the queue to be served is the one containing the packet that has experienced the maximum waiting time. This particular policy is FIFO within eligible packets and generally differs from a global FIFO which may suffer from head-of-the-line blocking. 3) Longest Queue First (LQF): the queue to be served is the longest queue in the set of eligible pairs; 4) Longest Virtual Waiting Time First (LVWTF): A modified LQF who weights the length of the queue by a factor inversely proportional to the flow rate used for dimensioning the system. LVWTF policy, for each node i, considers the variables size j /a ij for eligible queues, where a ij is the flow rate between nodes i and j (normalized to the wavelength capacity) and serves the queue for which this value is maximum. The above scheduling policies, except RAND, are MaxWeight scheduling disciplines [START_REF] Stolyar | MaxWeight Scheduling in a Generalized Switch: State Space Collapse and Workload Minimization in Heavy Traffic[END_REF].

C. Scheduling Policy performance Metrics

Scheduling policies are compared in terms of end-do-end queueing performance. End-to-end queueing latency is defined as the sum of insertion and extraction times. The insertion process leads to delays, because of the queueing process of data packets prior to insertion to the ring. The extraction time is due to the fact that a given node may receive packets on several wavelengths, whereas it delivers them to the client layers at a rate equal to a wavelength rate.

Fairness of a scheduling policy is defined by comparing the delay performance delivered by a given policy to the one delivered by the OPF policy. Indeed, OPF is the policy closest to a global FIFO, but not suffering from head-of-theline blocking. The fairness metric considered below, although rather primitive, is both well suited to a comparison carried out by simulation, and nevertheless very significant.

Assume a set of flows I 1 , I 2 , ..., I n . Let us note D I1 , D I2 , ..., D In and L I1 , L I2 , ..., L In , the mean delays (in number of time slots) of packets belonging to flows I 1 , I 2 , ..., I n , for OPF and some scheduling policy S, respectively.

The measure of the efficiency for the set of flows I 1 , I 2 , ..., I n , when using the scheduling rule S is defined by:

E S = n i=1 |D Ii -L Ii |.
(

) 1
The above metric is used in the next Section in order to compare the scheduling policies proposed here with OPF.

D. Scheduling Policy Complexity

As stated before, complexity is evaluated here by identifying the number of counters and states to maintain within a node in order to implement a particular scheduling policy. Hardware design is obviously much simpler when the number of counters to implement is fixed and if possible small.

Let us consider the complexity of the above scheduling disciplines:

1) RAND relies on a single occupancy indicator per queue.

2) In addition to the occupancy indicator, LQF and LVWTF necessitate maintaining a length indicator per queue. 3) In addition to the occupancy indicator, OPF necessitates storing arrival time for each packet. Obviously, RAND has a very low complexity. LQF and LVWTF, although slightly more complex, necessitate a fixed number of counters (one per possible destination). On the other hand, the complexity for OPF is high since it necessitates a (potentially unknown) number of counters, each carrying the arrival time of a given packet. Although OPF would be a natural choice in terms of performance, since it closely emulates a global FIFO queue, we wish to assess whether a less complex scheduling policy would present a good fairness performance (i.e. is "close" to OPF in terms of delay performance).

IV. ASSESSMENT OF DELAY PERFORMANCE

Using a ns-2 simulator developed for studying the ECOFRAME network, we have assessed the delay performance delivered by the four scheduling policies listed above.

All the simulation results in this work, are given with a confidence interval of 10% (at confidence level of 95%). In order to assess fairness as defined in (1), the zero measure of efficiency is determined by

0 i = 0.1 • D Ii .In other words, if |D Ii -L Ii | ≤ 0 i , then |D Ii -L Ii |
should be taken to be equal to 0, when calculating E S . The value of E S is then used to assess how close a given scheduling policy is to OPF.

A. Performance for a single-wavelength network

We first consider the simple case where a single wavelength is used by all traffic flows. We consider the network in Fig. 1 and we assess fairness in node A. Let us assume that nodes emit traffic flows as defined in scenario in Fig. 1, for k = 1.

Fig. 2 reports, for different scheduling policies, the average insertion times for flows AB and AC in two cases: when their rate is A → B = A → C = 0.2 and when their rate is A → B = 0.3, A → C = 0.1. Fig. 2(a) shows that when the rates for flows AB and AC are identical, the 4 mechanisms present identical behaviours. On the other hand, when the rates for flows AB and AC differ, the 4 mechanisms behave differently: some favor one flow above the other. Fig. 2(b) shows that OPF delivers the same insertion delay performance to the two flows, whereas LQF favors the larger flow (a result also observed in [START_REF] Mckeown | Achieving 100% Throughput in an Input-Queued Switch[END_REF]). On the other hand, RAND favors the smaller flow, while LVWTF only slightly disadvantages the larger flow.

Tab. I presents the efficiency scores for different scheduling rules and different values of k. Tab. I shows that the deviation of the optimal performance, increases with the value of k, for both LQF and RAND scheduling. This is not the case for LVWTF, which presents only a small deviation from OPF, for all values of k.

As specified above, we have assumed here that each node sends traffic as specified in the dimensioning process. This is not necessarily a plausible situation. In general, the dimensioning process considers "busy hour" estimates to dimension the network, and actual flow rates are lower than these estimates.

We have shows that LVWTF is still very close to OPF in this case; this is not shown here due to lack of space.

It may also happen that the dimensioning process has underestimated some actual flow rates. We argue here that this is not an issue that can be suitably handled by the scheduling policy because the dimensioning itself may not be stable anymore ! Supplementary policies, such as e.g. policing or conformance control mechanisms, are requested to deal with these cases, and to limit the actual flow rates to their agreed limits.

B. Performance for a multi-wavelength network

In order to illustrate this more general case, we consider a 6-node 2-wavelength ring where all stations can listen on both wavelengths. We label the ring nodes with A,B,C,D,E and F, in traffic direction. The traffic matrix is complete, symmetric and uniform of rate x with one exception: node A is the only node sending non-symmetric traffic. Node A sends traffic of load x to stations C,D,E and F, and traffic of load 2x to station B.

In Fig. 3, the expected waiting time for flows AB and AC are compared by simulation in cases of LVWTF and OPF policy, in function of traffic rate A → C. The differences in the latency between OPF and LVWTF remain limited in this scenario, which confirms that LVWTF is a good candidate for approximating OPF.

V. CONCLUSION

In this paper we have addressed the issue of selecting scheduling policies for a slotted WDM optical packet ring. We have first stated that, although desirable in terms of delay performance, FIFO scheduling is not appropriate as it presents head-of-the-line blocking. After pointing out the necessity of selecting a MaxWeight scheduling policy, we have compared the performance delivered by LQF and LVWTF to the one delivered by OPF which closely emulates FIFO.

We have shown that OPF is far more complex to implement than LQF and LVWTF, and that LVWTF, which is as simple as LQF, provides a queueing performance close to the one delivered by OPF. Actually, the differences between OPF, LQF and LVWTF can be shown to decrease when the number of wavelengths in the ring increases. This is due to the multiplexing gain brought by WDM as shown in a previous paper [START_REF] Uscumlic | Traffic Grooming in WDM Optical Packet Rings[END_REF].

Fig. 1 .

 1 Fig. 1. Network used to assess fairness in node A.

1 Fig. 2 .

 12 Fig. 2. Comparison of scheduling mechanisms in scenario from Fig. 1, for k = 1

Fig. 3 .

 3 Fig. 3. LVWTF versus OPF in 2-wavelength receiver ring:A → B = 2 • (A → C)

TABLE I EFFICIENCY

 I MEASURE VALUE FOR LQF, LVWTF AND RAND

		E LQF	E LV W T F	E RAND
	k=0.2	0	0	0
	k=0.4	0	0.27	0
	k=0.6	0	0.65	0
	k=0.8	1.40	1.46	0.88
	k=0.9	5.77	2.27	4.74
	k=0.95	14.66	2.51	12.03
	k=1.0	74.41	0	76.33

ACKNOWLEDGMENT

The work described in this paper was carried out with the support of the French "Agence Nationale de la Recherche" in the framework of the ECOFRAME project (grant number 2006 TCOM-002-06), and of the BONE Network of Excellence ("Building the Future Optical Network in Europe"), funded by the European Commission through the 7th ICT-Framework Programme.