Data Smashing 2.0: Sequence Likelihood (SL) Divergence For Fast Time Series Comparison
Yi Huang Ishanu Chattopadhyay
  
Range Adaptation for 3D Object Detection in LiDAR
    [
    December 30, 2023
=================================================




Recognizing subtle historical patterns is central to modeling and forecasting problems in time series analysis. Here we introduce and develop a new approach to quantify deviations in the underlying hidden generators of observed data streams, resulting in a new efficiently computable universal metric for time series. The proposed metric is universal in the sense that we can compare and contrast data streams regardless of where and how they are generated, and without any feature engineering step. The approach proposed in this paper is conceptually distinct from our previous work on data smashing <cit.>, and vastly improves discrimination performance and computing speed. The core idea here is the generalization of the notion of KL divergence often used to compare probability distributions to a notion of divergence in time series. We call this generalization the sequence likelihood (SL) divergence and show that it can be used to measure deviations within a well-defined class of discrete-valued stochastic processes. We devise efficient estimators of SL divergence from finite sample paths, and subsequently formulate a universal metric useful for computing distance between  time series produced by hidden stochastic generators. We illustrate the superior performance of the new  metric with synthetic data against the original data smashing algorithm and dynamic time warping (DTW) <cit.>. Pattern disambiguation in two distinct applications involving electroencephalogram data and gait recognition is also illustrated. We are hopeful that the  metric introduced here will become an important component of the standard toolbox used in classification, clustering and inference problems in time series analysis.



 Universal Metric; Data Smashing; Beyond Dynamic Time Warping; Probabilistic Finite Automata; Time Series Clustering 




§ INTRODUCTION

Efficiently learning stochastic processes is a key challenge in analyzing  time-dependency in domains where randomness cannot be ignored. For such learning to occur, we need to define either a distance metric or, more generally, a measurement of similarity to compare and contrast time series. Examples of such similarity measurement from the literature include the classical l_p distances and l_p distances with dimensionality reduction <cit.>, the short time series distance (STS)<cit.>, which takes into account of irregularity in sampling rates, the edit based distances<cit.> with generalizations to continuous sequences<cit.>, and the dynamic time warping (DTW)<cit.>, which is used extensively in the speech recognition community. However these measurement of similarity all have either one or both of the following limitations. First, dimensionality reduction and feature selection heavily relies on domain knowledge and inevitably incurs trade-off between precision and computability. Most importantly, it necessitates the attention of human experts and data scientists. Secondly, when dealing with data from non-trivial stochastic process dynamics, state of the art techniques might fail to correctly estimate the similarity or lack thereof between exemplars. For example, suppose two sequences recording n tosses of a fair coins, use 1 to represent a head and 0, tail. The two sequences are extremely unlikely to share any similarity on the face value, i.e. they have a large pointwise distance, but they are generated by the same process. A good measurement of similarity should strive to disambiguate the underlying processes. The  metric introduced here addresses both these limitations.

When presented with finite sample paths, the  algorithm is specifically designed to estimate a distance between the generating models of the time series samples.  The intuition for the  metric  follows from a basic result in information theory: If we know the true distribution p of the random variable, we could construct a code with average description length H(p), where H(·) is the entropy of a distribution. If, instead, we used the code for a distribution q, we would need H(p) + D(p||q) bits on the average to describe the random variable. Thus, deviation in the distributions show up as KL divergence. If we can generalize the notion of KL divergence to processes, then it might be possible to quantify deviations in process dynamics via an increase in the entropy rate by the corresponding divergence. 

Our ultimate goal is to design an algorithm that operates on a pair of data streams taking values in a finite alphabet. Nevertheless, to establish the correctness of our algorithm, we need to decide on a specific scheme for representing stochastic processes taking values in the alphabet. We further assume that our processes are ergodic and stationary. The specific modeling paradigm for modeling stochastic processes we use in this paper is called Probabilistic Finite-State Automaton, or PFSA for short, which has been studied in <cit.>. PFSA can model discrete-valued stochastic processes that are not Markov of any finite order<cit.>. It is also shown in <cit.> to be able to approximate any hidden Markov model (HMM) with arbitrary accuracy. Moreover, PFSA has the property that many key statistical quantities of the processes they generate, such as entropy rate<cit.> and KL-divergence<cit.>, have closed-form formulae. Here we want to point out the resemblance of the PFSA model to the variational autoencoder (VAE) <cit.> framework. The inference of PFSA from the input can be thought as the training of the encoder in a VAE, and the performance of both the VAE and the PFSA model are evaluated by the log-likelihood of input as being generated by the inferred models.

The work that has inspired the development of  is the data smashing algorithm () proposed in <cit.>.  is also based on PFSA modeling and designed directly to represent the similarity between the generating models rather than sample paths. However, as while as both  and  have the advantage of not requiring dimensionality reduction or domain knowledge for feature extraction,  is much more computationally efficient than . 

The remaining of the paper is organized as follows. In Sec. <ref>, we introduce basic concepts of stochastic processes and establish the correspondence between processes and labeled directed graphs via the core concept of causal state. The definition and basic properties of PFSA are introduced by the end of Sec. <ref>. In Sec. <ref>, we answer the question of when a stochastic process has a PFSA generator. An inference algorithm, , of PFSA is given in Sec. <ref>. In Sec. <ref>, we introduce the notion of irreducibility of PFSA and the closed-form formulae for entropy rate and KL divergence of the processes generated by irreducible PFSA. We conclude the section with log-likelihood convergence. In Sec. <ref> we introduce the definition of  together with quantization of continuous sequences. The comparison of  to  and  is given in Sec. <ref>. In Sec. <ref>, we apply  to two real world problems.



§ FOUNDATION




 §.§ Stochastic Processes and Causal States


In this paper we study the generative model for stationary ergodic stochastic processes <cit.> over a finite alphabet. Specifically, we consider a set of Σ-valued random variables X_t_t∈ℕ^+ indexed by positive integers representing time steps. By stationary, we mean strictly stationary, the finite-dimensional distributions <cit.> are invariant of time. By ergodic, we mean that all finite-dimensional distributions can be approximated with arbitrary accuracy with long enough realization. We are especially interested in processes in which the X_is are not independent. 

We denote the alphabet by Σ and use lower case Greeks (e.g. σ or τ) for symbols in Σ. We use lower case Latins (e.g. x or y) to denote sequences of symbols, x=σ_1σ_2…σ_n for example, with the empty sequence denoted by λ. The length of a sequence x is denoted by |x|. The set of sequences of length d is denoted by Σ^d, and the collection of sequences of finite length is denoted by Σ^⋆, Σ^⋆= ⋃_d=0^∞Σ^d. We use Σ^ω to denote the set of infinitely long sequences, and xΣ^ω to denote the collection of infinite sequences with x∈Σ^⋆ as prefix. We note that, since all sequences can be viewed as prefixed by λ, we have λΣ^ω=Σ^ω.

We note that = xΣ^ω| x ∈Σ^⋆. is a semiring over Σ^ω. Let PrX_1… X_n=σ_1…σ_n denote the probability of the process producing a realization with X_i = σ_i for i = 1, …, n, it is straightforward to verify that μ:→ [0, 1] defined by  

    μσ_1…σ_nΣ^ω = PrX_1… X_n=σ_1…σ_n,

is a premeasure on . By Charathéodory extension theorem, the σ-finite premeasure μ can be extended uniquely to a measure over =σ(), where σ() is the σ-field generated by . Denoting the measure also by μ, we see that every stochastic process induces a probability space Σ^ω, , μ over Σ^ω. In light of Eq. (<ref>) and also for notational brevity, we denote μxΣ^ω by Pr(x) when no confusion arises. We note that Pr(λ) = μλΣ^ω = μΣ^ω = 1. We refer to Chap. 1 of <cit.> as a more formal introduction to the measure-theory knowledge used here.

Taking one step further, and denoting the collection of all measures over Σ^ω, ℱ by ℳ_Σ, we see that we can get a family of measures in ℳ_Σ from a process in addition to μ.


For an observed sequence x∈Σ^⋆ with Pr(x)> 0, the measure μ_x is the extension to  of the premeasure defined on the semiring  given by

    μ_xyΣ^ω = Prxy/Pr(x),

for all y∈Σ^⋆.


Now we introduce the concept of Probabilistic Nerode Equivalence, which was first introduced in <cit.>. 

For any pair of sequences x, y∈Σ^⋆, x is equivalent to y, written as x∼ y, if and only if either Pr(x) = Pr(y) = 0, or μ_x = μ_y. 

One can verify that the relation defined above is indeed an equivalence relation which is also right-invariant in the sense that x∼ y ⇒ xz∼ yz, for all z∈Σ^⋆. We denote the equivalence class of sequence x by [x]. We note that μ_[x] is well-defined because μ_x=μ_y for x, y∈[x]. An equivalence class is also called a causal state <cit.> since the distribution of future events preceded by possibly distinct x,y∈[x] are both determined by μ_[x]. We denote μ_[x]yΣ^ω by Pr_[x](y) when no confusion arises. We note that μ_[λ] = μ.


Since the equivalence class x∈Σ^⋆|Pr(x) = 0 plays no role in our future discussion, we ignore it as a causal state from this point on.




For any d∈ℕ^+, the d-th order derivative of an equivalence class [x], written as ϕ^d_[x], is defined to be the marginal distribution of μ_[x] on Σ^d, with the entry indexed by y denoted by ϕ^d_[x](y). The first-order derivative is also called the symbolic derivative in <cit.> since Σ^1 = Σ, and is denoted by ϕ_[x] for short. The derivative of a sequence is that of its equivalence class, ϕ^d_x = ϕ^d_[x]. We note that ϕ^d_λ is the marginal distribution of μ on Σ^d, and is denoted by ϕ^d for short.




 §.§ From Causal States to Probabilistic Automaton


From now on, we denote the set of causal states of a process by Q when no confusion arises. We start this section by showing that there is a labeled directed graph <cit.> associated with any stochastic process. 

For any q∈ Q and σ∈Σ such that Pr_q(σ) > 0, by right-invariance of probabilistic Nerode equivalence, there exists a q'∈ Q, such that xσ∈ q' for all x∈ q. Whenever the scenario described happens, we can put a directed edge from q to q' and label it by σ and Pr_q(σ), and by doing this for all q∈ Q and σ∈Σ, we get a (possibly infinite) labeled directed graph with vertex set Q. 

[An Order-One Markov Process]

We now carry out the construction described above on an order-1 Markov process <cit.> over alphabet Σ = 0, 1, in which X_t+1 follows a Bernoulli distribution conditioned on the value of X_t. Specifically we have

    PrX_t+1 = 0 | X_t = 0 = .6,    PrX_t+1 = 0 | X_t = 1 = .4.

Together with the specification PrX_1 = 0 = .5, we can check that the process is stationary and ergodic. The reason that we choose this process as our first example is because it has a small set of causal states of size 3. We list the causal states of sequences up to length 3 in Tab. <ref>.


Since μ_x is defined on an infinite dimensional space, we only show the symbolic derivative ϕ_x in Tab. <ref>, but we can verify that μ_x = μ_y if and only if ϕ_x=ϕ_y for this process.



Now, we conceptualize the labeled directed graph obtained from analyzing the causal states by an automaton structure <cit.>, which we call probabilistic finite-state automaton <cit.>, and show how we can get a stochastic process from it.

A probabilistic finite-state automaton G is specified by a quadruple Σ, Q, δ,, where Σ is a finite alphabet, Q is a finite set of states, δ is a partial map from Q×Σ to Q called transition map, and , called observation probability, is a map from Q to 𝐏_Σ, where 𝐏_Σ is the space of probability distributions over Σ. The entry indexed by σ of (q) is written as (q, σ). 




We call the directed graph (not necessarily simple with possible loops and multiedges) with vertices in Q and edges specified by δ the graph of the PFSA and, unless stated otherwise, we assume it to be strongly connected <cit.>, which means for any pair q, q'∈ Q, there is a sequence σ_1σ_2⋯σ_k, such that δq_i-1, σ_i = q_i for i = 1, 2, …, k with q_0=q and q_k = q'. 

To generate a sequence of symbols, assuming G's current state is q, it then outputs symbol σ with probability q, σ, and moves to state δ(q, σ). We see that δ is partial because δ(q, σ) is undefined when (q, σ) = 0.



Given a PFSA Σ, Q, δ,, the observation matrix  is the Q×Σ matrix with the (q,σ)-entry given by (q, σ), and the transition matrix Π is the Q×Q matrix with the (q, q')-entry, written as π(q, q'), given by 

    π(q, q') = ∑_σ: δ(q, σ)=q'(q, σ).


It is straightforward to verify that both Π and  are stochastic, nonnegative with rows of sum 1.


We borrow the terms observation matrix and transition matrix from the study of HMM <cit.>. However, we need to point out here that our model differs from the HMM in that, in HMM, the transition from the current state to the next one is independent of the symbol generated by the current state, while in PFSA, the current state and symbol generated together determine the next state the PFSA will be in. 


Unless specified otherwise, we assume the initial distribution to be the stationary distribution <cit.> of Π. We denote the stationary distribution of G by 𝐩_G, or by 𝐩 if G is understood.



Stochastic process generated by a PFSA G with distribution on states initialized with 𝐩_G is stationary and ergodic.



proof omitted.



Example <ref> shows that we may derive a PFSA from a stationary ergodic process, and Thm. <ref> shows that the process generated by the PFSA thus obtained is also stationary and ergodic. This motivates us to seek a characterization for stochastic processes that gives rise to a PFSA. Since the process in Example <ref> is an order-1 Markov process, which is the simplest non-i.i.d. process, it is legitimate to ask whether a process has to be Markov to have a PFSA generator. This desired characterization is obtained from studying the properties of causal states, which we do in the next section.



Table <ref> compares the three generative models of stochastic processes mentioned in this paper: Markov chain(MC), PFSA, and hidden Markov model(HMM). We note that a Markov chain produces a sequence of states, while sequences produced by PFSA and hidden Markov model take values in their respect output alphabets. We can also see that HMM can be considered as an extension to MC by adding an output alphabet and observation probabilities while PFSA are not directly comparable to either MC or HMM.






§ STOCHASTIC PROCESSES WITH PFSA GENERATOR




 §.§ Persistent Causal States



Let Q be the set of causal states of a stationary ergodic process. For every q∈ Q and d∈ℕ, let p_d(q) = ϕ^d[x] = q, the probability of length-d sequences who are equivalent to q. A causal state q is persistent if lim inf_d→∞p_d(q) >0, and transient if otherwise. We denote the set of persistent causal states by Q^+.



Here we borrow the term transient state from Markov chains literature, for example <cit.>, but we should note that the two concepts are not identical. A Markov chain never revisits transient states as soon it hits a recurrent state. However, although a transient causal state could never be revisited, as the q_λ in Example <ref>, it could also be revisited for infinitely many times. The transient states in MC and PFSA are similar in that the probability of a Markov chain being in a transient state diminishes as time increases, and a transient causal states also has lim inf_d→∞p_d(q) = 0. Since transient states can recur, we name the counterpart to transient causal state in PFSA by persistent causal state, not by recurrent state as in MC.


For any pair q, q'∈ Q, let π_q, q' = ∑_σ: qσ = q'ϕ_q(σ), where the expression qσ = q' is a shorthand for [xσ] = q' for all [x] = q. The following proposition shows that π_q, q' captures the flow of probability over causal states as sequence length increases.  

We have p_d(q') = ∑_q∈ Qπ_q, q'p_d-1(q) for each q'∈ Q and d ∈ℕ^+. Furthermore, there is no flow from a persistent state to a transient one, π_q, q' = 0 for q∈ Q^+ and q'∈ Q∖ Q^+.



proof omitted.





Let Q^+ be the set of persistent causal states of a stationary ergodic process 𝒫. Then, p(q) = lim_d→∞p_d(q) exists for every q∈ Q^+. Furthermore,  if Q^+ is finite and ∑_q∈ Q^+ p(q) = 1, the process generated by the PFSA G = Σ, Q^+, δ, with δ(q, σ) = qσ and (q, σ) = ϕ_q(σ) is exactly 𝒫. In fact, we have .𝐩_G|_q = p(q) for q∈ Q^+.



proof omitted.



[An Order-Two Markov Process]

Now, let us consider an order-2 Markov process over alphabet Σ = 0, 1, in which X_t+2 follows a Bernoulli distribution conditioned on the value of X_tX_t+1. More specifically, denoting PrX_t+2 = 0 | X_tX_t+1 = ij by p_ij for i, j∈0,1, we have p_00 = .3, p_01 = .2, p_10=.8, p_11 = .7. Together with the specification X_2 = 0|X_1 = 0 = 8/15, PrX_2 = 0|X_1 = 1 = 7/15, and PrX_1 = 0 = .5, we can check that the process is stationary and ergodic. We list the causal states of sequences up to length 3 in Tab. <ref>.


Since μ_x is defined on an infinite dimensional space, we only show ϕ_x in Tab. <ref>, but we can check that ϕ_x = ϕ_y if and only if μ_x = μ_y for this process. Since q_λ, q_0, q_1 only show up once, while q_00, q_01, q_10, q_11 appear repeatedly, we have Q^+ = q_00, q_01, q_10, q_11. With more detailed calculation, we can show that pq_00 = 4/15, pq_01 = 7/30, pq_10 = 7/30, and pq_11 = 4/15, which sum up to 1. According to Thm. <ref>, we can construct a PFSA with state set Q^+ that generates exactly the same process. We demonstrate the labeled directed graph constructed on Q in Fig. <ref>, and the PFSA is exactly the induced subgraph <cit.> on Q^+, which is also the unique strongly connected component of the graph. We can show that the stationary distribution of the PFSA is exactly (4/15, 7/30, 7/30, 4/15). 






[A PFSA on Three States]

In this example, we analyze the stochastic process generated by the PFSA on the right of Fig. <ref>. We nickname the PFSA by T. We show that Q of this process is infinite, while Q^+ is of size 3. We first notice that, no matter what state the PFSA resides, the sequence 11, and hence any sequence ending in 11, will take it to state r, which generates symbol 0 with probability .8, and 1 with probability .2. We also note that, whenever there are two consecutive 1s in a given sequence in Σ^⋆, we know for sure the state the PFSA resides. For example, sequence 110 will take the PFSA to s, and 1101, to q. On the left of Fig. <ref>, we show the probabilities of causal states [11], [110], [1101], and the sum of probabilities of all other causal states for sequence length d = 0, …, 25. We see from the bar plots that the sum of concentrations of [11], [110], and [1101] approaches 1 as d increases. We also point out that, with all numbers rounded up to three decimal places, p_25[11]=0.182, p_25[110] = 0.474, and p_25[1101] = 0.279, while the stationary distributions of the states r, s, and q are 0.190, 0.506, 0.304, respectively.

However, we also note that Q of the process is actually infinite by observing the fact that ϕ_0^d, where σ^d means σ repeated d times, are all distinctive.

We note that the process generated by this PFSA is not Markov, as implied by the infinity of Q. However, the fact that there are only three persistent causal states whose sum of probabilities approaches 1  allows it to have  a PFSA generator.    


[A Stochastic Process with Empty Q^+]

In this example, we analyze the stochastic process generated by the PFSA on the left of Fig. <ref>. We nickname the PFSA by S. We show that Q of this process is infinite while Q^+ is empty. Without run into details of the computation, we point out the fact that causal states of this process are also uniquely characterized by their symbolic derivatives, and the set ϕ(q)|q∈ Q is in one-to-one correspondence with ℤ. More specifically, we have

    Q = q_n|ϕ_q_n∝.5^n + 1 + 1, .5^n + .5, n∈ℤ.,

where ∝ means being proportional to, and 

    π_q_n, q_n+1 = 2/3.5^n+1 + 1/.5^n + 1,  π_q_n, q_-n+1 = 1/3.5^n-1 + 1/.5^n + 1,

with π_q_n, q_n+1 + π_q_n, q_-n+1 = 1, for all n∈ℤ. We demonstrate on the left of Fig. <ref> the contour of p_dq_n against n for sequence length d = 10, 20, …, 150. It takes some more work to show rigorously, but we can speculate that, for any fixed n∈ℤ, p_dq_n approaches 0 as d approaches infinity, as the curves flatten out with increasing d.






 §.§ Accumulation Causal States


We see from Example <ref> that we can have a PFSA that generates stochastic process with empty Q^+. In such a case, can we still get the PFSA structure back by studying the the set of causal states of the process? The answer is yes. 


Denote the collection all measures on Σ^ω, by ℳ_Σ, and let ν∈ℳ_Σ, the ε-ball of order d centered at ν is defined by

    B_d, ε(ν) = ν'∈ℳ_Σ| ∑_x∈Σ^dν'xΣ^ω - νxΣ^ω < ε..

In another words, B_d, ε is the collection of all measures that is no more than ε away from ν with respect to total variation distance over Σ^d.



Let Q be the set of causal states of a stochastic process 𝒫, a measure ν∈ℳ_Σ is an accumulation causal state of 𝒫 if 

    p_l, d, ε(ν) = ϕ^lμ_[x]∈ B_d, ε(ν)

satisfies p_d,ε(ν) = lim inf_l→∞p_l, d, ε(ν) > 0 for all d∈ℕ^+ and ε > 0. That is, a measure ν is an accumulation causal state if, no matter how large d is and how small ε is, the sum of probabilities of length-l sequences falling in B_d, ε(ν) does not vanish as l approaches infinity.

The collection of accumulation causal states is denoted by Q. Since p_d, ε(ν) is monotonically decreasing as d →∞ and ε→ 0, p(ν)=lim_d→∞lim_ε→0p_d, ε(ν) is well-defined. A measure ν with p(ν) > 0 is called an atomic accumulation causal state, and the collection of all atomic accumulation causal states is denoted by Q^+. 



Let ν∈ℳ_Σ, the translation of ν by σ for νσΣ^ω > 0, denoted by ν_σ, is the extension to  of the premeasure on the semiring  given by

    ν_σxΣ^ω = νσ xΣ^ω/νσΣ^ω.




Q^+ is closed under translation.



proof omitted.




Let 𝒫 be a stationary ergodic stochastic process with finite Q^+ and ∑_ν∈Q^+p(ν) = 1. Then the process generated by the PFSA G = Σ, Q^+, δ, with δ(ν, σ) = ν_σ and ν, σ = νσΣ^ω is exactly 𝒫. In fact, we have .𝐩_G|_ν = p(ν). 



proof omitted.



[Example <ref> Revisited]

We demonstrate that Q^+ of the process in Example <ref> has two elements, again by observation. We plot the cumulative probability density functions of ϕ_q(0) for each sequence length d = 10, 20, 30, 40 in Fig. <ref>. More specifically, for each fixed d, the x-coordinates of the dots are in Φ^d =ϕ_[x](0)|x∈Σ^d., while the y-coordinate of a dot with x-coordinate h_0∈Φ^d equals ϕ^dϕ_[x](0) ≤ h_0. We can see clearly that, the cumulative function converges to a step function with steps at 1/3 and 2/3 as d increases. The fact implies that ν∈Q^+ satisfies that ν0Σ^ω is either 1/3 or 2/3. We see from (<ref>) that the two measures in Q^+ are exactly q_-∞ and q_∞. Fig. <ref> also implies that that pq_-∞ = 1/3 and pq_∞ = 2/3, which is exactly the stationary distribution on the state set of S.





§ INFERENCE ALGORITHM OF PFSA


From the discussion in Sec. <ref>, we see that a stochastic process has a PFSA generator if either finitely many causal states get all the probability in the limit, as described in Sec. <ref>, or there exist finitely many measures in ℳ_Σ whose arbitrarily small neighborhoods are populated by almost all the causal states in the limit, as described in Sec. <ref>. The implication of these observations goes beyond the theory of PFSA, and guide us through the designing of inference algorithms of the model. In fact, a valid heuristic of the inference algorithm of PFSA would be to apply any clustering algorithm to the set of causal states corresponding to sequences up to a certain length, and use the center of the clusters to serve as estimates to the states. However, this primitive heuristic has a drawback since the cluster structure of [x]|x∈Σ^d may not be clear enough to facilitate a clustering algorithm. In order to get better estimates of the states, we need to fine tune our view into the set of causal states using the notion of ε-synchronizing sequence <cit.>.



 §.§ Epsilon-synchronizing Sequences

Before introducing ε-synchronizing sequence, we first introduce the concept of observation induced distributions over the state set. Let G be a PFSA, we know that the initial distribution over states is exactly the stationary distribution 𝐩_G. Let us assume that the first symbol generated by G is σ, denote by 𝐩_G(σ) the distribution over states after G producing σ, we have

    .𝐩_G(σ)|_q = 1/Z∑_q'|δq',σ = q(q', σ)𝐩_G|_q',

where 

    Z = ∑_q ∈ Q∑_q'|δq',σ = q(q', σ)𝐩_G|_q',

is the normalizer. 


Let x=σ_1…σ_n be a sequence observed, the distribution over states induced by x is defined inductively by

    .𝐩_Gσ_1…σ_i|_q = 1/Z∑_q'|δq',σ_i = q(q', σ)𝐩_Gσ_1…σ_i-1|_q',

where 

    Z = ∑_q ∈ Q∑_q'|δq',σ_i = q(q', σ)𝐩_Gσ_1…σ_i-1|_q',

for i = 1, …, n, with the base case 𝐩_G(λ) = 𝐩_G.



Let G be a strongly connected PFSA on state set Q over alphabet Σ. A sequence x∈Σ^⋆ is called an ε-synchronizing sequence for some ε>0 if there exists a q ∈ Q such that 𝐩_G(x) - 𝐞_q_∞ <ε, where 𝐞_q is the base probability vector with the entry indexed by q equalling 1. 


The reason that the ε-synchronizing sequences are important to inference is that [x_εx]| x∈Σ^d tends to have a much clearer cluster structure than [x]|x∈Σ^d for an ε-synchronizing sequence x_ε. 



 §.§  Algorithm

We give a brief review to the algorithm called  proposed in <cit.> in this section. By a sub-sequence, we mean a consecutive sub-sequence.

Let x∈Σ^⋆, the empirical symbolic derivative ϕ̂^x_y of a  sub-sequence y of x is given by 

    ϕ̂^x_y(σ) = number of sub-sequence yσ in x/number of sub-sequence y in x,

for all σ∈Σ.

Our inference algorithm is called  for Generator Extraction Using Self-Similar Semantics,  With the input of a long enough observed sequence x,  takes the following three steps to infer a PFSA: 



Step one: Approximate ε-synchronizing sequence: Calculate 

    𝒟^x_ε = ϕ̂^x_y | y is a sub-sequence of x with  |y|≤log_Σ1/ε.,

Then, select a sequence x_ε with ϕ̂^x_x_ε being a vertex of the convex hull of 𝒟^x_ε.



Step Two: Identify transition structure: For each state q, we associate a sequence identifier x_q ∈ x_εΣ^⋆, and a probability distribution d_q on Σ. We extend the structure recursively: Initialize the state set as Q = q_0, find x_q_0 and set d_q_0 = ϕ̂^x_x_q_0; Calculate the empirical symbolic derivative of x_qσ for each state q ∈ Q and σ∈Σ. If ϕ̂^x_x_qσ - d_q'_∞≤ε for some q' ∈ Q, then define δ(q,σ) = q'. However, if no such q' exists in Q, add a new state q' to Q, and define x_q' = x_qσ, and d_q' = ϕ̂^x_x_qσ. The process terminates when no more states can be added to Q. The inferred PFSA is the strongly connected component of the directed graph thus obtained.



Step Three: Identify observation probabilities: Initialize counter N^q_σ for each state q and symbol σ; choose an arbitrary initial state in the graph obtained in step two and run sequence x through it, if current state is q, and the next symbol from x is σ, then move to δ(q,σ), and add 1 to counter N^q_σ; finally, calculate the observation probability map by (q) = N^q_σ_σ∈Σ.




§ ENTROPY RATE AND KL DIVERGENCE




 §.§ Irreducibility of PFSA


We first discuss the concept of irreducibility for PFSA.

A PFSA G is irreducible if there is no other PFSA with strictly fewer number of states that generates the same stochastic process as G does.  

The definition of PFSA itself doesn't ensure irreducibility, as shown by example <ref>.

[Reducible PFSA]

In Fig. <ref>, we show two reducible PFSA. The PFSA on left generates the same process as the PFSA on the right Fig. <ref> does, while the PFSA on right generates the same procces as the PFSA on the right of Fig. <ref> does, but both with one more state than their respective irreducible versions.




Let a PFSA be specified by the quadruple Σ, Q, δ, and the measure μ_q be defined by μ_qxΣ^ω = (q, x). Two states q, q'∈ Q are equivalent if and only μ_q = μ_q'. 


We note that q_0^1 and q_0^2 in both PFSA in Fig. <ref> are equivalent. We also see that we can get the corresponding irreducible PFSA back by collapsing equivalent states to a single state.


A PFSA is irreducible if and only if it has no equivalent states. Furthermore, a irreducible PFSA is unique in the sense that, if two irreducible PFSA G_1 = Σ, Q_1, δ_1, _1 and G_2 = Σ, Q_2, δ_2, _2 generate the same stochastic process, there must be a one-to-one correspondence f:Q_1→ Q_2 such that fδ_1(q,σ)=δ_2(f(q), σ) and _2(f(q)) = _1(q).



The PFSA constructed on the set of persistent states Q^+ and the set of atomic accumulation states Q^+ are irreducible.





 §.§ Entropy Rate and KL Divergence



The entropy rate of a PFSA G is the entropy rate of the stochastic process G generates <cit.>. Similarly, the KL divergence of a PFSA G' from the PFSA G is the KL divergence of the process generated by the G' from that of G.
More precisely, we have the

    ℋ(G) = -lim_d→∞1/d∑_x∈Σ^dPr(x)log Pr(x),

and the KL divergence

    𝒟_KLGG' = lim_d→∞1/d∑_x∈Σ^dPr_G(x)logPr_G(x)/Pr_G'(x),

whenever the limits exist.



The entropy rate of a PFSA G = Σ, Q, δ, is given by 

    ℋ(G) = ∑_q∈Q𝐩_G(q)· h((q)),

where h(·) is the entropy of a probability distribution.



proof omitted.




Let G = Σ, Q, δ, and G' = Σ, Q', δ', ' be two PFSA, and let 𝐩_G(q, q') be the joint G-probability ofjoint state (q, q')
[The formal definition of joint G-probability needs long and technical derivation, which is outside the main focus of this paper. We can interpret 𝐩_G(q, q') as follows. Suppose we have a sample path generated by G, and we run the sample path on both G and G' (from arbitrary initial states) and calculate the frequency of the event “G is in state q and G' is in state q'” as a function of sequence length d. The frequency can be shown to converges as d approaches infinity and the limit is 𝐩_G(q, q').], 
then we have the KL divergence of G' from G is given by 

    ∑_(q, q')∈ Q×Q'p_G(q, q')D_KL(q)(q'),

where D_KL·· is the KL divergence between two probability distributions.



proof omitted.





 §.§ Log-likelihood



Let x∈Σ^d, the log-likelihood <cit.> of a PFSA G generating x is given by

    L(x, G) = -1/dlog Pr_G(x).




Let G and G' be two irreducible PFSA, and let x∈Σ^d be a sequence generated by G. Then we have

    L(x, G')→ℋ(G) + 𝒟_KLGG',

in probability as d→∞.

[sketch of proof]
We first notice that

    ∑_x∈Σ^dPr_G(x)logPr_G(x)/Pr_G'(x)
    
        =   ∑_x∈Σ^d-1∑_σ∈ΣPr_G(x)𝐩_G(x)._G|_σlogPr_G(x)𝐩_G(x)._G|_σ/Pr_H(x)𝐩_G'(x)._G'|_σ
    
        =   ∑_x∈Σ^d-1Pr_G(x)logPr_G(x)/Pr_H(x)
       + ∑_x∈Σ^d-1Pr_G(x)∑_σ∈Σ𝐩_G(x)._G|_σlog𝐩_G(x)._G|_σ/𝐩_G'(x)._G'|_σ_D^d.

By induction, we have 𝒟_KLGG' = lim_d→∞1/d∑_i=1^dD^i, 
and hence by Cesàro summation theorem, 𝒟_KLGG' = lim_d→∞ D^d whenever the limit exists.

Let x=σ_1σ_2...σ_n be a sequence generated by G. Let x^[i] be the truncation of x at the i-th symbols, we have

    -1/n∑_i=1^nlog𝐩_G'x^[i-1]._G'|_σ_i
    
        =   1/n∑_i=1^nlog𝐩_Gx^[i-1]._G|_σ_i/𝐩_G'x^[i-1]._G'|_σ_i_A_x, n-1/n∑_i=1^nlog𝐩_Gx^[i-1]._G|_σ_i_B_x, n.

Since the stochastic process G generates is ergodic, we have

    lim_n→∞A_x,n = lim_d→∞D^d = 𝒟_KLGG',

and lim_n→∞B_x,n = ℋ(G).


In this example we show the convergence of log-likelihood using the PFSA G on the left of Fig. <ref> and the PFSA H that is the induced subgraph on q_00, q_01, q_10, and q_11 in Fig. <ref>. We have

    ℋ(G)≈ 0.9710,  ℋ(H) ≈ 0.8069,

and

    𝒟_KLGH≈ 0.2266,  𝒟_KLHG≈ 0.2030.

Let us use G→x as the short hand for x is generated by G, we show in Fig. <ref> the log-likelihood of G producing a sequence x generated by G (top left), the log-likelihood of H producing a sequence x generated by H (top right), the log-likelihood of H producing a sequence x generated by G (bottom left), and the log-likelihood of G producing a sequence x generated by H (bottom right). We can clear see that the convergence of log-likelihood from the plots.








§ 


With the assumption of discrete-valued input, we first show in Sec. <ref> how to use log-likelihood convergence to define a pairwise distance between sequences. Because PFSA is a model for sequences on finite alphabet, continuous-valued input should first be quantized to discrete ones before being modeled by PFSA. So we discuss in Sec. <ref> ways of doing quantization and how their fitness can be evaluated. In Sec. <ref>, we compare  to  and fastDTW in both performance and efficiency.



 §.§ : Distance between Time Series


The way we calculate distance between two sequences is as follows. We first choose a set of PFSA 𝒢 = G_0,…, G_k as base, and the coordinate for a sequence x is defined to be 

    Lx, G_0, …, Lx, G_k,

where L(x, G) is the log-likelihood of G generating x, as defined in Sec. <ref>. The distance between a pair of sequences can then be any valid distance between their coordinates. 

For example, for the two numerical experiments in Sec.<ref>, we use 𝒢 that contains four PFSA

    0.    δq_0, 0 = q_0, δq_0, 1 = q_1, δq_1, 0 = q_0, δq_1, 1 = q_1,   
       q_-1=(.3, .7), q_1=(.7, .3),   
    
        1.    δq_0, 0 = q_0, δq_0, 1 = q_1, δq_1, 0 = q_1, δq_1, 1 = q_0,   
       q_-1=(.3, .7), q_1=(.7, .3),   
    
        2.    δq_0, 0 = q_1, δq_0, 1 = q_2, δq_1, 0 = q_2, δq_1, 1 = q_0,   
       δq_1, 0 = q_0, δq_2, 1 = q_1,   
       q_-1=(.3, .7), q_1=(.7, .3), q_2=(.6, .4),   
    
        3.    δq_0, 0 = q_0, δq_0, 1 = q_1, δq_1, 0 = q_2, δq_1, 1 = q_3,   
       δq_1, 0 = q_0, δq_2, 1 = q_1, δq_3, 0 = q_2, δq_3, 1 = q_3,   
       q_-1=(.3, .7), q_1=(.7, .3),   
       q_1=(.8, .2), q_3=(.2, .8).

and the distance between coordinates to be the total variation distance (l_0 distance). One thing we'd like to point out is that, we use these four PFSA as base primarily for simplicity. A better way of forming 𝒢, especially for supervised problems, is to use PFSA inferred by, for example,  proposed in Sec. <ref>, from the training set.



 §.§ Quantization of Continuous Sequence


The simplest approach to turn a continuous sequence to a symbolic one with alphabet size k is by choosing k-1 cut-off points p_1 < p_2 < ⋯ <p_k-1. With the additional assumption that p_0 = -∞ and p_k = +∞, we can replace a data point p in the continuous sequence with symbol i if p∈[p_i  p_i+1). We call the set of cut-off points a partition. The most common practice of choosing a partition is to apply the entropy maximization principle, in which the p_is are chosen so that we have as equal as possible numbers of data points in each interval. However, we can also perturb the cut-off points a little bit so that the distribution of symbols in the quantized sequences have smaller entropy. 

We can foresee the simple quantization above may work sub-optimally in the following two scenarios: first, the input sequences share a common trend; second, the input sequences have shifts and re-scalings. An example of the first case is the stock prices dataset, and examples for the second case are common among datasets of voltages that have different reference point or sound that have difference volumes. In some of these cases, We may want to first take derivatives (or detrend) and normalize the sequences to have empirical mean 0 and variance 1 before they are translated to symbolic sequences using a partition.

In our study of PFSA modeling, we try to use as many different quantization schemes as possible in the hope of exploring the dataset to a fuller extent. We do so by combining different decisions for whether to detrend and normalize and different choice of partitions to form a pool of quantization schemes. To record the parameters for a quantization scheme, we develop the following shorthand: 1) 𝖣𝖽, with 𝖽 being an non-negative integer, means detrending d times; 2) 0 means not to apply normalization, and 1, apply normalization; 3) [p_1 p_2 ⋯ p_k-1] means a k-part partition with cut-off points p_is. As an example, 𝖣1𝖭1[3.] means that we first detrend once, then normalize, and finally replace all data points that are less than 3. with a symbol 0, and those greater than or equal to 3., symbol 1.

To evaluate the quality of a quantization scheme, we first need to calculate the distance matrix D on the training set, with D_i,j being the distance between the i-th sequence and j-th sequence. Assuming that there are n sequences and l_i is the class label of the i-th sequence, we define the average inter-class distances to be

    s(D) = ∑_i= 1^n∑_j= 1^nδ_l_il_jD_i,j/∑_i= 1^n∑_j= 1^nδ_l_il_j,

and average intra-class distance to be

    d(D) = ∑_i= 1^n∑_j= 1^n1 - δ_l_il_jD_i,j/∑_i= 1^n∑_j= 1^n1 - δ_l_il_j,

where δ_ab = 1 if a = b and 0 if otherwise. Now we can evaluate the quality of the quantization scheme by the ratio r(D) = s(D)/d(D). The intuition behind this quality measurement is that a good quantization scheme should make the average distances between sequences from the same class small and the average distances between sequences from different classes large.



 §.§ , , and 


In this section we compare three similarity measurement of time series,  proposed in <cit.>,  proposed in <cit.>, and . We show using a synthetic dataset that  outperforms the other two algorithms in both performance and efficiency. The dataset we use contains two classes, each with 20 sequences of length 500 generated by the two PFSA, G and H, as shown in Fig. <ref>. The two PFSA are functionally different as shown by their their KL divergence. We have 𝒟_KLGH≈ 0.19677 and 𝒟_KLHG≈ 0.20756. For , we use the Python package  (<https://pypi.org/project/fastdtw/>) with default value 1 for radius. We use 2 for number of reruns for  because of its probabilistic nature. In Fig. <ref>(a-c) we show the heatmaps of the distance matrices given by the three algorithms. In Fig. <ref> (d), we show running time of the three algorithms on dataset constructed the same way as above but with sequence length ranging from 200 to 4000, with 200 increment. The computer we used to do the calculation has Intel E5-2680v4 2.4GHz CPU and 64GB memory. 







§ APPLICATIONS


In this section, we use  to study two real world problems. 


 §.§ Dataset 1: Motor Movement Imagery Dataset


This dataset is an excerpt of a dataset from  PysioNet<cit.>[<http://www.physionet.org/pn4/eegmmidb/>]. The dataset contains 64-channel 160Hz EEG recorded by the BCI2000 system<cit.>[<http://www.bci2000.org>] while subjects performed different motor imagery tasks. There are four tasks in the dataset and we focus on two of them: 

    
  : A target appears on either the left or the right side of the screen. The subject opens and closes the corresponding fist until the target disappears. Then the subject relaxes.
    
  : The same as the first task, except the subject imagines opening and closing the corresponding fist but doesn't really move.


For each subject, three 2-minute EEG recordings are taken for each task and we use the first two recordings to get the results listed below. During each recording, an object appears on the screen for 4 seconds and disappear for 4 second, and hence a subject supposedly moves or imagines to move his or her fists for 4 seconds, then rest for 4 seconds, and repeat. In Fig. <ref> and <ref>, the top two EEG recordings are for task  and the bottom two, task . We color the rest sections blue, while the movement/imaginary movement sections orange. 

For each subject we form a dataset with 56 sequence for each task, whose composition is detailed in Tab. <ref>. We drop the first section and the last section from each recordings since they tend to be more noisy.





The heatmap on the left of Fig. <ref> demonstrates the distance matrix calculated from channel 25 of the subject  on task . The quantization scheme of the sequences is 𝖣0𝖭1[-0.4526], with r(D) = .572. The heatmap on the right of Fig. <ref> demonstrates the distance matrix calculated from channel 21 of the same subject on task . The quantization scheme is 𝖣0𝖭0[-15.], with r(D) = .451. We can see that subject  has drastically different patterns in EEG between rest and (imaginary)movement sections, both from the wave and from the heatmaps of the distance matrices. The rest-movement difference of EEG is persistent across recordings as the distance between both the rest and movement sections from the first and second recordings are relatively insignificant.




The heatmap on the left of Fig. <ref> demonstrates the distance matrix calculated from channel 41 of subject  on task . The quantization scheme of the sequences is 𝖣0𝖭0[14.] with r(D) = .825. The heatmap on the right Fig. <ref> demonstrates the distance matrix calculated from channel 59 of the same subject on task . The quantization scheme is 𝖣1𝖭1[.4617], with r(D) = .830. We can see that subject  does not have significant difference in the EEG between rest and (imaginary) movement sections. Instead the subject seem to be in very different brain states during the two recordings for the same task.



 §.§ Dataset 2: User Identification from Walking Activity


This dataset from UCI machine learning repository <cit.> contains 22 subject walking along a predefined trail in the wild. Accelerometer measurements in x, y, z directions were taken with an Android smartphone placed in the chest pocket of each participant. The challenge is to identify a user using his or her pattern of motion. To form a training dataset for each subject, we get 10 sequence from the beginning of the measurement, each of 500 time steps long (each time step is about 0.03 second) and with 250 time step overlap between two consecutive sequences. There are 12 participants who has long enough measurement to form training datasets described above. 

In Fig. <ref>, <ref>, and <ref>, we plot the heatmaps of the best two distances for acceleration measurement in the x, y, and z directions, respectively. Although the distance matrices calculated on all 12 participants, we demonstrate the heatmaps for the first 5 participants for clarity.






From this example we can see that although it may be difficult to tell two participant apart by using measurement from one directions, using measurement from multiple channels combined with the freedom of our approach to choose quantization schemes customized for each individual channel, we may boost the classification performance. As for an example, while the distance in the x direction between the first two participants (sequences 0-9, and sequences 10-19) may not be great enough to tell them apart, their patterns of motion in the z direction do have much big distinction as indicated by the the first 20 by 20 diagonal block of the two heatmaps in Fig. <ref>.




§ CONCLUSION


In this paper, we propose a distance metric  between time series based on PFSA modeling and sequence likelihood divergence. We give a self-contained introduction to the mathematical foundation of PFSA as a time series model, the quantification of entropy rate and KL divergence of the stochastic process generated by PFSA, and finally, log-likelihood convergence of sample paths. We show how to infer PFSA from sequences and how to evaluate sequence likelihood divergence using log-likelihood convergence. We define the distance metric  using sequence likelihood divergence, and with the help of quantization algorithm of continuous data streams, we demonstrate how to apply  to the analysis of time series datasets arising from real world scenarios.

Possible future research effort includes 1) finding a better way to choose base PFSA for the  in unsupervised settings; 2) finding a base-free way to calculate  distance by, for example, inferring a PFSA model from each time series in the dataset.

siam



