Hadrien Glaude

Cyrille Enderli

Olivier Pietquin

SPECTRAL LEARNING WITH PROPER PROBABILITIES FOR FINITE STATE AUTOMATON

Keywords: Baum-welch, learning automata, spectral learning, language models, non-negative matrix factorization

Probabilistic Finite Automaton (PFA), Probabilistic Finite State Transducers (PFST) and Hidden Markov Models (HMM) are widely used in Automatic Speech Recognition (ASR), Text-to-Speech (TTS) systems and Part Of Speech (POS) tagging for language modeling. Traditionally, unsupervised learning of these latent variable models is done by Expectation-Maximization (EM)-like algorithms, as the Baum-Welch algorithm. In a recent alternative line of work, learning algorithms based on spectral properties of some low order moments matrices or tensors were proposed. In comparison to EM, they are orders of magnitude faster and come with theoretical convergence guarantees. However, returned models are not ensured to compute proper distributions. They often return negative values that do not sum to one, limiting their applicability and preventing them to serve as an initialization to EM-like algorithms. In this paper, we propose a new spectral algorithm able to learn a large range of models constrained to return proper distributions. We assess its performances on synthetic problems from the PAutomaC challenge and real datasets extracted from Wikipedia. Experiments show that it outperforms previous spectral approaches as well as the Baum-Welch algorithm with random restarts, in addition to serve as an efficient initialization step to EM-like algorithms.

INTRODUCTION

Natural language processing (NLP) makes extensive use of automata and relies on training data to learn their parameters. The focus of this paper is a novel algorithm for learning parameters of a general family of automata, capturing most of the ones that are being used for NLP, i.e. Probabilistic (non-deterministic) Finite Automata (PFA). A PFA is a generative model for sequential data. Intuitively a PFA is similar to a Hidden Markov Model (HMM) in that it consists of a set of states, each of which when visited emits a symbol according to an emission probability distribution. Probabilistic Finite State Transducers (PFST) are also similar to PFA defined over pairs of input and output symbols. These models have been used, for example, to train POS tagging models [START_REF] Johnson | Why doesn't em find good hmm postaggers?[END_REF], Automatic Speech Recognition [START_REF] Lawrence R Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] and other applications [START_REF] Mohri | Finite-state transducers in language and speech processing[END_REF][START_REF] Merialdo | Tagging english text with a probabilistic model[END_REF]. Learning such automata parameters often rely on Expectation-Maximization (EM)-like algorithms trying to maximize the joint likelihood like Baum-Welch, Viterbi, Variational Bayes or Gibbs Sampling. However, these algorithms are iterative procedures and converge to local minima. In addition to being prone to get stuck into local optima, these algorithms are computationally expensive, to the point where obtaining good solutions for large models becomes intractable.

A recent alternative line of work consists in designing learning algorithms for latent variable models exploiting the so-called Method of Moments (MoM). The MoM leverages the fact that low order moments of distributions contain most of the distribution information and are typically easy to estimate. The MoM have several pros over iterative methods: it can be applied to a large variety of models [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF][START_REF] Bailly | A spectral approach for probabilistic grammatical inference on trees[END_REF][START_REF] Thon | Links between multiplicity automata, observable operator models and predictive state representations-a unified learning framework[END_REF][START_REF] Glaude | Subspace identification for predictive state representation by nuclear norm minimization[END_REF] and MoM-based algorithms are extremely fast as estimated moments can be computed in a linear time w.r.t. the number of samples. In addition, they are often consistent with theoretical guarantees in form of finite-sample bounds. This is particularly appealing because the resulting parameters can be used to initialize EM-like algorithms in a two-stage learning procedure, resulting in the best of both worlds: parameters efficiently computed by a MoM-based algorithm provide a good initialization for a maximum likelihood approach.

Beyond all the appealing traits of MoM-based algorithms, a well-known concern is the Negative Probability Problem (NPP) [START_REF] Thon | Links between multiplicity automata, observable operator models and predictive state representations-a unified learning framework[END_REF]. When learnt from samples, MoM estimators are not guaranteed to find a valid set of parameters. Although errors in the estimated parameters are bounded, the parameters themselves may lie outside the class of models that define proper distributions. NPP prevents the use of MoM-based algorithms to initialize a local search algorithm like EM. As mentioned in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF][START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF], this is a longstanding issue. Most of the recent attempts try to project the learnt model onto the space of valid model parameters [START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF]. While these heuristics produce a usable model, the resulting model may no longer be close to the true parameters. In this paper, we adopt a different approach by identifying a subclass of models realizing proper distributions and learnable through the MoM. Although this subclass does not encompass all the generative models for sequential data, it has an expressiveness large enough to closely approximate many of the distributions used in practice. In addition, learnt models from that subclass can serve as a good initialization to EM-like algorithms. The paper is organized as follows: in Section 2, we recall the definition of a PFA and develop the basic Spectral learning algorithm; in Section 3 we define a particular form of automaton capturing proper distributions and a polynomial-time learning algorithm that we call SepPFA (for Separable PFA); Section 4 analyses classes of automaton that can and cannot be learnt by SepPFA; finally, we assess the performance of SepPFA on synthetic problems and real large datasets that cannot be handled by traditional methods like EM.

BACKGROUND

Probabilistic Finite Automaton

PFA are graphical models constrained to represent distributions over sequences of symbols. Let Σ be a set of symbols, also called an al-phabet. We denote by Σ , the set of all finite words made of symbols of Σ, including the empty word ε. Words of length k form the set Σ k . Let u and v ∈ Σ , uv is the concatenation of the two words and uΣ is the set of finite words starting by u. We are interested in capturing a distribution over Σ . Let p be such a distribution, for a set of words S, we define p(S) = u∈S p(u), in particular we have that p(Σ ) = 1. In addition, for any word u we define p such that p(u) = p(uΣ ). Thus, p defines distributions over prefixes of fix length : ∀n, u∈Σ n p(u) = 1.

Some of these distributions can be modeled by graphical models called Probabilistic Finite Automaton (PFA).

Definition. A PFA is a tuple Σ, Q, {Ao} o∈Σ , α0, α∞ , where Σ is an alphabet and Q is a finite set of states. Matrices Ao ∈ IR + |Q|×|Q| contain the transition weights. The vectors α∞ ∈ IR + |Q| and α0 ∈ IR + |Q| contain respectively the terminal and initial weights. These weights should verify,

1 α0 = 1 α∞ + o∈Σ Ao1 = 1 (1)
A PFA realizes a distribution over Σ defined by

p(u) = p(o1 . . . o k ) = α 0 Auα0 = α 0 Ao 1 . . . Ao k α∞. ( 2
)
Because of the constraints defined in Equation ( 1), the weights belong to [0, 1] and can be viewed as probabilities over initial states, terminal states and transitions with symbol emission. For a word, we define a path as a sequence of states starting in an initial state, transiting from state to state, emitting symbols of the word in each state and exiting in a final state. The probability of a path is the product of the weights along the path including initial and final weights. Hence, the probability of a word is given by the sum of all paths probabilities, as written in Equation [START_REF] Lawrence R Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]. A path with a positive probability is called an accepting path. PFA define a particular kind of Multiplicity Automaton (MA). A MA is also a tuple Σ, Q, {Ao} o∈Σ , α0, α∞ but without any constraints on the weights, which can be negative and lose their probabilistic meaning. We call a formal power series a function that associates to any word of Σ a real. Thus, MA can realize formal power series that are not distributions. A MA that realizes a distribution is called a Stochastic MA (SMA). SMA are more general than PFA because they do not require the non-negativity of the weights. The Spectral algorithm presented in the next section relies on the Hankel matrix representation of a series to learn MA.

Spectral Learning

Let p : Σ → IR be a formal power series, we define H ∈ IR Σ ×Σ the bi-infinite Hankel matrix whose rows and columns are indexed by When p is a distribution, H contains occurrence probabilities that can be estimated from samples by empirical frequencies. Let for all o ∈ Σ, Ho ∈ IR Σ ×Σ , hS ∈ IR Σ and hP ∈ IR Σ be such that Ho(u, v) = p(uov), hS(u) = hP (u) = p(u). These vectors and matrices can be extracted from H. The Hankel representation of formal series lies at the heart of all MoM-based learning algorithms, because of the following fundamental theorem.

Σ such that H[u, v] = p(uv), H =        ε a
Theorem 1 (See [START_REF] Carlyle | Realizations by stochastic finite automata[END_REF]). Let p be a formal power series realized by a MA with n states, then rank(Hp) ≤ n. Conversely, if the Hankel matrix Hp of a formal power series p has a finite rank n, then p can be realized by a MA with exactly n states but not less.

For a MA with n states, observe that H[u, v] = (α 0 Au)(Avα∞). Let P ∈ K Σ ×n and S ∈ K n×Σ be matrices defined as follows,

P = ((α 0 Au) ) u∈Σ , S = (Avα∞)v∈Σ ,
then H = P S. Moreover, we have that,

Ho = P AoS, h S = α 0 S, hP = P α∞. (3)
So the MA parameters can be recovered by solving Equation ( 3). Hopefully, we do not need to consider the bi-infinite Hankel matrix to recover the underlying MA. Given a basis B = (P, S) of prefixes and suffixes, we denote by HB the sub-block of H. Similarly, HB,o is a sub-block of Ho. A basis B is complete if HB has the same rank than H. In [START_REF] Balle | Learning finite-state machines: algorithmic and statistical aspects[END_REF], the author shows that if B = (P, S) is a complete basis, by defining P over P, S over S and H, Ho over B, we can recover a MA using Equation ( 3). In the literature, several methods are used to build a complete basis from data. For example, one can use all prefixes and suffixes that appear in the training set. Once a basis is chosen, the Spectral algorithm first estimates the probabilities in ĤB by empirical counts. Then, it recovers a factorized form of ĤB = Û D V through a truncated Singular Value Decomposition (SVD). Finally, setting P = Û D and Ŝ = V , the algorithm solves Equation (3),

Âo = D-1 Û ĤB,o V , α 0 = ĥ S V , α∞ = D-1 Û ĥP .
In the experiments, following the advice of [START_REF] Shay B Cohen | Experiments with spectral learning of latent-variable pcfgs[END_REF], we normalized the feature-variance of the coefficients of the Hankel matrix by independently scaling each row and column by a factor cu = |S| /(#u + 5), where #u is the number of occurrences of u. In addition, depending on the problem, it can be better to work with other series derived from p. For example, the substring-based series p substring (u) = w,v∈Σ p(wuv) is related to p string . According to [START_REF] Balle | Learning finite-state machines: algorithmic and statistical aspects[END_REF], if p string is realized by a SMA, then p substring is too. In addition, [START_REF] Balle | Learning finite-state machines: algorithmic and statistical aspects[END_REF] provides an explicit conversion between string-based SMA and substring-based SMA preserving the number of states. For all MoM-based algorithms compared on the PAutomaC, we used the series leading to the best results for each problem.

Although the Spectral algorithm can return a MA arbitrary close to a SMA or a PFA that realizes the target distribution, it does not ensure that the returned MA will be a SMA or a PFA. This causes the so-called Negative Probability Problem and prevents to use the learned MA as an initialization to EM. Indeed, when the Spectral algorithm computes a low rank approximation of the sampled Hankel matrix, it does not preserve the Hankel structure. In [START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF], authors proposed to compute a Structured Low-Rank Approximation (SLRA). However, the problem is strongly non-convex and thus cannot be applied directly on large scale problems.

RESIDUAL REPRESENTATION

This section defines a particular representation of Multiplicity Automaton called residual realizing a formal power series p. As in the Spectral algorithm, our method assumes that a complete basis B = (P, S) is given as input. Let P be a subset of P, such that rank(HB) = rank(H B ) where B = (P , S). In addition, let's assume there exists for all u ∈ P and all o ∈ Σ, two sets of variables {a v u,o } v∈P and {a v ε } v∈P verifying, for all w ∈ Σ ,

p(uow) p(uo) = v∈P a v u,o p(vw) p(v) and p(w) = v∈P a v ε p(vw) p(v) (4)
In IR-MA, states are associated to elements of P and the weights are defined as follows,

α 0 = (a u ε ) u∈P (5) α∞ = p(u) p(u) u∈P (6) ∀u, v ∈ P , Ao[u, v] = a v u,o p(uo) p(u) . ( 7
)
The next proposition gives a necessary condition for a formal power series p to admit a residual representation. Proposition 2. If there exists P and {a v u,o } v∈P , {a v ε } v∈P verifying Equation (4) for all w ∈ Σ , then for any word w = o1...on, we have,

Awα∞ = Ao 1 ...Ao n α∞ = p(uw) p(u) u∈P , (8)
and so, p(w) = α 0 Awα∞.

Proof. The proof is done by induction on |w|. For the empty word, Equation ( 8) follows from the definition of α∞. Now, for w = o1...on, the induction hypothesis gives,

Ao 2 ...Ao n α∞ = p(uo2...on) p(u) u∈P , So we have, Awα∞ = v∈P p(vo2...on) p(v) × a v u,o 1 p(uo) p(u) u∈P = p(uo) p(u) v∈P p(vo2...on) p(v) a v u,o 1 u∈P = p(uo1o2...on)) p(u) u∈P = p(uw) p(u) u∈P ,
where the last step uses Equation (4), with w = o2...on. In addition, Equation (4) gives,

α 0 Awα∞ = u∈P a u ε p(uw) p(u) = p(w)
Note that, for a MA with k states, the infinite set of linear equations given by Equation ( 4) for all w ∈ Σ forms a system of rank k and admits a solution which can be found by solving the set of Equation ( 4) for all w ∈ S, if B is complete.

There exists many methods [START_REF] Thon | Links between multiplicity automata, observable operator models and predictive state representations-a unified learning framework[END_REF][START_REF] Gybels | Some improvements of the spectral learning approach for probabilistic grammatical inference[END_REF] to compute a minimal subset P such that B is complete. Once B is found, one has to solve Equation ( 4) by linear regression using empirical estimations p of p, and then compute the weights using p and the a v u,o s. Recalling that our goal is to learn proper distributions, one would like to add constraints to ensure that the residual form MA learned by regression realizes a proper distribution. Unfortunately, this would require two things : 1) the non-negativity of series for any word (∀ ∈ Σ , p(u) ≥ 0), 2) the convergence of the series to one ( u∈Σ p(u) = 1). Although 2) can be checked in polynomial time, 1) requires adding an infinite set of constraints during the linear regression step. That is why, in general, verifying if a MA is a SMA is undecidable [START_REF] Franc | On rational stochastic languages[END_REF]. So, we restrict ourselves to the learning of PFA in their residual form.

Learning Residual PFA

In this section, we detail how to build a residual representation of a PFA that realizes a distribution (under some conditions on this distribution). As in the previous section, we start with a basis B = (P, S). In contrast to SMA, for PFA the non-negativity of the series is ensured by the non-negativity of the initial, final and transition weights. In addition, the convergence of the series is shown in [START_REF] Franc | On rational stochastic languages[END_REF] to be equivalent to the finite set of linear constraints given in Equation ( 1). Thus, constraints will be added to find both the a v u,o s and a minimal subset P of P, that ensures the non-negativity and the convergence of the series.

Firstly, the convergence of a series is shown in [START_REF] Franc | On rational stochastic languages[END_REF] to be equivalent to the finite set of linear constraints given in Equation ( 1). Thus, in its residual form, a PFA must verifies for all u ∈ P that p(u) p(u)

+ v∈P o∈Σ a v u,o p(uo) p(u) = 1 and v∈P a v ε = 1. (9)
Note that, when working with the true distribution p, Equation ( 9) is always satisfied, if Equation ( 4) holds for all w ∈ Σ . Indeed, we have

p(u) + v∈P o∈Σ a v u,o p(uo) = p(u) + o∈Σ p(uo) = p(u),
because, summing Equation (4) over w ∈ Σ gives v∈P a v u,o = 1 and v∈P a v ε = 1. However, if p and p are estimated, to ensure the convergence of the series, the constraints defined in Equation ( 9) has to be added when solving Equation (4).

Secondly, as the non-negativity of the weights in the residual form of a PFA is equivalent to the non-negativity of the a v u,o s and a v ε s, we must find a subset P verifying for all w ∈ Σ Equation ( 4) where a v u,o ≥ 0 and a v ε ≥ 0. This constraint can be rewritten in a vector form where pu is the infinite vector p(uw) p(u) w∈Σ where coni denote the conical hull of a set. Finding a subset P like that is not trivial in the general case. As an approximation, we propose to assume the following condition :

Condition 1 (k-Separability). There are k prefixes u1, ..., u k of P such that coni{pu | u ∈ Σ } = coni{pu 1 , ...pu k }, and

coni{p S u | u ∈ P} = coni{p S u 1 , ...p S u k },
where, p S u is the restriction of pu on the basis S.

Section 4 details the models satisfying this condition. Now, assuming the k-separability condition, let P be {u1, ..., u k }, then Equation ( 4) holds trivially for all w ∈ Σ .

Algorithm 1 The SepPFA algorithm

Input: A basis B satisfying the k-separability condition. Output: A PFA with a residual form Σ, P , {Ao} o∈Σ , α0, α∞ .

1: Build empirical estimates of p and p. We now detail Algorithm 1. In the second step, finding the k edges can be cast as a well-known problem: near-separable Non-negative Matrix Factorization (NMF). In the literature, many algorithms for NMF have been proposed. Although in its general form NMF is NP-Hard and ill-posed [START_REF] Gillis | The Why and How of Nonnegative Matrix Factorization[END_REF], in the near-separable case the solution is unique and can be found in a time proportional to O(k |S| |P|). State-of-the-art algorithms for near-separable NMF comes with convergence guarantees and robustness analysis, making the near-separable assumption appealing. For Algorithm 1, we used the Successive Projection Algorithm (SPA) with a post-processing step described and analyzed in [START_REF] Gillis | Semidefinite programming based preconditioning for more robust nearseparable nonnegative matrix factorization[END_REF]. As rescaling the p S u i s does not change the solution, we multiplied them by p(ui) in order to normalize the estimation noise, before applying the SPA. Indeed, p(ui) is an approximation of the inverse of the variance.

In the fourth and fifth steps of Algorithm 1 the minimization problem is actually an instance of a quadratic optimization problem under linear constraints, that can be solved in polynomial time using a solver like MOSEK. As the problem is strongly convex, it admits a unique stationary point to which the solver converges.

Finally, if the target distribution p and the given basis B satisfy the k-separability assumption, Algorithm 1 returns a PFA with k states in a residual form that realizes p. This automaton can then serve as an initialization to a Baum-Welch algorithm.

THE SEPARABILITY CONDITION

After defining the family of SPA that can be estimated by Algorithm 1 given a suitable basis, we show that k-separable PFA encompass a large variety of models, including the Residual PFA (RPFA) defined in [START_REF] Esposito | Learning probabilistic residual finite state automata[END_REF]. We start with few definitions.

Definition. A k-separable PFA is a PFA such that, there exists k prefixes u1, ..., u k satisfying

coni{pu | u ∈ Σ } = coni{pu 1 , ...pu k }.
Thus, if p is realized by a k-separable PFA M and a suitable basis B is given, Algorithm 1 launched with parameters k and B is able to recover the residual form of M from the true values of p. Next, we show that RPFA with k states are k-separable. First, we need few definitions. Let M = Σ, Q, {Ao} o∈Σ , α0, α∞ be a PFA, for each state q ∈ Q we define Mq = Σ, Q, {Ao} o∈Σ , β 0 , α∞ be the PFA such that β 0 = 1q. We denote by pq the distribution realized by Mq.

Definition. Let M = Σ, Q, {Ao} o∈Σ , α0, α∞ be a PFA, then M is a Residual PFA (RPFA) iff ∀q ∈ Q, ∃uq ∈ Σ , ∀w ∈ Σ pq(w) = p (uqw) p(uq) .
Theorem 3. A RPFA with k states is k-separable.

Proof. Let M = Σ, Q, {Ao} o∈Σ , α0, α∞ be a RPFA. Let qw be a vector such that qw[q] = pq(w) = p(uq w) p(uq ) . We denote by bq the q-th coefficient of the vector α 0 Au/p(u). Note that we have bq ≥ 0. By definition of RPFA, we have for all u, w ∈ Σ that,

p(uw) p(u) = 1 p(u) α 0 Auqw = q∈Q bq p (uqw) p(uq)
In terms of vectors, we have pu = q∈Q bqpu q with bq ≥ 0

The expressiveness of RPFA has been thoroughly studied in [START_REF] Esposito | Learning probabilistic residual finite state automata[END_REF]. Authors have shown that all distributions realized by a PDFA can be realized by a RPFA but there are some distributions a RPFA can realized that cannot be realized by PDFA. Thus, RPFA strictly generalize PDFA. In [START_REF] Pfau | Probabilistic deterministic infinite automata[END_REF], authors have shown that PDFA strictly generalize n-th order Markov chains for any n. On the bad side, PFA are strictly more general than RPFA [START_REF] Esposito | Learning probabilistic residual finite state automata[END_REF]. In [START_REF] Dupont | Links between probabilistic automata and hidden markov models: probability distributions, learning models and induction algorithms[END_REF], authors gave explicit conversions between HMMs and PFA, showing that a PFA can represent an HMM with the same number of states and a HMM can represent a PFA using an increased number of states. Thus, PFA and HMM have the same expressiveness, although PFAs seems slightly more compact.

EXPERIMENTAL RESULTS

PAutomaC Challenge

The Probabilistic Automata learning Competition (PAutomaC) deals with the problem of learning probabilistic distributions from strings drawn from finite-state automata. From the 48 problems available, we have selected the same twelve problems than in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF], to provide a fair comparison with other algorithms. The generating model can be of three kinds: PFA, HMMs or PDFA. Four models have been selected from each class. A detailed description of each problem can be found in [START_REF] Verwer | Results of the pautomac probabilistic automaton learning competition[END_REF]. Table 1 compares SepPFA, BW-SepPFA (Baum-Welch initialized with SepPFA) and BW (with random restarts initialization) to best results (denoted by MoM) among the following MoM-based algorithms : CO [START_REF] Balle | Local loss optimization in operator models: A new insight into spectral learning[END_REF] using strings statistics, Tensor [START_REF] Anandkumar | Tensor decompositions for learning latent variable models[END_REF] using strings statistics and Spectral using strings and substrings statistics. A description and comparison of these algorithms can be found in [START_REF] Balle | Methods of moments for learning stochastic languages: Unified presentation and empirical comparison[END_REF]. For SepPFA we used 500 prefixes and suffixes. The quality of a model can be measured by the quality of the probability distribution it realizes. The objective is to learn a MA realizing a series p close to the distribution p , which generated the training set T . The quality of p is measured by the perplexity that corresponds to the average number of bits needed to represent a word using the optimal code given by p . A grid search was performed to find the optimal rank for each of the performance metric. For each problem, the best between BW-SepPFA and SepPFA is indicated by a bold number and the best score between SepPFA and the best of other MoM-based algorithms is underlined. The score and the size of the true model is reported for comparison. First, SepPFA achieves the better score than MoM on the majority of the problems showing that RPFA can provide very good approximations of more complex models. When using BW on this dataset, without any stopping conditions, lot of over-fitting was observed. This explains why few iterations of BW degrades the performances of SepPFA on problem 42 and 43. The results in Table 1 have been obtained by stopping the BW algorithm when it does not improve the perplexity anymore. BW-SepPFA in comparison to BW achieves a lowest perplexity for 8 problems among 12. For the other 4 problems, the perplexity of both BW-SepPFA and BW are very close to the true one. This clearly argues for the use of Algorithm 1 as an initialization step to BW.

Perplexity(M) = 2 -u∈T p (u) log(p M (u))

Wikipedia

We evaluate Algorithm 1, on a large corpus of real-world Wikipedia text treated as chunked of 250 characters randomly extracted from the 2GB corpus used in [START_REF] Sutskever | Generating text with recurrent neural networks[END_REF]. Each character stands for an observation. There are 85 different characters in the data set. For the training phase, we used a small database of 500 chunks and a medium one of 50000 characters. For testing, an independent data set of 5000 chunks is used. We evaluate three algorithms : Baum-Welch (BW) with 5 random restarts of 4 iterations (then the best run is continued for 20 iterations), SepPFA (Sep) and Baum-Welch initialized with SepPFA (BW-Sep) with 20 iterations. Performance is measured by the average likelihood of the next observation given the past. More precisely, on chunk of 250 characters, the 50 firsts serve to initialize the belief, then the next observation is predicted based on the characters observed so far. After predicting the next observation, this is added to the past and a prediction is made for the next one until the end of the chunk is reached. The likelihood of all these predictions is then averaged over observations and chunks. Many conclusions can be drawn from Figures 1 and2. First, on the two datasets, Baum-Welch initialized using SepPFA performs best. Although BW and BW-SepPFA achieve almost the same likelihood on the joint distribution in our experiments, we can see on figs. 1 and 2 that BW achieves a very small likelihood on the conditional distribution. This means that the probabilities of rare events are incorrectly estimated by BW, where MoM-based algorithms are much more efficient. Actually, BW tends to over-fit the distribution, but as RPFA defines a simpler class of model it brings generalization properties and so is a good model to initialized BW. Comparing SepPFA and Spectral, SepPFA seems to have better generalization properties on small datasets, whereas Spectral is able to learn more complex models by better fitting the target distribution when more data is available. This can be explained by the fact that SepPFA is able to learn only a particular constrained subclass of MA. About the size of the basis, it denotes the number of prefixes or suffixes we used. In our experiments, we used the same set for prefixes and suffixes. This set contains the most frequents substrings of length up to 3 in the training set. SepPFA behaves better on small basis, whereas Spectral uses very large basis to capture the maximum of information. Indeed, the larger the basis is, the more infrequent substrings will be added and so considered as potential edges of the conical hull. When the vectors pu are not well estimated, they prevent the near-separable NMF algorithm identifying the right edges, leading to potential high errors. Besides the improved performances, being able to work with small basis is a good gain for computation speed and scalability. As expected, applying BW to the output of SepFPA has always improved the solution. This again advocates for the combination of SepPFA and BW. Finally, using two datasets of different sizes shows how BW cannot handle large data sets. The running time of BW between the two datasets has changed from few minutes to at least a day, whereas the one of SepPFA has only slightly increased of one or two seconds, i.e. just the time needed to compute the empirical estimates on a larger dataset.

CONCLUSIONS

In this paper, we proposed a new algorithm based on a near-separable NMF and constraint quadratic optimization that can learn in polynomial time a separable PFA from the distribution p it realizes. In addition, even if p is not realized by a separable PFA or is empirically estimated, our algorithm returns a separable PFA that realizes a proper distribution. Then, we empirically demonstrated both, its good performances in comparison to other MoM-based algorithms, and, its scalability in comparison to BW. Finally, experiments have shown that initializing EM with SepPFA has almost always improved the performances of both EM and SepPFA. For future work, one could try to establish the consistency of Algorithm 1 and some finite sample bounds based on the robustness analysis of the near-separable NMF. In addition, extending Algorithm 1 to handle controlled processes such that (Partially Observable) Markov Decision Processes would allows to efficiently learn dialogue systems as proposed in [START_REF] Lemon | Data-Driven Methods for Adaptive Spoken Dialogue Systems: Computational Learning for Conversational Interfaces[END_REF].

,

  ∀u ∈ P , o ∈ Σ, puo = v∈P a v u,o pv with a v u,o ≥ 0 and pε = v∈P a v ε pv with a v u,o ≥ 0 . meaning that ∀u ∈ P , o ∈ Σ, puo ∈ coni{pu | u ∈ P } and pε ∈ coni{pu | u ∈ P } .

Fig. 1 .

 1 Fig. 1. Performances on the small Wikipedia training set.

Fig. 2 .

 2 Fig.2. Performances on the medium Wikipedia training set. For BW and BW-SepPFA, the computation was too long to be evaluated for all model sizes.

Table 1 .

 1 Comparison with other algorithms for perplexity. Model sizes are listed in parentheses.

	ID	BW-Sep	Sep	BW	MoM True model
	HMM 1	38.68(5)	68.91 500.10	44.77	29.90(63)
	14 116.92(15) 125.09 116.84 128.53	116.79(15)
	33	32.11(20)	38.44	32.14	49.22	31.87(13)
	45	25.92(5)	28.73 107.75	31.87	24.04(14)
	PDFA 6	67.09(40)	83.71	67.32	95.12	66.98(19)
	7	51.30(15)	66.51	51.27	62.74	51.22(12)
	27	59.39(20)	74.48	94.40 102.85	42.43(19)
	42	20.47(35)	18.13 168.52	23.91	16.00(6)
	PFA 29	24.71(30) 275.21	25.09	34.57	24.03(36)
	39	14.30(5)	18.96	10.43	11.24	10.00(6)
	43	36.08(5)	34.93 461.23	36.61	32.64(67)
	46	12.07(30) 105.62	12.02	25.28	11.98(19)


Min Zhu

Shilin Xiao
email: slxiao@sjtu.edu.cn

Wei Guo

He Chen

Anne Wei

Yaohui Jin

Weisheng Hu

Benoit Geller

Overlay of Multicast Service in WDM-PON Based on Dynamic Wavelength Reflection Scheme

Keywords: WDM-PON, multicast, dynamic wavelength reflection, fiber bragg grating (FBG)

We propose a novel scheme to overlay multicast service over a wavelength-division-multiplexed (WDM) passive optical network (PON) based on dynamic wavelength reflection scheme. At the optical line terminal (OLT), for each WDM channel, a dynamic wavelength reflector is used to selectively bypass the corresponding optical carrier to avoid multicast data superimposition, while the remaining optical carriers go through the dynamic wavelength reflector and simultaneously carry the multicast data. The proposed scheme is experimentally demonstrated with 5-Gb/s downstream unicast, multicast data and 625-Mb/s upstream re-modulated data.

INTRODUCTION

WDM-PON is regarded as a promising solution for next generation broadband access networks because of its almost unlimited bandwidth, network security, simple management and upgradeability. To deploy multicast service over WDM-PON, some efforts have already been reported [START_REF] Deng | A WDM Passive Optical Network with Centralized Light Sources and Multicast Overlay[END_REF][START_REF] Cai | Video-Service-Overlaid Wavelength-Division-Multiplexed Passive Optical Network[END_REF][START_REF]Selective multicasting of Digital Video Signals over a WDM Passive Optical Network[END_REF][START_REF] Moon | Overlay of multicasting Signal in a WDM-PON[END_REF][START_REF] Ohm | Quaternary Optical ASK-DPSK and Receivers with Direct Detection[END_REF] , and fall into three categories: 1) In-band transmission [START_REF] Deng | A WDM Passive Optical Network with Centralized Light Sources and Multicast Overlay[END_REF][START_REF] Cai | Video-Service-Overlaid Wavelength-Division-Multiplexed Passive Optical Network[END_REF] , the multicast service shares the downstream wavelength bandwidth with unicast data service via time-division multiplexing scheme or orthogonal modulation format, which requires complicated timing or synchronization control; 2) Out-of-band transmission [START_REF]Selective multicasting of Digital Video Signals over a WDM Passive Optical Network[END_REF][START_REF] Moon | Overlay of multicasting Signal in a WDM-PON[END_REF] , multicast service is overlaid on the same wavelength channel by means of sub-carrier multiplexing (SCM), which requires high-frequency electronic components at both the transmitter and receiver sides; 3) Out-ofwavelength transmission [START_REF] Ohm | Quaternary Optical ASK-DPSK and Receivers with Direct Detection[END_REF] , a additional specific wavelength channel is used to transmit the multicast service with simple optical power splitter and the arrayed waveguide grating (AWG). However, to accommodate the dynamic change of members in the multicast group, the dynamic deployment of broadcast/multicast service over WDM-PON needs to be further studied.

In this paper, a simple scheme to overlay multicast service is proposed. It uses a set of dynamic wavelength reflector added on each wavelength channel. When the dynamic wavelength reflector is in the bar (straight-through) state, the multicast data is orthogonal-overlaid with downstream unicast data on the corresponding optical carrier. The multicast data could not be superimposed when the dynamic wavelength reflector is in the reflection state to bypass the corresponding downstream optical carrier. By dynamically adjusting the state of the dynamic wavelength reflector, the multicast data can be enabled on a unicast wavelength channel to achieve selective multicast and thus realizing flexible network functions. We experimentally demonstrate the proposed WDM-PON with multicast overlay at 5-Gb/s and investigate its feasibility and system performance. The upstream data re-modulation based on the downstream nonreturn-to-zero (NRZ) format with a low extinction ratio (ER) is also demonstrated.

PROPOSED SCHEME AND PRINCIPLE

The proposed WDM-PON architecture with multicast overlay is shown in Fig. 1. In the optical line terminal (OLT), the downstream carrier of each wavelength channel is generated by a continuous-wave (CW) laser, and modulated by an optical intensity modulators (IMs), to generate NRZ unicast data, with an extinction ratio (ER) of around 3dB. The lower ER makes it easier that the successful superimposition of the differential phase-shift keying (DPSK) multicast data and upstream data re-modulation 5 . To realize flexible multicast, the dynamic wavelength reflector can be added in between IM and arrayed waveguide grating (AWG) on each wavelength channel. In the schematic diagram as shown in the Fig. 1, the dynamic wavelength reflector consists of an optical circulator and a controllable Fiber Bragg Grating (FBG).

When an optical network unit (ONU) on a wavelength makes a request for the multicast service, the controllable FBG in the corresponding wavelength channel is adjusted via the control circuit to deviate from center wavelength and to allow unicast data carrier go straight to the AWG. Thus multicast signal can be simultaneously modulated on the downstream unicast carrier by an optical phase modulator (PM). To disable the multicast signal for an ONU, the FBG is dynamically tuned and reflects the corresponding wavelength to circulator and another AWG, so as to bypass multicast signal modulation. Hence the multicast signal can be selectively and dynamically superimposed onto downstream unicast signal on arbitrary wavelength channels. Then these two multiplexed signals I and II are combined by coupler and carried to remote node (RN) and finally routed to individual ONU. It should be noted that all the multicast control intelligence is deployed at OLT by adjusting a set of the dynamic wavelength reflector (FBG and circulator), which is transparent to all ONUs and can be a cost-effective solution. After transmission over fiber, a 2×16 AWG is used to de-multiplex the optical signals and route them to individual ONUs. In ONU, one part of downstream power is tapped off by an optical splitter for NRZ downstream unicast data and differential-phase-shift-keying (DPSK) multicast data, respectively. For users who subscribe to multicast service, the optical carrier not only carries the downstream unicast data, but also contains the multicast data information. The remaining part is re-modulated by upstream data with higher ER and sent back to the OLT. In this way, colorless ONU without the need of light source could also be implemented. It is an attractive feature for PON deployment.

EXPERIMENTAL SETUP AND RESULTS

To verify the feasibility of the proposed scheme, we demonstrate downstream unicast/multicast transmission in two different wavelength channels and upstream data re-modulation in Fig. 2.

In the OLT, continues wave (CW) in the 1551.95 nm wavelength is fed into a Mach-Zehnder modulator (MZM) which is biased at the quadrature point and is driven by a 5-Gb/s 2 31 -1 pseudorandom binary sequence (PRBS) to generate downstream unicast NRZ signal, with 3dB ER. After an EDFA and a tunable band-pass filter (BPF) with a bandwidth of 0.4 nm, the unicast NRZ signal is orthogonally phase-modulated by a PM to superimpose multicast data. To obtain better signal performance, downstream unicast data and multicast data signals need to be bit synchronized, which is adjusted by an electrical phase shifter.

Another CW at 1548.84nm is just modulated by unicast data with 5-Gb/s complementary PRBS and then is bypassed by an optical circulator and FBG with a 3-dB bandwidth of 0.15 nm and a 90% reflection ratio. Thus, the downstream optical carrier just contains unicast data. Then, both of the above wavelengths are combined and fed into two pieces of 12.5km single mode fiber (SMF) to reach respective ONU. An AWG located in the Remote Node (RN) links the two pieces of SMF. Fig. 2 Experimental setup of the proposed scheme After transmission over 25-km SMF, at the ONU, an optical splitter is used to tap optical carrier power into two optical paths. The part of the received downstream signal power is again split into two parts. One half is detected for NRZ downstream unicast data and another half is demodulated by a delay interferometer (DI) for DPSK multicast data. The remaining part is re-modulated by 625-Mb/s upstream data with a higher ER of ~10dB and sent back to the OLT. The bit-error-rate (BER) curves and eye diagrams of the downstream unicast NRZ signal and multicast DPSK signal are provided in Fig. 3 a-b). The power penalties of downstream unicast signal caused by the fiber transmission are ~0.8 dB with multicast data and ~1.3dB without multicast data, respectively. The power penalty is ~0.5 dB for the DPSK multicast signal.

The four eye diagrams are observed in Fig. 3 b). The inset figure 1) shows the eye diagram of the demodulated multicast DPSK signal with downstream unicast NRZ, which has three folds of up-eyelids. The inset figures 2) and 3) display the eye diagrams for the demodulated NRZ with DPSK signal and the demodulated NRZ without DPSK signal, respectively. We observe that the performance degradation of the above eye diagrams induced by the superimposition of multicast DPSK signal is represented by the thicker eyelid. The eye diagram for the downstream unicast NRZ demodulated by DI without DPSK is shown in the inset figure 4).

1)

2)

3)

4)

CONCLUSIONS

We have proposed and experimentally demonstrated a novel scheme to simultaneously deliver downstream unicast and multicast signal based on dynamic wavelength reflection scheme over WDM-PON architecture. In each WDM channel, a dynamic wavelength reflector is used to selectively bypass the corresponding optical carrier to avoid multicast data superimposition, while the remaining optical carriers go through the dynamic wavelength reflector and simultaneously carry the multicast data. The multicast function can be quickly and dynamically reconfigured with simple and centralized control in the OLT and is transparent for all ONUs. The technique just adds a little complexity to the existing WDM-PON structure, and its feasibility is demonstrated with 5-Gb/s downstream unicast, multicast signal and 625-Mb/s upstream re-modulated data. By employing re-modulation technique, no light sources and colored components are needed in the ONUs, which effectively reduce the cost and complexity.

Fig. 1

 1 Fig. 1 Schematic diagram of the proposed WDM-PON architecture with multicast overlay

Fig. 3 (

 3 Fig. 3 (a) BER curves for the downstream NRZ/DPSK

Fig. 3 (

 3 Fig. 3 (b) Eye diagrams of 1) the downstream demodulated DPSK with NRZ, 2) the demodulated NRZ with DPSK, 3) the NRZ demodulated by DI without DPSK, 4) the demodulated NRZ without DPSK

Acknowledgements: The work was jointly supported by the National Nature Science Fund of China (No. 60972032 and No.60632010) and the National "863" Hi-tech Project of China (No. 2006AA01Z251 and No. 2007AA01Z271).


S Bay

C Herzet

J M Brossier

J P Barbot

B Geller

Analytic and Asymptotic Analysis of Bayesian Cramér-Rao Bound for Dynamical Phase Offset Estimation

In this paper, we present a closed-form expression of a Bayesian Cramér-Rao lower bound for the estimation of a dynamical phase offset in a non-data-aided BPSK transmitting context. This kind of bound is derived considering two different scenarios: a first expression is obtained in an off-line context and then, a second expression in an on-line context logically follows. The SNR-asymptotic expressions of this bound drive us to introduce a new asymptotic bound, namely the Asymptotic Bayesian Cramér-Rao Bound. This bound is close to the classical Bayesian bound but is easier to evaluate.

I. INTRODUCTION

During the last decade, synchronization has become one of the most challenging task of digital receivers. On the one hand, since the advent of very powerful error correcting codes (see e.g., [START_REF] Gallager | Low-Density Parity-Check Codes[END_REF]- [START_REF] Pyndiah | Near-optimum decoding of product codes: block turbo codes[END_REF]), the synchronizers have to properly operate at Signal-to-Noise Ratios (SNRs) lower than ever before. On the other hand, the high data-rate requirements of future communication standards imply rapidly varying synchronization parameters at the receiver. In order to cope with such stressful environments, Maximum-Likelihood (ML) synchronization performance are often a desirable feature. Unfortunately, the ML estimator is often unfeasible in practice. Hence, sub-optimal synchronizers have been proposed, see e.g., [START_REF] Meyr | Digital Communication Receivers: Synchronization, Channel Estimation, and Signal Processing[END_REF].

In order to assess the performance of such synchronizers, lower bounds on mean-squared estimation error are needed. In particular, the family of Cramér-Rao Bounds (CRBs) [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF], [START_REF] Cramér | Mathematical Methods of Statistics[END_REF] has been shown to give tight estimation lower bounds in a number of practical scenarios, (see e.g., [START_REF] Rife | Single tone parameter estimation from discrete-time observations[END_REF]- [START_REF] Rice | Cramér-Rao Lower Bounds for QAM Phase and Frequency Estimation[END_REF] for CRB contributions associated to synchronization problems). Rife et al. [START_REF] Rife | Single tone parameter estimation from discrete-time observations[END_REF] and Cowley [START_REF] Cowley | Phase and frequency estimation for PSK packets: bounds and algorithms[END_REF] have derived CRB closedform expressions for constant phase-offset estimation in the socalled Data-Aided (DA) and Non-Data-Aided (NDA) scenarios, respectively. In [START_REF] Noels | The Cramér-Rao bound for phase estimation from coded linearly modulated signals[END_REF], the authors proposed a semi-analytical method enabling to efficiently evaluate CRBs in Code-Aided (CA) scenarios. A Modified CRB (MCRB), easier to evaluate than the standard CRB, has been introduced in [START_REF] D'andrea | The modified Cramér-Rao bound and its application to synchronization problems[END_REF], [START_REF] Gini | The modified Cramér-Rao bound in vector parameter estimation[END_REF]. The MCRB is usually looser than the CRB but it is equivalent to the CRB at high SNR [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF]. More recently, the problem of deriving CRBs suited to time-varying parameters has been This work is partially supported by the European Community contract NO. 507325, NEWCOM. addressed throughout the Bayesian context. In [START_REF] Tichavsky | Posterior Cramér-Rao bounds for discrete-time nonlinear filtering[END_REF], the authors propose a general framework for deriving analytical expression of on-line or also called filtering CRBs. In [START_REF] Dauwels | A numerical method to compute Cramér-Raotype bounds for challenging estimation problems[END_REF], the authors developed a numerical graph-based algorithm to evaluate the CRB in time-varying scenarios.

In this contribution, we address the open problem of deriving an analytical expression of the off-line or also called smoothing Bayesian CRB (BCRB) for time-varying phase estimation in NDA scenarios.

Explicit expressions of the BCRB and its modified version are provided in an off-line scenario. This bound helps us to evaluate and also to predict the estimator performance without any particular assumption or simulation restriction. Moreover, an on-line interpretation follows. The asymptotic cases at low and high SNR are presented. In particular it is shown that the Modified BCRB (MBCRB) is equal to the derived asymptote at high SNR. A new asymptote-based bound is then introduced simplifying the general expressions of the BCRB.

This paper is organized as follows. In Section II we set the system model. After recalling in Section III the different kinds of Cramér-Rao bounds, we derive in Section IV the off-line BCRB. In Section V the asymptotic cases are considered and we present the asymptote-based lower bound. Next, the closed-form expression of the BCRB in the on-line context is given (Section VI). Finally, the different results are illustrated and interpreted in Section VII.

The notational convention adopted is as follows: italic indicates a scalar quantity, as in a; boldface indicates a vector quantity, as in a and capital boldface indicates a matrix quantity as in A. The (m, n) th entry of matrix A is denoted [A] m,n . The matrix transpose operator is indicated by a superscript T as in A T . |.| denotes the modulus of a scalar quantity or the determinant of a matrix quantity. a n m represents the vector

[a m • • • a n ]
T , where m and n are positive integers (m < n). {a} and {a} are respectively the real and the imaginary part of a. E xy [.] denotes the expectation over x and y. k is equal to the integer part of k. δ k,l is the Kronecker symbol. ∇ θ and ∆ θ ψ represent the first and second-order partial derivatives operator i.e.,

∇ θ = ∂ ∂θ 1 • • • ∂ ∂θ K T and ∆ θ ψ = ∇ ψ ∇ T θ . The error function is defined as erf (x) = 2 √ π
x 0 e -t 2 dt. The actual value of an unknown parameter is indicated by the superscript as in θ .

II. SYSTEM MODEL

We consider the transmission of a BPSK modulated sequence

a = [a 1 • • • a K ]
T over an Additive White Gaussian Noise (AWGN) channel affected by some carrier phase offsets stacked in a vector θ

= [θ 1 • • • θ K ]
T . Assuming that the received signal has been ideally filtered and sampled at the optimum sampling instant, the discrete-time baseband signal

y = [y 1 • • • y K ]
T is given by

y k = a k e jθ k + n k with k = 1 . . . K, ( 1
)
where a k is the k th unknown transmitted BPSK symbol (a k = ±1) and n k is a zero-mean circular Gaussian noise with known variance σ 2 n . We suppose that the system operates in a NDA synchronization mode, i.e., the transmitted symbols are independent and identically distributed (i.i.d.) with p(a k = a) = 1 2 where a = ±1. We assume that the energy per symbol is normalized, so that the noise variance is related to the E s /N 0 -ratio as

E s N 0 = 1 σ 2 n .
We consider the case of a Wiener phase-offset evolution, i.e.,

θ k = θ k-1 + w k , ( 2
)
where w k is an i.i.d. zero-mean Gaussian noise with known variance σ 2 w . This model is commonly used [START_REF] Mcneill | Jitter in ring oscillators[END_REF], [START_REF] Demir | Phase noise in oscillators: a unifying theory and numerical methods for characterization[END_REF] to describe the behavior of practical oscillators for which the frequency is randomly perturbated. As the phase is the integral of the frequency, the phase variance increases with time index k. This phenomenon is thus non stationary and one readily obtains that var (θ k ) = var (θ k-1 ) + σ 2 w . In order to recover the transmitted symbols, practical receivers have first to estimate the carrier phase offsets corrupting the received observations and perturbing the whole data recovery. Phase estimation can actually be considered following two main scenarios: i) Off-line synchronization: the receiver waits until the whole observation frame, i.e., y = [y

= [y 1 y 2 • • • y k ] T .
In the sequel, the Cramér-Rao bound evaluation will be considered within the context of both the off-line and the online scenarios.

III. CRAM ÉR-RAO BOUNDS: SOME USEFUL EXPRESSIONS

A natural question which arises when designing estimators is the ultimate accuracy that can be achieved in the estimation operation. The family of CRBs (partially) answers this question by providing a lower bound on the Mean Square Error (MSE) achievable by any unbiased estimator. Different kinds of CRBs can be considered. In the following, we briefly recall general expressions of these bounds.

A. The Standard Cramér-Rao Bound and its Modified Version

Let θ(y) denote an unbiased estimator of θ. Then the estimator θ(y) satisfies the following inequality

E y|θ=θ ( θ(y) -θ ) ( θ(y) -θ ) T ≥ SCRB (θ ), (3)
where SCRB is the Standard CRB (SCRB) and θ is the actual parameter value. The SCRB is equal to the inverse of the Fisher Information Matrix (FIM) F(θ ) defined by [START_REF] Cramér | Mathematical Methods of Statistics[END_REF]

F(θ ) = E y|θ=θ -∆ θ θ log p(y|θ) |θ=θ . ( 4
)
The practical evaluation of the FIM for modulated signals is generally quite tedious because of the symbols a which are nuisance parameters. In order to circumvent this problem a "Modified" CRB has been proposed [START_REF] D'andrea | The modified Cramér-Rao bound and its application to synchronization problems[END_REF] and comes from the inversion of the following matrix:

G(θ ) = E a E y|a,θ=θ -∆ θ θ log p(y|a, θ) |θ=θ .
(5) Unlike the SCRB, it is usually possible to derive a closed form expression of the MCRB. In particular, in [START_REF] D'andrea | The modified Cramér-Rao bound and its application to synchronization problems[END_REF], the analytical expression associated to a constant phase offset is derived:

MCRB = 1 2K E s N 0 -1 , ( 6
)
where K is the number of observations. Note that although easier to evaluate than the SCRB, the MCRB is also generally looser. In the particular case of high-SNR transmissions, it has however been shown [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF] that the MCRB corresponds to the SCRB. We can also notice that equation ( 6) coincides with the CRB provided in [START_REF] Rife | Single tone parameter estimation from discrete-time observations[END_REF]. This bound is computed in the case of an unmodulated scenario or, equivalently, of a DA scenario.

B. The Bayesian Cramér-Rao Bound and its Modified Version

The standard and the modified CRBs are not suited to time-varying parameter estimation. In particular, they do not take into account the statistical dependence which may exist between phase offsets at different instants.

This dependence is naturally considered within the Bayesian framework: on the one hand, prior distribution p(θ) implicitly models the time dependence between stochastic phase offsets; on the other hand, Bayesian framework is intrinsically based on the knowledge of a vector prior distribution p(θ).

Within this context, a Bayesian CRB (BCRB) has been proposed in [START_REF] Trees | Detection, Estimation and Modulation Theory[END_REF] such that

E y,θ ( θ(y) -θ) ( θ(y) -θ) T ≥ BCRB. (7)
Unlike the SCRB and the MCRB, the BCRB does not depend on a particular value θ . The BCRB is the inverse of the Bayesian Information Matrix (BIM), which from [START_REF] Trees | Detection, Estimation and Modulation Theory[END_REF] can be written as

B = E θ F(θ) + E θ -∆ θ θ log p(θ) , (8)
where F(θ) is the FIM defined in [START_REF] Meyr | Digital Communication Receivers: Synchronization, Channel Estimation, and Signal Processing[END_REF]. Each term of the BIM (8) plays a different role: the first term can be interpreted as the average information about θ brought by the observations y whereas the second term can be regarded as the information available from the prior knowledge of θ, i.e., p(θ). This term allows us to take into account the time dependence between phase offsets at different instants.

Similarly to the MCRB, a modified version of the BCRB has been considered in [START_REF] Bobrovsky | Some classes of global Cramér-Rao bounds[END_REF]. This Modified BCRB (MBCRB) is the inverse of the following information matrix

C = E θ G(θ) + E θ -∆ θ θ log p(θ) , (9)
where G(θ) is the modified Fisher information matrix defined in [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF]. The MBCRB is usually easier to compute than the BCRB but it is also looser [START_REF] Bobrovsky | Some classes of global Cramér-Rao bounds[END_REF].

IV. THE OFF-LINE BCRB AND ITS MODIFIED VERSION FOR DYNAMICAL PHASE ESTIMATION

In this section, we derive an analytical expression of the BCRB associated to an off-line carrier-phase-offset estimation. We proceed as follows. We first obtain analytical expressions of the two terms contributing to the BIM [START_REF] Cowley | Phase and frequency estimation for PSK packets: bounds and algorithms[END_REF]. Then, the expression of the diagonal elements of the inverse of the BIM are derived. In addition, we also give an analytical expression of the MBCRB.

A. Computation of E θ [F(θ)]

E θ [F(θ)] corresponds to the first term in the right-hand side of equation ( 8). This term evaluation requires the computation of the Fisher information matrix F(θ), which in turn requires the evaluation of the Hessian of the log-likelihood function log p(y|θ). Using the observation model defined in section II, one has that the log-likelihood function can be expanded as

log p(y|θ) = log a p(y|a, θ) p(a). ( 10
)
Using the whiteness of the noise and the independence of the transmitted symbols, one then obtains that

∆ θ θ log p (y|θ) = K k=1 ∆ θ θ log p (y k |θ k ) . ( 11
)
It is important to note that each term of the sum ( 11) is a matrix with only one non-zero element at most, namely,

∆ θ θ log p (y k |θ k ) k,k = ∂ 2 ∂θ 2 k log p (y k |θ k ) . ( 12
)
As a direct consequence, the Hessian ∆ θ θ log p (y|θ) is a diagonal matrix with the k th diagonal element given by equation [START_REF] D'andrea | The modified Cramér-Rao bound and its application to synchronization problems[END_REF]. Moreover, because of the circularity of the observation noise, the expectation of equation ( 12) with respect to p (y k |θ k ) does not depend on θ k . One then obtains

E θ [F (θ)] = J D I K , ( 13
)
where I K is the K × K identity matrix and J D is defined as follows

J D E y,θ - ∂ 2 log p (y k |θ k ) ∂θ 2 k . ( 14
) B. Computation of E θ [∆ θ θ log(p(θ))]
Due to the Wiener structure of the phase model ( 2), ∆ θ θ log p (θ) can be expanded as

∆ θ θ log p (θ) = ∆ θ θ log p (θ 1 ) + K k=2 ∆ θ θ log p (θ k |θ k-1 ) . ( 15
)
Let us detail the expressions of each term contributing to equation [START_REF] Tichavsky | Posterior Cramér-Rao bounds for discrete-time nonlinear filtering[END_REF]. The first term is a matrix with only one nonzero element, namely, the entry (1, 1) which is equal to

∆ θ θ log p (θ 1 ) 1,1 = ∂ 2 log p (θ 1 ) ∂θ 2 1 . ( 16
)
The other terms in [START_REF] Tichavsky | Posterior Cramér-Rao bounds for discrete-time nonlinear filtering[END_REF] are matrices with only four non-zero elements, namely, the entries

(k-1, k-1), (k-1, k), (k, k-1)
and (k, k). Due to the Gaussian nature of the noise, one finds

∆ θ θ log p (θ k |θ k-1 ) k,k = ∆ θ θ log p (θ k |θ k-1 ) k-1,k-1 = -1 σ 2 w , ( 17
)
∆ θ θ log p (θ k |θ k-1 ) k,k-1 = ∆ θ θ log p (θ k |θ k-1 ) k-1,k = 1 σ 2 w . ( 18
)
Taking ( 16), ( 17) and ( 18) into account, one finally obtains

-E θ [∆ θ θ log p(θ)] =            1 σ 2 w -E θ1 ∂ 2 ∂θ 2 1 log p(θ 1 ) -1 σ 2 w 0 . . . 0 -1 σ 2 w 2 σ 2 w -1 σ 2 w . . . . . . . . . . . . . . . . . . 0 . . . -1 σ 2 w 2 σ 2 w -1 σ 2 w 0 . . . 0 -1 σ 2 w 1 σ 2 w            (19)
In the sequel, for the sake of conciseness, we will set

E θ 1 ∂ 2 log p(θ1) ∂θ 2 1
= 0. This corresponds to the case of a noninformative prior about θ 1 (see [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF]). The reasoning which is made in the following could however be extended to the general case, at the expense of some more difficult derivations and would not modify our results interpretation.

C. Analytical Expression of the Off-Line BCRB

In this subsection, an analytical expression of the diagonal elements of the inverse of the BIM is derived. From ( 13) and ( 19) one readily obtains that the BIM has a particular mathematical structure, i.e.,

B K = b        A + 1 1 1 A 1 . . . . . . . . . 1 A 1 1 A + 1        , ( 20
)
where A and b are defined as A -σ 2 w J D -2 and b -1/σ 2 w . In particular, B K is a symmetric sparse matrix having an (almost) periodical structure. Based on these observations, we show in Appendix I that the k th diagonal element of B -1 K can be expressed as

B -1 K k,k = 1 |B K | ρ 2 1 (b + r 1 ) 2 r K-3 1 + ρ 2 2 (b + r 2 ) 2 r K-3 2 - b 2 A -2 (r k-2 1 r K-k-1 2 + r K-k-1 1 r k-2 2 ) , (21)
where

r 1 1 σ 2 w + J D 2 1 -1 + 4 J D σ 2 w , ( 22
)
r 2 1 σ 2 w + J D 2 1 + 1 + 4 J D σ 2 w , ( 23
)
ρ 1 1 + 4 σ 2 w J D -1 -2 σ 2 w J D 2 1 + 4 σ 2 w J D , ( 24
)
ρ 2 1 + 4 σ 2 w J D + 1 + 2 σ 2 w J D 2 1 + 4 σ 2 w J D , ( 25
)
Note that (21) together with ( 22)-( 25) provides an analytical expression of the off-line BCRB associated to the estimation of θ k . This function of index k depends only on three parameters: the number of observations K, the phase noise variance σ 2 w and the observation noise variance σ 2 n (or equivalently J D , see equation ( 14)). The evaluation of J D is discussed in Section V.

D. The Off-line MBCRB

We now consider the MBCRB, (see equation ( 9)). The second term in the right-hand side of (9) has already been computed in section IV-B. The first term E θ G(θ) requires the evaluation of G(θ) which corresponds to the modified Fisher information matrix defined in [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF]. Using the observation model, one has that ∆ θ θ log p(y|a, θ) is a diagonal matrix where

[∆ θ θ log p(y|a, θ)] k,k = ∂ 2 ∂θ 2 k log p(y k |a k , θ k ). (26)
Due to the Gaussian distribution of the noise, one further finds that

∂ 2 ∂θ 2 k log p(y k |a k , θ k ) = 2 σ 2 n -a * k y k e -jθ k . ( 27
)
Consequently, one obtains that

[G(θ)] k,k = E a E y|a,θ - 2 σ 2 n -a * k y k e -jθ k (28) = E a [ 2 σ 2 n |a k | 2 ] (29) = 2 σ 2 n . ( 30
) One can observe that [G(θ)] -1 k,k = 1 2 Es N0 -1
, which corresponds to the MCRB (6) for K = 1 observation. As straightforwardly E θ G(θ) = 2/σ 2 n I K , the modified BIM is obtained exactly like the BIM of subsection IV-C with J M 2 σ 2 n playing the role of J D in equation [START_REF] Bobrovsky | Some classes of global Cramér-Rao bounds[END_REF]. Note, however, that the MBCRB is usually looser than the BCRB.

V. EVALUATION OF J D AND ASYMPTOTE-BASED LOWER

BOUND

In this section we first point out that, in the general case, the evaluation of J D implies to resort to numerical integration. We then derive the high-SNR and low-SNR asymptotes of the BCRB, and emphasize that their evaluation is easy. Finally, we show that these asymptotes are themselves lower bounds on the MSE. In particular, we illustrate that the combination of the low and high SNR asymptotes leads to an alternative (tight) lower bound, whose evaluation is straightforward.

A. Evaluation of J D

In this part, we calculate an expression of J D defined in [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF]. First, using the Gaussian nature of the noise and the equiprobability of the data symbols, one finds that

log p(y k |θ k ) = log 1 πσ 2 n e - 1+|y k | 2 σ 2 n cosh 2 σ 2 n y k e -jθ k . ( 31
)
Taking the second derivative of (31), one easily obtains that

∂ 2 log p(y k |θ k ) ∂θ 2 k = - 2 σ 2 n (x k ) tanh 2 σ 2 n (x k ) + 4 σ 4 n 2 (x k ) 1 -tanh 2 2 σ 2 n (x k ) , (32)
where x k y k e -jθ k . In the general case, the expectation of (32) with respect to p(y|θ) does not have any simple analytical solution. Hence, in practice, we have to resort to either numerical integration methods or some approximations.

In the following, we present both the high-SNR and the low-SNR approximations of the BCRB.

B. High-SNR BCRB asymptote

In this part, we investigate the BCRB behavior at high SNR. From the definition of the BIM [START_REF] Cowley | Phase and frequency estimation for PSK packets: bounds and algorithms[END_REF], only the first term, i.e., E θ [F(θ)], depends on the SNR. Moreover, we have from ( 13) that E θ [F(θ)] is fully characterized by J D [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF]. Hence, in the remainder of this section, we focus on the behavior of J D .

At high SNR (i.e., σ 2 n → 0), the tanh-function in (32) can be approximated as

tanh 2 σ 2 n (x k ) ≈ sign (x k ) . ( 33
)
Hence, replacing tanh in (32) by its approximation (33) and using the definition of J D [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF], one has (see Appendix III)

J D ≈ E y k ,θ k - 2 σ 2 n (x k )sign (x k ) (34) = 2 σ n √ π e -1 σ 2 n + 2 σ 2 n erf 1 σ n . ( 35
)
So, taking the limit of (35) for σ 2 n tending to zero, one has that

J Dh lim σ 2 n →0 J D = 2 σ 2 n . ( 36
)
We see therefore from (30) and (36), that the high-SNR asymptote of the BCRB (i.e., J -1 Dh ) is equal to the MBCRB. This corroborates the result derived by Moeneclaey [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF] in the non-Bayesian case for a scalar parameter.

C. Low-SNR BCRB asymptote

We now consider the low-SNR asymptote of the BCRB in the NDA BPSK context. Following the same reasoning as before, we equivalently focus on the behavior of J D for σ 2 n tending to infinity.

From (32), using the fact that tanh(z) ≈ z around z = 0, we obtain

∂ 2 log p (y k |θ k ) ∂θ 2 k ≈ - 2 σ 2 n (x k ) 2 + 2 σ 2 n (x k ) 2 . ( 37
)
Plugging (37) into ( 14), one readily obtains that the BCRB low-SNR asymptote is

J Dl lim σ 2 n →∞ J D = 4 σ 4 n . ( 38
)
In [START_REF] Steendam | Low-SNR limit of the Cramér-Rao bound for estimating the time delay of a PSK, QAM, or PAM waveform[END_REF], the authors obtain a result similar to (38) in a non-Bayesian scenario for a scalar parameter.

D. Asymptote-based Lower Bound

In this subsection, we show that the combination of the low and high-SNR BCRB asymptotes, presented in the previous subsections, still leads to a lower bound on the MSE. In particular, we define the "Asymptotic" BCRB (ABCRB) as ABCRB min(J Dl , J Dh )

I K + E θ [-∆ θ θ log p(θ)] -1 , ( 39
)
which can be evaluated with equations ( 19), (36) and (38). We now show that

MBCRB ≤ ABCRB ≤ BCRB. ( 40
)
The proof of the first inequality directly comes from the definition of the ABCRB and of the MBCRB. The second inequality is straightforward to show when J Dh ≤ J Dl because in such a case, we have

ABCRB = J Dh I K + E θ [-∆ θ θ log p(θ)] -1 , = MBCRB. (41)
Since we know that BCRB ≥ MBCRB (see [START_REF] Bobrovsky | Some classes of global Cramér-Rao bounds[END_REF]), ( 40) is proven for J Dh ≤ J Dl . The proof for J Dl < J Dh is more lengthy and is detailed in Appendix IV. The ABCRB is thus an easy-to-evaluate lower bound on the MSE. This Bayesian bound which does not require any Monte-Carlo simulation is compared to the BCRB in Section VII. This kind of reasoning only depends on the observation model: the ABCRB is relevant in the case of a BPSK transmission. However one can adapt this asymptote-based lower bound on any phase model system.

VI. THE ON-LINE BCRB : SEQUENTIAL CRAM ÉR-RAO BOUND

Up to this section, we have focused on the off-line scenario. We now show how the previous results can be used in the case of an on-line synchronization mode. In this mode, the receiver updates the observation vector y k 1 in order to estimate θ k : only the past and the current observations are available.

In order to evaluate the performance of this kind of estimator, a Posterior (on-line) Cramér-Rao Bound was derived in [START_REF] Tichavsky | Posterior Cramér-Rao bounds for discrete-time nonlinear filtering[END_REF]. Tichavsky et al. provide a method for updating the BIM from the time index k to the time index k + 1. This method was already applied in [START_REF] Amblard | Phase tracking: what do we gain from optimality? Particle filtering versus phase-locked loops[END_REF] on the same model described by equations ( 1) and ( 2): the lower bound for the last state parameter θ k is given by the following recursive sequence

C k+1 = σ 2 w + C k J D σ 2 w + 1 + J D C k and C 1 = 1 J D . ( 42
)
On the other hand, using our previous derivations, we can provide an alternative expression of the on-line BCRB. Indeed, the on-line BCRB associated to observation vector y K 1 is clearly equal to entry (K, K) of the inverse of the BIM, i.e., B -1 K K,K . Using expression ( 21), we therefore have

C K = B -1 K K,K = 1 |B K | ρ 2 1 (b + r 1 ) 2 r K-3 1 + ρ 2 2 (b + r 2 ) 2 r K-3 2 - b 2 A -2 (r K-2 1 r -1 2 + r -1 1 r K-2 2 ) . ( 43
)
The behavior of the off-line ( 21) and of the on-line (43) BCRBs will be studied and compared in the next section.

In addition, it is shown in Appendix II that this bound decreases with time index K to the following limit:

lim K→∞ B -1 K K,K = -σ 2 w + σ 4 w + 4 σ 2 w J D 2 . ( 44
)
Asymptotically with the observation size, the estimation performance is limited both by the phase noise and the observation noise (due to J D ). Therefore, unlike the non-Bayesian CRB [START_REF] Cowley | Phase and frequency estimation for PSK packets: bounds and algorithms[END_REF], the BCRB does not necessarily tend to zero when the number of observations goes to infinity.

VII. DISCUSSION

In this section we bring to the fore the behavior of the previous bounds, namely the off-line and the on-line BCRBs and ABCRB, by considering different scenarios.

We first consider a transmission disturbed by an AWGN with variance σ 2 n = 4 and phase noise with variance σ 2 w = 0.16 rad 2 . Figure 1 superimposes versus time index, the on-line BCRB (see equation ( 43)), its asymptote described by equation (44), and the off-line BCRBs for different block-observation lengths K (see equation ( 21)). We then obtain the lower bound for each phase offset θ k according to the considered scenario.

In the off-line context, one can see that the best phase estimate is achieved at mid-block, whereas the estimates are likely to be poorer at the block boarder (the proof of this property is detailed in Appendix I-D). This stems from the fact that in the center position of the phase vector θ one has more adjacent (past or future) and strongly correlated variables than at the boarder of the phase vector. Stated more precisely, from equation ( 2), the correlation coefficient between θ k and θ k+l (where l is a positive integer) is equal to

E θ [θ 2 1 ]+(k-1) σ 2 w E θ [θ 2 1 ]+(k+l-1) σ 2 w
, and thus even if a boarder phase is estimated with the same number of observations than a mid-size position phase, the correlation coefficient at the boarder are globally poorer than at the center and the information cannot be exploited as well.

Concerning the on-line bound, at the beginning when the number of observations increases, the estimator takes into account more and more information and the estimation is improved; the bound thus decreases and converges to its asymptote: the estimation performance is then limited by the phase noise and the observation noise independently of the number of observations taken into account. We now analyze the bound behavior versus the SNR over a block of K = 20 BPSK transmitted symbols. We illustrate some results for two distinct phase-noise variances: σ 2 w = 0.01 rad 2 and σ 2 w = 0.16 rad 2 . Figure 2 superimposes for the two previous phase-noise variances, the BCRBs (see equation ( 21)) evaluated over 10 7 Monte-Carlo trials and the ABCRBs (see equation ( 39)) with k = K+1 2 versus the SNR. One can recover on Figure 2 for σ 2 w = 0.16 rad 2 and SNR = -6 dB, the minimum of the off-line BCRB displayed on Figure 1 for K = 20. We distinguish three SNR ranges.

• At high SNR, we notice that the ABCRBs in dashed lines and the corresponding BCRBs logically merge. Moreover, the ABCRB is close to the BCRB in the whole range of SNR for the two phase noise variances considered. At high SNR, both BCRBs converge to the deterministic MCRB evaluated for K = 1 observation [START_REF] D'andrea | The modified Cramér-Rao bound and its application to synchronization problems[END_REF], [START_REF] Gini | The modified Cramér-Rao bound in vector parameter estimation[END_REF] i.e., M CRB = 1/J Dh . In this range of SNR, the information provided by each observation y k is preponderant over the a priori knowledge of θ. Consequently in this case the observations are reliable enough to only take into account the present observation y k in order to estimate θ k and this is why the bounds converge to the well-known MCRB for any variance σ 2 w . As the a priori distribution of θ has no influence, the Bayesian problem tends to a deterministic phase estimation problem where we estimate K independent phases θ k with K independent observations. • In mid-range SNR, when the phase variance σ 2 w decreases, the BCRB gets closer to the deterministic MCRB evaluated for K = 20 observations. This is because in the limiting case for which the phase variance tends to zero, the time-varying phase estimation problem can be simplified to a deterministic constant-phase estimation problem. The MCRB is then a lower bound of the BCRBs for the same number of observations K = 20.

• At low SNR, the BCRB and the ABCRB coincide. In this SNR-range, the observation noise masks the a priori dynamic phase evolution. The Bayesian bound tends to the deterministic asymptote 1

K J Dl = 4 K Es N 0 -1
given by [START_REF] Steendam | Low-SNR limit of the Cramér-Rao bound for estimating the time delay of a PSK, QAM, or PAM waveform[END_REF]. One can note that the phase evolution is not stationary i.e., var (θ k ) = var (θ k-1 ) + σ 2 w , and consequently that uncertainty over the phase offset grows with time. Moreover as the phase model ( 2) is not restricted to a finite horizon, the MSE is not upper bounded. This is why there is no observed saturation as one could expect with a traditional Bayesian bound.

VIII. CONCLUSION

In this contribution, we have derived an analytical expression of a Bayesian Cramér-Rao Bound for the estimation of a realistic dynamical phase. We provide a closed-form expression in the case of a BPSK transmission disturbed by AWGN in a NDA context. The reasoning can be extended to any M-ary Phase Shift Keying transmissions. Furthermore, we can readily predict the bound behavior from its low and high-SNR asymptotes. We combine these asymptotes and we give a new Bayesian bound which is easy to evaluate at the price of being slightly lower than the BCRB. These bounds are useful when analyzing the performance of actual phase-tracking estimators in stressing environments for both the offline and on-line scenarios.

APPENDIX I EXPRESSION OF DIAGONAL ELEMENTS OF THE BCRB

In this appendix, we detail the calculus to obtain the diagonal elements [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF] of the inverse of the BIM [START_REF] Bobrovsky | Some classes of global Cramér-Rao bounds[END_REF]. We first use the classical matrix-inversion formula

[B -1 K ] k,k = C k,k |B K | , ( 45

A. Preliminary calculus : Determinant d k

In this subsection, we derive an analytical expression of the determinant d k of a k × k matrix D k defined as

D k = b        A 1 1 A 1 . . . . . . . . . 1 A 1 1 A        (46)
Expanding d k along the first and the second column, one obtains the recursive equation

d k = Abd k-1 -b 2 d k-2 with d 0 = 1 and d 1 = bA. ( 47
)
{d k } ∞ k=0 is thus a linear recurrent sequence and its characteristic polynomial p (x) = x 2 -Abx + b 2 has real roots r 1 and r 2 . Then, using the initial terms, an analytical expression of d k is given by

d k = ρ 1 (r 1 ) k + ρ 2 (r 2 ) k for k ∈ N, (48)
where

r 1 = b 2 A + A 2 -4 , ρ 1 = √ A 2 -4 + A 2 √ A 2 -4 , r 2 = b 2 A -A 2 -4 , ρ 2 = √ A 2 -4 -A 2 √ A 2 -4 . B. Determinant |B K |
In this subsection, we derive an expression of |B K | by using the preliminary result derived Appendix I-A. In a first step, expanding this determinant along the first column, we obtain the sum of two cofactors only. Then, expanding in turn these cofactors along the last column and using (47), we have

|B K | = (A + 2) b d K-1 , ( 49
)
or with (48),

|B K | = (A + 2) b ρ 1 (r 1 ) K-1 + ρ 2 (r 2 ) K-1 . ( 50
) C. Cofactor C k,k
In order to calculate the cofactor C k,k , one has to delete the k th row and the k th column of the matrix B K and one obtains a two-block diagonal matrix. The upper (respectively lower) block is noted U B K (respectively L B K ). The cofactor which is obtained from the determinant of the previous matrix is thus the product of two determinants:

C k,k = det(U B K )det(L B K ) (51) = (bd k-2 + d k-1 ) (bd K-k-1 + d K-k , ) , (52)
where d k is calculated in Appendix I-A.

D. Expression of the diagonal elements of B -1

K

An analytical expression of the diagonal elements of the inverse of ( 20) is now derived. Then, the behavior of these elements along their index k and the matrix size K × K is studied.

Rewriting (45) and using (52), one has

B -1 K k,k = (bd k-2 + d k-1 ) (bd K-k-1 + d K-k ) |B K | . ( 53
)
Then, using (48), B -1 K k,k can be expanded as follows

B -1 K k,k = 1 |B K | ρ 2 1 (b + r 1 ) 2 r K-3 1 + ρ 2 2 (b + r 2 ) 2 r K-3 2 - b 2 A -2 (r k-2 1 r K-k-1 2 + r K-k-1 1 r k-2 2 ) ,
which is equivalent to [START_REF] Kay | Fundamentals of statistical signal processing: estimation theory[END_REF]. For a fixed block length K from ( 21), note that B -1

K k,k
depends on the index k through the expression g K (k) where

g K (x) = r x-2 1 r K-x-1 2 + r K-x-1 1 r x-2 2 . ( 54
)
After analyzing this function, by classically studying the sign of dg K dx , one readily finds that • the symmetry of g K (x) with respect to K+1 2 implies that

B -1 K k,k = B -1 K l,l with l = K + 1 -k, • the minimal diagonal element is the mid-coefficient B -1 K K+1 2 , K+1 2 ,
• the maximal diagonal elements are B -1 K 1,1 and B -1 K K,K .

APPENDIX II DERIVATION OF EXPRESSION (44)

Using (53), B -1 K K,K can be written as follows:

B -1 K K,K = A + 1 b (A + 2) - 1 b 2 (A + 2) d K d K-1 , ( 55
)
where d K is defined in Appendix I-A. Since B -1 K K,K is only a function of d K d K-1 , we consider the fixed points of a sequence {u n } ∞ n=0 defined as u n d n bd n-1 = Au n-1 -1 u n-1 , (56)
with u 1 = A. Clearly, this sequence is strictly increasing and converges to u ∞ = 1 2 A -√ A 2 -4 . Combining this result with (55), we have that

B -1 K K,K = A + 1 b (A + 2) - 1 b (A + 2) u K , ( 57
)
is a strictly decreasing sequence with the following limit:

lim K→∞ B -1 K K,K = A + 1 b (A + 2) - A - √ A 2 -4 2b (A + 2) . ( 58
)
Using the definitions of A and b in subsection IV-C, one obtains

lim K→∞ B -1 K K,K = -σ 2 w + σ 4 w + 4 σ 2 w J D 2 . ( 59
) APPENDIX III HIGH-SNR ASYMPTOTE OF J D Defining F H (θ k ) E y k |θ k 2 σ 2 n (x k )sign (x k ) , (34) becomes J D ≈ E θ k [F H (θ)] . ( 60
)
Since the noise affecting the observation is circular, x k = a k + ñk where ñk = n k e -jθ k is a circular Gaussian noise with variance σ 2 n , therefore, we find

F H (θ k ) = E a k ,ñ I k |θ k 2 σ 2 n a k + ñI k sign a k + ñI k , where ñI k {ñ k }.
Using first the definition of the expectation and then the definition of the error function erf(•)

F H (θ k ) = a k ∈{-1,1} 1 σ 2 n +∞ -∞ a k + ñI k × sign a k + ñI k p(ñ I k )dñ I k = 1 σ 2 n 2 +1 -1 p(ñ I k )dñ I k + 4 +∞ 1 ñI k p(ñ I k )dñ I k , = 2 σ n √ π e -1 σ 2 n + 2 σ 2 n erf 1 σ n . ( 61
)
Since (61) does not depend on θ k , we finally obtain

J D ≈ E θ k [F H (θ k )] = 2 σ n √ π e -1 σ 2 n + 2 σ 2 n erf 1 σ n . ( 62
)
APPENDIX IV PROOF THAT BCRB > ABCRB FOR J Dl < J Dh To prove (40), we show that

J D -J Dl < 0, (63)
Using the definition of J D and J Dl , i.e., [START_REF] Moeneclaey | On the true and the modified Cramér-Rao bounds for the estimation of a scalar parameter in the presence of nuisance parameters[END_REF] and (38), we have

J D -J Dl =E y k ,θ k 4 σ 4 n 2 (x k ) tanh 2 2 σ 2 n (x k ) (64)
+ 2 σ 2 n (x k ) tanh 2 σ 2 n (x k ) - 2 σ 2 n (x k ) .
For the sake of conciseness, we use the following notations: X (x k ) and Y (x k ). Taking then into account that the expectation in (64) is independent of θ k , we have

J D -J Dl = E X,Y 4 σ 4 n Y 2 tanh 2 2 σ 2 n X (65) + 2 σ 2 n X tanh 2 σ 2 n X - 2 σ 2 n X .
Since the noise affecting the observations is circular, X and Y are independent. Moreover Y is a zero-mean random Gaussian variable with variance σ 2 n /2. Hence,

J D -J Dl = E X [f (X)] , ( 66
)
where f (X) Then we rewrite (66) as

2 σ 2 n X tanh 2 σ 2 n X -2 σ 2 n X 2 + tanh 2 2
J D -J Dl = 1 2 E X|a=+1 [f (X)] + E X|a=-1 [f (X)] .
(67)

Noting that f (X) = f (-X), we deduce from a change of variable (X = -X) that the two terms in the right-hand side of (67) are equal. This equation thus becomes

J D -J Dl = E X|a=+1 [f (X)] . ( 68
)
We now prove that E X|a=+1 [f (X)] is negative. We use Stein's lemma [START_REF] Stein | Estimation of the mean of a multivariate normal distribution[END_REF]: suppose Z is a normally distributed random variable with expectation µ and variance v. (69)

Applying this lemma to the random variable X with expectation µ = 1 and variance σ 2 n /2 and to g(X) = tanh 2 σ 2 n X , we obtain

σ 2 n 2 E X|a=+1 2 σ 2 n 1 -tanh 2 2 σ 2 n X (70) = E X|a=+1 (X -1) tanh 2 σ 2 n X .
Moreover, we directly have

E X|a=+1 4 σ 4 n X 2 = 2 σ 2 n + 4 σ 4 n . ( 71
)
Then, plugging equations (70) and (71) into equation (68)

J D -J Dl = E X|a=+1 2 σ 2 n tanh 2 σ 2 n X - 4 σ 4 n . (72)
Considering that X is a Gaussian variable with expectation µ = 1, we now show that E X|a=+1 tanh 2

σ 2 n X -2 σ 2 n X < 0.
Suppose Z is a normally distributed random variable with expectation µ = 1 and variance v. Denoting its pdf by p(Z), we have for any λ > 0 As p(Z) -p(-Z) > 0 for any Z > 0, and as tanh(λZ) < λZ for any Z > 0 and λ > 0, we actually have E Z [tanh (λZ) -λZ] < 0, so that

E Z [tanh (λZ)] < λ, λ > 0. ( 75
)
As X is a normally distributed random variable with expectation µ = 1, applying this result with λ = 2 σ 2 n , we obtain

E X|a=+1 2 σ 2 n tanh 2 σ 2 n X - 4 σ 4 n < 0, (76)
and finally from equations ( 72) and (76) we have proved equation ( 63)

J D -J Dl < 0, σ n > 0. (77)

Fig. 1 .

 1 Fig. 1. BCRBs versus the number of observations for σ 2 n = 4, σ 2 w = 0.16 rad 2 (J D is evaluated over 10 5 Monte-Carlo trials).

Fig. 2 .

 2 Fig. 2. Comparison between Bayesian CRBs and Asymptote-based CRBs for two different phase noise variances σ 2 w = 0.01 rad 2 and σ 2 w = 0.16 rad 2 .

  ) where C k,k is the cofactor of the element [B K ] k,k and |B K | is the determinant of matrix B K . In the sequel, we focus on the derivation of particular expressions of C k,k and |B K |. We first present some preliminary calculus (sections I-A to I-C) aiming at simplifying the expansion of |B K | and C k,k . Then, in section I-D, we derive a closed-form expression of the diagonal elements of B -1 K .

σ 2 nX

 2 and where X is a random variable with the following probability density function (pdf)

  Furthermore suppose g is a function for which the two expectations E[(Z -µ)g(Z)] and E[g (Z)] both exist. Then vE Z [g (Z)] = E Z [(Z -µ)g(Z)] .

E

  Z [tanh (λZ) -λZ] λZ) -λZ] (p(Z) -p(-Z)) dZ. (74)