{\bf e}gin{document}

\title {
{\LARGE \bf
	Symbolic Analysis-based Reduced Order 
\\ Markov Modeling of Time Series Data
	}
}

\author{Devesh K. Jha$^{a,1}$,  Nurali Virani$^{a,2}$, Jan Reimann$^{b}$, Abhishek Srivastav$^{c}$, Asok Ray$^{a, b}$
	\thanks {$^a$ Devesh K. Jha, N. Virani and Asok Ray are with Mechanical \& Nuclear Engineering Department, Pennsylvania State University, University Park, PA 16802, USA, and are partially supported by the U.S. Air Force Office of Scientific Research under Grant No. FA9550-15-1-0400; {\tt\{dkj5042,nnv105,axr2\}@psu.edu}}
	\thanks{$^b$ Jan Reimann and Asok Ray are with Department of Mathematics, Pennsylvania State University, University Park, PA, 16802, USA; {\tt\{jan.reimann,axr2\}@psu.edu}. Jan Reimann was partially supported by NSF Grant DMS-1201263}
	\thanks{$^c$ Abhishek Srivastav is with AI and Machine Learning Lab, GE Global Research Center, San Ramon, CA, USA; {\tt srivastav@ge.com}}
	\thanks{$^1$ Currently with Mitsubishi Electric Research Laboratories, Cambridge, MA 02139}
	\thanks{$^2$ Currently with AI and Machine Learning Lab, GE Global Research Center, Niskayuna, NY} \\        

 
\small \textbf{Keywords}: Symbolic Analysis, Markov Modeling, Order reduction, Combustion Instability
}

\maketitle

{\partial a}gestyle{plain}

{\bf e}gin{abstract}
This paper presents a technique for reduced-order Markov modeling for compact representation of time-series data. In this work, symbolic dynamics-based tools have been used to infer an approximate generative Markov model. The time series data are first symbolized by partitioning the continuous measurement space of the signal and then, the discrete sequential data are modeled using symbolic dynamics. In the proposed approach, the size of temporal memory of the symbol sequence is estimated from spectral properties of the resulting stochastic matrix corresponding to a first-order Markov model of the symbol sequence. Then, hierarchical clustering is used to represent the states of the corresponding full-state Markov model to construct a reduced-order (or size) Markov model with a non-deterministic algebraic structure. Subsequently, the parameters of the reduced-order Markov model are identified from the original model by making use of a Bayesian inference rule. The final model is selected using information-theoretic criteria. The proposed concept is elucidated and validated on two different data sets as examples. The first example analyzes a set of pressure data from a swirl-stabilized combustor, where controlled protocols are used to induce flame instabilities. Variations in the complexity of the derived Markov model represent how the system operating condition changes from a stable to an unstable combustion regime. In the second example, the  data set is taken from NASA's data repository for prognostics of bearings on rotating shafts. We show that, even with a very small state-space, the reduced-order models are able to achieve comparable performance and that the proposed approach provides flexibility in the selection of a final model for representation and learning.
\end{abstract}

\section{Motivation and Introduction}

Hidden Markov model (HMM) is a widely used statistical learning tool for modeling uncertain dynamical systems~\cite{B06}, where the associated temporal data are used to infer a Markov chain with unobserved states. In this setting, the learning task is to infer the states and the corresponding parameters of the Markov chain. In addition to HMM, several other nonlinear techniques have been proposed for Markov modeling of time-series data. Symbolic time-series analysis-based Markov modeling is a recently proposed technique~\cite{R04} where the states of a Markov chain are represented as a collection of words (i.e., symbol blocks, also referred to as memory words) of different lengths, which can be identified from the time-series data on a discrete space with finite cardinality~\cite{SS04, R04, MR14, CL13}. The symbols are created from the continuously varying time-series data by projecting the data to a set with finite cardinality. A common ground among all these tools of Markov modeling as discrete sequences, is that the Markov chain is induced by probabilistic representation of a deterministic finite state auotmaton (DFSA), often called probabilistic finite state automata (PFSA)~\cite{VTDCC05}. While the PFSA-based inference provides a consistent, deterministic graph structure for learning, the deterministic algebraic structure is generally not a very compact representation and may often lead to large number of states in the induced Markov model. To circumvent this problem attempts have been made to reduce the state-space by merging statistically similar states of the model~\cite{MR14}. The problem is, however, that as these models are constructed by partitioning of phase space of the dynamical system, merging states that are statistically similar leads to algebraic inconsistency. On the other hand, if the states are merged to preserve the algebraic consistency, it leads to statistical impurity in the final models (i.e., states which have different statistics could be merged together). Other approaches for state aggregation in Markov chains could be found in~\cite{GPKK15, V12, XSB14}. However, these papers do not consider inference of the Markov model from the data which may not be suitable for analysis of data-driven systems~\cite{D05}. 

The state space for Markov models, created by using symbolic analysis, increases exponentially with increase in memory or order of the symbolic sequence.  Estimating the right memory is critical for temporal modeling of patterns observed in the sequential data. However, some of the states may be statistically similar and thus merging them can reduce the size of state-space. This paper presents reduced-order Markov modeling of time-series data to capture temporal patterns, where we estimate the size of temporal memory of the symbolic data using the spectral properties of a PFSA whose states are words of length one~\cite{Srivastav2014, JSMR15}. The constraint of deterministic algebraic structure is not imposed by the end objective, but due to the choice of the data representation model. Thus we propose to merge the states and remove the constraint of deterministic algebraic properties associated with PFSA, where the states of the Markov chain are now collection of words from its alphabet of length estimated in the last step. This state aggregation induces a non-determinism in the finite state model. The parameters of the reduced-order Markov model are estimated by a Bayesian inference technique from the parameters associated with the higher-order Markov model. The final model for data representation is selected using information-theoretic criteria, and thus, we get a unique stopping point to terminate the state-merging procedure. We also present a bound on the distortion of the predictive capability of the models up on reduction in the size of the state-space. The final model obtained is a generative model for the data; however, some predictive capability is lost as we remove the deterministic algebraic structure of a DFSA.

The proposed technique of state merging is inspired by time-critical applications where it is imperative to arrive at a reliable decision quickly as the dynamics of the process being monitored is really fast. In such applications, there are strict constraints on accuracy as well as the time needed to come to a decision. In this paper, we illustrate the concepts using two different datasets. We discuss in detail the example of combustion instability which is a highly nonlinear and complex phenomena and results in severe structural degradation in jet turbine engines. Some good surveys on the current understanding of the mechanisms for the combustion instability phenomena could be found in~\cite{OAL15, SSDC03, CDSBM14, HY09, MBDSC12}. Active combustion instability control (ACIC) with fuel modulation has proven to be an effective approach for reducing pressure oscillations in combustors~\cite{BMJK06, BMH07}. Based on the work available in literature, one can conclude that the performance of ACIC is primarily limited by the large delay in the feedback loop and the limited actuator bandwidth ~\cite{BMJK06, BMH07}. Early detection of combustion instability can potentially alleviate the problems with delay in the ACIC feedback loop and thus possibly improve the performance. Some recent work for detection and prediction of combustion instabilities could be found in~\cite{JSR16, VJR16, SCRR16, NTS14, MS15}. While the results in these papers are encouraging, there is no interpretation of the expected changes in the data-driven model that could be observed during changes in the operating regime of the underlying process. In contrast to the work reported in literature, we have presented an overall idea of changes in the underlying stochastic model structure and parameters during the complex instability phenomenon.

\textbf{Contributions.} This paper presents a technique for Markov modeling of time series data using a PFSA with nondeterministic algebraic structure. Nondeterminism is induced by merging states of a PFSA with deterministic algebraic structure inferred from discrete sequential data, which in turn allows very compact representation of temporal data. In contrast to the approach in~\cite{MR14}, we present a method to use information-theoretic criteria to arrive at a consistent stopping criterion for model selection. The resulting reduced-order model has fewer parameters to estimate; this is turn leads to faster convergence rates and thus faster decisions during test (or operation). We also present a bound on the distortion in the predictive capability of the models due to state-space reduction using Hamming distance between the sequences generated by the original and final model. The algorithms presented in the paper are validated on two different datasets-- pressure data obtained from a swirl-stabilized combustor to monitor thermo-acoustic instability and a public data set for bearing prognostics. We show changes in the complexity of the pressure data as the process moves from stable to unstable through the transient phase which is then used to arrive at a criterion that provides perfect class separability. Apart from the results on Markov modeling, the results on combustion instability could be of independent interest in combustion community.

\section{Background and Mathematical Preliminaries}

Symbolic analysis of time-series data is a recent approach where continuous
sensor data are converted to symbol sequences via partitioning of the continuous
domain~\cite{SAX07, R04}. The dynamics of the symbols sequences are
then modeled as a probabilistic finite state automaton (PFSA), which is defined
as follows:

{\bf e}gin{definition}[PFSA]\label{defn:PFSA}
	A probabilistic finite state automaton (PFSA) is a tuple $G =( qSet, \alphaphabetSet, \delta, \bm{M})$ where
	{\bf e}gin{itemize}
		\item $qSet$ is a finite set of states of the automata;
		\item $\alphaphabetSet$ is a finite alphabet set of symbols $\alphaphabet \in \alphaphabetSet$;
		\item $\delta: qSet \times \alphaphabetSet \rightarrow qSet$ is the state transition function;
		\item {$\bm{M}: qSet \times \alphaphabetSet \rightarrow [0, 1]$ is the
			$\card{qSet}\times\card{\alphaphabetSet}$ emission matrix. The matrix
			$\bm{M} = [m_{ij}]$ is row stochastic such that $m_{ij}$ is the
			probability of generating symbol $\alphaphabet_{j}$ from state $q_{i}$}.
	\end{itemize}
\end{definition}
{\bf e}gin{remark} The PFSA defined above has a deterministic algebraic structure which is governed by the transition function $\delta$; thus a symbol emission from a particular state will lead to a fixed state. However, the symbol emissions are probabilistic (represented by the emission matrix). On the other hand, the transition function for a non-deterministic finite state automaton is given by a map, $\delta: qSet \times \alphaphabetSet \rightarrow 2^qSet$ where, $2^qSet$ denotes the power set of $qSet$ and includes all subsets of $qSet$. The idea is also presented in Figure~\ref{fig:NonDeterminism} where we show that the same symbol can lead to multiple states, however in a probabilistic fashion. This allows more flexibility in modeling at the expense of some predictive accuracy.
\end{remark}

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.9\textwidth]{figures/fig1_eps.eps}
	\caption{Graphical model showing non-determinism in a PFSA. The symbol $1$ emitted from state $q_1$ leads to different states with fixed probabilities indicating non-deterministic behavior.}
	\label{fig:NonDeterminism}
\end{figure}

For symbolic analysis of time-series data, a class of PFSAs called the
$D$-Markov machine have been proposed~\cite{R04} as a sub-optimal but
computationally efficient approach to encode the dynamics of symbol sequences as
a finite state machine.
{\bf e}gin{definition} (\textbf{$D$-Markov Machine}~\cite{R04, MR14}) \label{def:D-Markov} A $D$-Markov machine is a statistically stationary stochastic process $S= \cdots a_{-1}  a_{0} a_{1} \cdots $ (modeled by a PFSA in which each state is represented by a finite history of $D$ symbols), where the probability of occurrence of a new symbol depends only on the last $D$ symbols, i.e.,
	{\bf e}gin{equation}\label{eq:D-Markov}
		{\partial r}ob(s_n \mid \cdots  s_{n-D} \cdots s_{n-1} ) = {\partial r}ob(s_n \mid s_{n-D} \cdots s_{n-1}) \nonumber

	\end{equation}
	where $D$ is called the depth of the Markov machine.

\end{definition}
A $D$-Markov machine is thus a $D^{th}$-order Markov approximation of the discrete symbolic process. For most stable and controlled engineering systems that tend to forget
their initial conditions, a finite length memory assumption is reasonable. The $D$-Markov machine is represented as a PFSA and states of this PFSA are words over alphabet $\alphaphabetSet$ of length $D$ (or less);
the state transitions are described by a sliding block code of memory $D$
and anticipation length of one~\cite{LM95}.

For systems with fading memory it is expected that the predictive influence of a
symbol progressively diminishes. In this context, depth is defined as follows.
{\bf e}gin{definition}[Depth]\label{def:depth}
	Let $sSeq = s_{1}\dotss_{k}s_{k+1}s_{k+2}\dots$ be the
	observed symbol sequence  where each $s_{j}\in\alphaphabetSet\;\forall\; j \in
	\mathds{N}$. Then, the depth of the process generating $sSeq$ is defined as the
	length $D$ such that:
	{\bf e}gin{equation}\label{eq:depthTrueDefn}
		{\partial r}ob(s_{k}|s_{k - 1},\dots,s_{1}) = {\partial r}ob(s_{k}|s_{k - 1},\dots,s_{k - D})
	\end{equation}
	
\end{definition}
An accurate estimation of depth for the symbolic dynamical process is required for the precise modeling of the underlying dynamics of the discrete sequence. Next we introduce an information-theoretic metric which is used for merging the states of the Markov model later in next section.
{\bf e}gin{definition}[Kullback-Leibler Divergence]\label{def:KLD}~\cite{G90}
	The Kullback-Leibler (K-L) divergence of a discrete probability distribution $P$ from another distribution $\tilde{P}$ is defined as follows.
	{\bf e}gin{equation}
		D_{\textrm {KL}}(P\|\tilde{P})=\sum_{x\in X} {p}_X(x)\log\bigg(\frac{{p}_X(x)}{\tilde{p}_X(x)}\bigg) \nonumber
	\end{equation}
	It is noted that K-L divergence is not a proper distance as it is not symmetric. However, to treat it as a distance it is generally converted into symmetric divergence as follows, $d(P,\tilde{P})= D_{\textrm {KL}}(P\|\tilde{P})+D_{\textrm {KL}}(\tilde{P}\|P)$. This is defined as the K-L distance between the distributions $P$ and $\tilde{P}$.
\end{definition}

This distance is used to find out the structure in the set of the states of the PFSA-based Markov model whose states are words, over the alphabet of the PFSA, of length equal to the depth estimated for the discretized sequence.

\section{Technical Approach}
In this section, we present the details of the proposed approach for inferring a Markov model from the time series data. As discussed earlier, the first step is the discretization of the time-series data to generate a discrete symbol sequence. While it is possible to optimize the symbolization of time-series using some optimization criterion, we do not discuss such a technique here. The data is discretized using the unbiased principle of entropy maximization of the discrete sequence using maximum entropy partitioning (MEP)~\cite{RR06}. The proposed approach for Markov modeling then consists of the following four critical steps
{\bf e}gin{itemize}
	\item Estimate the approximate size of temporal memory (or order) of the symbol sequence.
	\item Cluster the states of the high-order Markov model.
	\item Estimate the parameters of the reduced-order Markov model (i.e., the transition matrix).
	\item Select the final model using information theoretic scores (described below, Section~\ref{subsec:MDL}).
\end{itemize}
Memory of the discrete sequence is estimated using a recently introduced method based on the spectral analysis of the Markov model with depth $1$. induced by a PFSA~\cite{Srivastav2014, JSMR15}. It is noted that these steps are followed during training to estimate the approximate model for data and during test, the parameters are estimated for the reduced-order model.

The key ideas behind these steps are explained in the next section.

\subsection{Estimation of Reduced-Order Markov Model}\label{subsec:reducedorder}
Depth $D$ of a symbol sequence has been redefined in~\cite{Srivastav2014}
as the number of time steps after which probability of current symbol is
independent of any past symbol i.e.:
{\bf e}gin{equation}\label{eq:depthDefn}
	{\partial r}ob(s_{k}|s_{k - n}) = {\partial r}ob(s_{k})  \ \forall n>D
\end{equation}
Note that dependence in the proposed definition (eq.~\ref{eq:depthDefn}) is
evaluated on individual past symbols using ${\partial r}ob(s_{k}|s_{k - n})$ as
opposed to the assessing dependence on words of length $D$ using  $
{\partial r}ob(s_{k}|s_{k - 1},\dots,s_{k - D})$. It is shown that if
the observed process is {\it forward causal} then observing any additional
intermediate symbols  $s_{k- 1},\dots,s_{k -  n + 1}$ cannot induce a
dependence between $s_{k }$ and $s_{k - n}$ if it did not exist on
individual level~\cite{Srivastav2014}.

Let $\bm{\Pi} = [\pi^{(1)}_{ij}]$ be the one-step transition probability
matrix of the PFSA $G$ constructed from this symbol sequence i.e.
{\bf e}gin{equation}\label{eq:stateTransition1step}
	\bm{\Pi} = {\partial r}ob(s_{k}|s_{k - 1})
\end{equation}
Then using the distance of the transition matrix after steps from the stationary
point, depth can be defined as a length $D$ such that
{\bf e}gin{equation}\label{eq:depthTrace}
	\abs{\trace{\bm{\Pi}^{n}} - \trace{\bm{\Pi}^{\infty}}} \leq \sum_{j = 2}^{J}
	\abs{\lambda_{j}}^n < \epsilonsilon \ \forall n > D
\end{equation}
where $J$ is number of non-zero eigenvalues of $\bm{\Pi}$. Thus, the depth $D$ of the symbol sequence is estimated for a choice of $\epsilonsilon$ by estimating the stochastic matrix for the one-step PFSA. Next, another pass of data is done to estimate the PFSA parameters whose states are words over $\alphaphabetSet$ of length $D$, i.e., $\bm{\Pi} = {\partial r}ob(s_{k}|s_{k - 1},\dots, s_{k - D})$. It is noted that this step is critical for modeling accuracy.

The states of the reduced-order Markov model are then estimated by partitioning the set of words over $\alphaphabetSet$ of length $D$ estimated in the last step. This is done by using an agglomerative hierarchical clustering approach. The advantage of using the hierarchical clustering approach is that it helps visualize the structure of the set of the original states using an appropriate metric. Agglomerative hierarchical clustering is a bottom-up clustering approach~\cite{XW05} that generates a sparse network (e.g., a binary tree) of the state set $qSet$ (where $|Q|=|\alphaphabetSet|^D$) by successive addition of edges between the elements of $qSet$. Initially, each of the states $q_1,q_2,\dots,q_n$  is in its own cluster $C_1,C_2,\dots, C_n$ where $C_i\in \bm{M}hcal{C}$, which is the set of all clusters for the hierarchical cluster tree. The distance between any two states in $qSet$ is measured using the K-L distance between the symbol emission probabilities conditioned on them, i.e.,
{\bf e}gin{align}\label{eq:kldistance}
	d(q_i,q_j)& =  D_{\textrm{KL}}({\partial r}ob(\alphaphabetSet|q_i)\|{\partial r}ob(\alphaphabetSet|q_j))\nonumber \\
	&+D_{\textrm{KL}}({\partial r}ob(\alphaphabetSet|q_j)\|{\partial r}ob(\alphaphabetSet|q_i))
\end{align}
where the terms on the right have the following meaning.
{\bf e}gin{align}
& D_{\textrm{KL}}({\partial r}ob(\alphaphabetSet|q_i)\|{\partial r}ob(\alphaphabetSet|q_j)) \nonumber \\
& = \sum_{s \in \alphaphabetSet}{\partial r}ob(s|q_i)\log\bigg( \frac{{\partial r}ob(s|q_i)}{{\partial r}ob(s|q_j)}\bigg) \nonumber
\end{align}
In terms of the distance measured by eq.~\eqref{eq:kldistance}, the pair of clusters that are nearest to each other are merged and this step is repeated till only one cluster is left. The tree structure displays the order of splits in the state set of the higher-order Markov model and is used to aggregate the states close to each other. For clarification of presentation, we show an example of a Markov chain with $27$ states and $3$ symbols on a simplex plane in Figure~\ref{fig:Simplexcluster}, where each \mbox{col}or{red} red pentagon \mbox{col}or{black} on the simplex represents one row of the symbol emission matrix. The hierarchical clustering is used to find the structure of the state set on the simplex place using the K-L distance. The set of states clustered together could be obtained based on the number of final states required in the final Markov model.

The overall algorithm is presented as a pseudo-code in Algorithm~\ref{algorithm:Modeling}. This algorithms is used to find the parameters of the models during training. The parameters during test are estimated using the clustering map $f_{N_{\textrm{max}}}$ and is further discussed in next section. In the later sections we show how an information theoretic criterion could be used to select the appropriate model to terminate the state merging algorithm or select a final model from the set of reduced-order models. Through numerical experiments using two different data-sets we also illustrate the main motivation of this work that although the right memory is required for accurate modeling of the symbolic process, the state-space not necessarily consist of all words corresponding to the estimated memory and we can achieve sufficiently-high predictive accuracy even with a smaller state-space. We are able to achieve this trade-off between the model complexity and predictive modeling accuracy using the information-theoretic criteria.

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/simplex_fig1_1.eps}
	\caption{The symbol emission probabilities for a Markov chain with $3$ symbols are shown on a simplex. Symmetric K-L distance is used to find the structure in the state-set in the information space and the states are clustered based on the revealed structure.}
	\label{fig:Simplexcluster}
\end{figure}

{\bf e}gin{algorithm}[h] \small
\SetKwInput{KwIn}{Input}
\SetKwInput{KwOut}{Output}

	\KwIn {The observed symbol sequence $sSeq=\{\dots s_1s_2s_3\dots|s_i \in \alphaphabetSet\}$}
	
	\KwOut {The final Markov model, $\bm{M}hcal{M}=(\tilde{qSet},\tilde{\bm{M}},\tilde{\bm{\Pi}})$}	
	
	Estimate the $\bm{\Pi}$ matrix for 1-step Markov model using frequency counting with an uniform prior\; 
	
	Estimate the size of temporal memory, $D(\epsilonsilon)$ for $sSeq$ using equation~\eqref{eq:depthTrace}\;
	
	Estimate $\bm{M}$ and $\bm{\Pi}$ for the $D(\epsilonsilon)$-Markov model using frequency counting with an uniform prior\;
	
	$\bm{M}hcal{C}_{\midqSet\mid}=\{q_i\mid q_i \in qSet\}$\;
	
	\For{$i = \midqSet\mid-1,\dots,1$}{{find distinct clusters $A,B \in \bm{M}hcal{C}_{i+1}$  minimizing $d(A\cup B)$\;}{$\bm{M}hcal{C}_i:=(\bm{M}hcal{C}_{i+1}\setminus\{A,B\})\cup \{A\cup B\}$}}
	\Return{$\bm{M}hcal{C}_1,\dots,\bm{M}hcal{C}_{\midqSet\mid}$ and $f_i:qSet \rightarrow \bm{M}hcal{C}_i$} $\forall i \in \{1,\dots,\midqSet\mid\}$
	
	Calculate the parameters of reduced model using $\tilde{qSet}=\bm{M}hcal{C}_{N_{\textrm{max}}}$, $f_{N_{\textrm{max}}}$ and equations~\eqref{eq:ParameterEst} through~\eqref{eq:pq1q2}\;
	
	Calculate the Log-likelihood for models with Equation~\eqref{eq:depth_loglikelihood}\;
	
	The final model is selected using the AIC or BIC criteria explained in Section~\ref{subsec:MDL}\;

	\caption{Reduced Order Markov Modeling}
	\label{algorithm:Modeling}
\end{algorithm}

\subsection{Parameter Estimation of the Reduced-Order Markov Model}\label{subsec:DBN}
The parameters of the Markov model obtained after clustering the states of the original PFSA with $|\alphaphabetSet|^D$ states is obtained using a Bayesian inference technique using the parameters estimated for the PFSA. In this proposed approach, the state transition matrix $\bm{\Pi}$, the emission matrix $\bm{M}$, and the state probability vector $pVec$ of the original PFSA model $G$ are available, along with the deterministic assignment map $f:qSet \rightarrow \widetilde{qSet}$ of the state in $qSet$ (i.e., state set of original model) to one of the state in $\widetilde{qSet}$ (i.e., state set of the reduced order model).
Since the reduced order model can represented by the tuple $\widetilde{G} = (\widetilde{qSet}, \widetilde{\bm{\Pi}})$, where $\widetilde{\bm{\Pi}} = [\tilde{\pi}_{ij}]$ is the state transition matrix, we employ a Bayesian inference technique to infer the individual values of transition probabilities $\tilde{\pi}_{ij} = {\partial r}ob(\tilde{q}_{k+1} = j \mid \tilde{q}_{k} = i)$ for all $i, j \in \widetilde{qSet}$.

Let $qVar_{k}$ be the random variable denoting the state of PFSA model at some time step $k \in \mathds{N}$ and $sVar_{k}$ denotes the symbol emitted from that state, this probabilistic emission process is governed by the emission matrix $\bm{M}$. The state of the reduced order model is obtained from a deterministic mapping of the state of the PFSA model, thus the state of this model is also a random variable, which is denoted by $\widetilde{qVar}_{k} = f(qVar_{k})$. The Bayesian network representing the dependencies between these variables is shown in the recursive as well as unrolled form in the Figure~\ref{fig:DBN}.

{\bf e}gin{figure*} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/DBN2.eps}
	\caption{Graphical models representing the dependencies between the random variables}
	\label{fig:DBN}
\end{figure*}
The conditional density ${\partial r}ob(\widetilde{qVar}_{k} = \tilde{q} \mid qVar_{k} = q)$ can be evaluated by checking if state $q$ belongs to the state cluster $\tilde{q}$ and assigning the value of 1 if true, else assign it the value of 0. Since we know that $\widetilde{qSet}$ partitions the set $qSet$, the conditional density is well-defined. Thus, it can be written as
{\bf e}gin{align}\label{eq:pc1q1}
	{\partial r}ob(\widetilde{qVar}_{k} = \tilde{q} \mid qVar_{k} = q) = \indicator{\tilde{q}}{q},
\end{align}
where $\operatorname{I}$ is the indicator function with $\indicator{\tilde{q}}{q} = 1$, if element $q$ belongs to the set $\tilde{q}$, else it is $0$. The derivation of the Markov model ${\partial r}ob(\widetilde{qVar}_{k+1}\mid \widetilde{qVar}_{k})$ using ${\partial r}ob(qVar_{k+1}\mid qVar_{k})$, stationary probability vector $pVec$, and assignment map $f$ is shown ahead.
{\bf e}gin{align}\label{eq:ParameterEst}

	&{\partial r}ob(\widetilde{qVar}_{k+1}\mid \widetilde{qVar}_{k}) = \sum_{q \in qSet} {\partial r}ob(\widetilde{qVar}_{k+1}, qVar_{k+1} = q \mid \widetilde{qVar}_{k})  \\
	& \text{(Marginalization)}\notag\\

	&{\partial h}antom{{\partial r}ob(\widetilde{qVar}} = \sum_{q \in qSet} {\partial r}ob(qVar_{k+1} = q \mid \widetilde{qVar}_{k}) {\partial r}ob(\widetilde{qVar}_{k+1} \mid qVar_{k+1} = q)\\
	& \text{(Factorization using Figure~\ref{fig:DBN})}\notag\\
	&{\partial h}antom{{\partial r}ob(\widetilde{qVar}} = \sum_{q \in qSet} {\partial r}ob(qVar_{k+1} = q \mid \widetilde{qVar}_{k}) \indicator{\widetilde{qVar}_{k+1}}{q} \\
	& \text{(using~\eqref{eq:pc1q1})}\notag\\
	&{\partial h}antom{{\partial r}ob(\widetilde{qVar}} = \sum_{q \in \widetilde{qVar}_{k+1}} {\partial r}ob(qVar_{k+1} = q \mid \widetilde{qVar}_{k}) \label{eq:pc2c1}.
\end{align}
We can obtain ${\partial r}ob(qVar_{k+1} \mid \widetilde{qVar}_{k})$ from Bayes' rule as
{\bf e}gin{align}\label{eq:pq2c1}
	{\partial r}ob(qVar_{k+1} \mid \widetilde{qVar}_{k}) = \dfrac{{\partial r}ob(\widetilde{qVar}_{k} \mid qVar_{k+1}){\partial r}ob(qVar_{k+1})}{\sum_{q \in qSet}{\partial r}ob(\widetilde{qVar}_{k} \mid qVar_{k+1}=q){\partial r}ob(qVar_{k+1} = q)}.
\end{align}
Following the steps to obtain~\eqref{eq:pc2c1}, we also derive
{\bf e}gin{align}\label{eq:pc1q2}
	{\partial r}ob(\widetilde{qVar}_{k} \mid qVar_{k+1}) = \sum_{q \in \widetilde{qVar}_{k}} {\partial r}ob(qVar_{k} = q \mid qVar_{k+1}) .
\end{align}
We can obtain ${\partial r}ob(qVar_{k} \mid qVar_{k+1})$ from Bayes' rule as
{\bf e}gin{align}\label{eq:pq1q2}
	{\partial r}ob(qVar_{k} \mid qVar_{k+1}) = \dfrac{{\partial r}ob(qVar_{k+1} \mid qVar_{k}){\partial r}ob(qVar_{k})}{\sum_{q \in qSet}{\partial r}ob(qVar_{k+1} \mid qVar_{k}=q){\partial r}ob(qVar_{k} = q)}.
\end{align}
Note that, for the distribution ${\partial r}ob(qVar_{k})$ and ${\partial r}ob(qVar_{k+1})$, we use the stationary probability $pVec$. Using the equations \eqref{eq:pc2c1}, \eqref{eq:pq2c1},
\eqref{eq:pc1q2}, and \eqref{eq:pq1q2} together, one can easily obtain the desired state transition matrix $\widetilde{\bm{\Pi}}$ of the reduced order model. Once the state cluster set $\widetilde{qSet}$ and state transition matrix $\widetilde{\bm{\Pi}}$ are available, the reduced order model is completely defined.

\subsection{Model Selection using information theoretic criteria}\label{subsec:MDL}
In this section, we describe the model selection process during the underlying state merging process for model inference. We compute ``penalized'' likelihood estimates for different models. Then, the model with the lowest score is selection as the optimal model.

The (unpenalized) log-likelihood of a symbol sequence $sSeq$ given a Markov model $G$ is computed as follows:
{\bf e}gin{align}\label{eq:loglikelihood_Markov}

\mathcal{L}(sSeq | G) \operatorname{co}ng
\sum_{k=1}^{N}\log {\partial r}ob\left( s_{k} | q_k \right)
\end{align}
where the effects of the initial state are ignored because they become negligible for long statistically stationary symbol sequences. It is noted that with a finite symbol sequence, the log-likelihood is always finite. Furthermore, with the Markov models considered in this paper, the sum is simplified to the following form.
{\bf e}gin{align}\label{eq:loglikelihood_DMarkov}

\mathcal{L}(sSeq | G) \operatorname{co}ng
\sum_{k=D+1}^{N}\log {\partial r}ob\left( s_{k} | s_{k-1},\dots,s_{k-D} \right)
\end{align}

As discussed earlier, the states are merged using hierarchical clustering and thus, for every desired number of final states we get the deterministic map $f_{N_{\textrm{max}}}$ which determines how the original states are partitioned using the hierarchical clustering. This map is known for every terminal number of states and thus, we can find the log-likelihood of the symbol sequence using the following relationship.
{\bf e}gin{align}\label{eq:depth_loglikelihood}

\mathcal{L}(sSeq | \tilde{G}) \operatorname{co}ng
\sum_{k=D+1}^{N}\log {\partial r}ob\left( s_{k} | \tilde{q}_k=f_{N_{\textrm{max}}}(q_k) \right)
\end{align}
where, $\tilde{q}_k$ is the state of the reduced model and $q_k$ is the state of the original full-order model.

In the next step of the model selection process, a ``complexity penalty'' is added to the log-likelihood estimates, thereby balancing goodness of fit against the complexity of the model (and hence trying to prevent overfitting). We apply two widely-used such model selection functions, namely the Akaike information criterion (AIC)~\cite{Akaike:1974a} and the Bayesian information criterion (BIC)~\cite{Schwarz:1978a}:
{\bf e}gin{enumerate}
\item $\bm{M}hcal{M}_{\textrm{BIC}}=-2\mathcal{L}(sSeq | \tilde{G})+K\log(N)$, where $K$ is the number of free parameters and $N$ is the number of observations.
\item $\bm{M}hcal{M}_{\textrm{AIC}}=-2\mathcal{L}(sSeq | \tilde{G})+2K$, where $K$ is the number of free parameters.
\end{enumerate}
The number of free parameters to be estimated from the data is the parameters of the symbol emission parameters, i.e., $K=\mid\alphaphabetSet\mid \mid\tilde{qSet}\mid$. It is noted that this allows model selection for individual symbol sequences. The criterion here allows a terminal condition for state merging; however, different symbol sequences can have different models. The model with the minimum score is selected as the best model. Through the results presented in next sections we illustrate the fact that most of the temporal and predictive capabilities can be preserved for the models with a very small number of states when compared to the original model.
{\bf e}gin{remark}
	The final Markov model is a finite depth approximation of the original time-series data. However, compared to the PFSA-based D-Markov machines in~\cite{R04, MR14}, the current aggregated model has a non-deterministic algebraic structure, i.e., the same symbol emissions from a state can lead to different states. While this leads to some loss in predictive capability as compared to the models in~\cite{R04, MR14}, this allows us to compress the size of the model as per the requirement at hand. This allows faster convergence rates for the symbol emission probabilities as we only require fewer parameters to estimate from data, which might lead to faster decisions during testing.
\end{remark}
In the rest of the paper, we will present a Hamming distance-based bound for distortion in the predictive capabilities of reduced models and demonstrate the utility of these models in practical problems of fault/anomaly detection from time-series data.
\section{Analysis of the Proposed Algorithm}\label{sec:analysis}
In this section, we will present a bound on the distortion of the model due to the reduction of state-space of the Markov model using Hamming distance between two symbol sequences. We first present the Pinsker's inequality~\cite{G90} which relates the information divergence with the variational distance between probability measures defined on arbitrary spaces. This is followed by another theorem which can be used to derive Hamming distance bounds using the informational divergence.
{\bf e}gin{theorem}[Pinsker's inequality]~\cite{G90}
Let $P$ and $Q$ be two probability distributions on a measurable space $(\bm{M}hds{X},\Sigma)$. Then, the following is true
{\bf e}gin{eqnarray}
d_{TV}(P,Q)\leq \sqrt{\frac{1}{2}D_{\textrm {KL}}(P\|Q)}
\end{eqnarray}
where $d_{TV}(P,Q)=\sup\limits_{A\in \Sigma} \{|P(A)-Q(A)|\}$ is the total variation distance.
\end{theorem}
{\bf e}gin{theorem}~\cite{M96}
Let $\bm{M}hds X$ be a countable set and let us denote by $x^n$ the sequence $(x_1,x_2,\dots,x_n)\in \bm{M}hds{X}^n$. Let $q^n$ be a Markov measure on $\bm{M}hds X^n$, that is, $q(x^n)=q(x_1){\partial r}od\limits_{i=2}^n q_i(x_i|x_{i-1})$. Then for any probability measure $p^n$ on $\bm{M}hds X^n$, the following is true
{\bf e}gin{eqnarray}
\bar{d}(p^n,q^n)\leq \bigg[\frac{1}{2n}D_{\textrm {KL}}(p^n\|q^n)\bigg]^{1/2}
\end{eqnarray}
where, $\bar{d}$ denotes the normed Hamming distance on $\bm{M}hds X^n \times \bm{M}hds X^n:$
{\bf e}gin{equation}
\bar{d}(x^n,y^n)=n^{-1}\sum\limits_{n=1}^n d(x_i,y_i),
\end{equation}
where $d(x_i,y_i)=1$ if $x_i\neq y_i$ and $0$ otherwise. The $\bar{d}$-distance between $p^n$ and $q^n$ is
{\bf e}gin{eqnarray}
\bar{d}(p^n,q^n)=\min E  \bar{d}(\hat{X}^n,X^n),
\end{eqnarray}
where $\min$ is taken over all joint distributions with marginals $p^n=\bm{M}hrm{dist} \hat{X}^n$ and $q^n=\bm{M}hrm{dist} {X}^n$ and $E$ denotes the expectation operator.
\end{theorem}

The above theorem provides us a way to bound Hamming distance between sequences generated by two different distributions. Thus, using the above theorem, we find a bound on the Hamming distance between the symbol sequences generated by the reduced-order Markov model and the original model by estimating the K-L distance between the measure on symbol sequences induced by these models. An approximate estimate of the K-L distance between the original and a reduced model could be expressed and estimated as shown in the following.

 Let the original D-Markov model be denoted by $\bm{M}hcal{M}$ and the reduced-order model by $\hat{\bm{M}hcal{M}}$. The Markov measure on the probability space $(S^n,\bm{M}hcal{E},P)$ where the set $S^n$ consists of sequences of length $n$ from an alphabet $\alphaphabetSet$ could be estimated using the symbol emission probabilities. More explicitly, the Markov measure of a sequence $S_n$ on $S^n$ induced by $\bm{M}hcal{M}$ is given by $P_{\bm{M}hcal{M}}(S_n)={\partial r}ob(q_1){\partial r}od\limits_{i=D+1}^n{\partial r}ob(s_i\mid q_i)$ (where $D$ is the depth of the model). Then, the K-L divergence between $\bm{M}hcal{M}$ and $\hat{\bm{M}hcal{M}}$ is given by the following expression.
{\bf e}gin{equation}~\label{HammingBound}
D_{\rm {KL}} (P^n_{\bm{M}hcal{M}}\|P^n_{\hat{\bm{M}hcal{M}}})=\sum\limits_{S_n \in S^n}P_{\bm{M}hcal{M}}(S_n) \log\bigg(\frac{P_{\bm{M}hcal{M}}(S_n)}{P_{\hat{\bm{M}hcal{M}}}(S_n)}\bigg)
\end{equation}
Then, the above expression can be simplified as follows.
{\bf e}gin{equation}
\log\bigg(\frac{P_{\bm{M}hcal{M}}(S_n)}{P_{\hat{\bm{M}hcal{M}}}(S_n)}\bigg)=\sum\limits_{i=D+1}^n\log({\partial r}ob(s_i\mid q_i))-\log({\partial r}ob(s_i\mid \hat{q}_i)), \nonumber
\end{equation}
where, $\hat{q}$ is the merged state and $q$ is the original state. Then the expression on the right could be further bounded using the Lipschitz constant for the logarithm function and under the assumption that $\log({\partial r}ob(s_j\mid q_i))\neq 0$ $\forall q_i \in qSet$ and all $s_j \in \alphaphabetSet$.
{\bf e}gin{align}
&\sum\limits_{i=D+1}^n\log({\partial r}ob(s_i\mid q_i))-\log({\partial r}ob(s_i\mid \hat{q}_i)) \label{eqn:logsum}\\
&\leq \sum\limits_{i=D+1}^n(\frac{{\partial r}ob(s_i\mid q_i)-{\partial r}ob(s_i\mid \hat{q}_i)}{{\partial r}ob(s_i\mid q_i)})\label{eqn:lipschitz}\\
&\leq (n-D-1)\kappa
\end{align}
where, $\kappa=\max\limits_{q\in Q, s\in \alphaphabetSet}\frac{{\partial r}ob(s\mid q)-{\partial r}ob(s\mid \hat{q})}{{\partial r}ob(s\mid q)}$. In the above inequalities, equation~\eqref{eqn:lipschitz} is obtained from equation~\eqref{eqn:logsum} by using the observation that ${\partial r}ob(s_i\mid \hat{q}_i)={\partial r}ob(s_i\mid {q}_i)+\eta$, where $\eta$ is the perturbation in the symbol emission probability from $q_i$ when it is clustered into a new state $\hat{q}_i$. Hence, the K-L distance in equation~\eqref{HammingBound} could be bounded by the following term.
{\bf e}gin{align}
D_{\rm {KL}} (P^n_{\bm{M}hcal{M}}\|P^n_{\hat{\bm{M}hcal{M}}})& \leq \sum\limits_{S_n \in S^n} P_{\bm{M}hcal{M}}(S_n) (n-D-1)\kappa \nonumber \\
& = (n-D-1)\kappa \sum\limits_{S_n \in S^n} P_{\bm{M}hcal{M}}(S_n) \nonumber \\
& =(n-D-1)\kappa
\end{align}
 Thus, a uniform bound on the Hamming distance between the original and the final model could then be obtained as follows,
{\bf e}gin{equation}\label{eqn:bound}
\bar{d}(P_\bm{M}hcal{M}(S_n),P_{\bm{M}hcal{\hat{M}}}(S_n))\leq \sqrt{\frac{(n-D-1)\kappa}{2n}}
\end{equation}
The above inequality thus, allows us to compare models with different state-space based on the predictive accuracy of a reduced model when compared to the original model. As compared to the earlier information theoretic criteria, which were based on the efficiency of data compression by different models, the inequality in~\eqref{eqn:bound} allows to compare them based on their symbol emission statistics and thus, is computationally efficient. It is possible to find a rather tighter bound in an expected sense by using the stationary distribution of the two Markov chains to find an expected bound on Hamming distance. However, finding the same is left as an exercise for future work. Using the above bound for selection of models could be more efficient than the information theoretic metrics (as it can estimated by using the symbol emission probabilities instead of the penalized likelihoods); however, finding a penalized version of the bound for model selection is also left as a future exercise.

\section{Description of Experimentation and Data Sets}\label{sec:experiment}
In this section, we  briefly describe the two different data-sets which have been used in this paper to illustrate and validate the proposed concepts. Specifically, we will describe the experiments done at Penn State to investigate instability in lean-premixed combustion and another benchmark data-set for anomaly detection in bearings. An important point to be noted here is that the numerical experiments we present in the following sections is to justify the fact that the reduced-order models obtained by the proposed algorithms are able to achieve the trade-off between predictive accuracy and model complexity. Further results for classification and anomaly detection are to illustrate that this proposed approach of model learning can still achieve good performance for machine learning objectives of class separability and anomaly detection.
\subsection{Combustion}\label{subsec:combustionexperiment}
A swirl-stabilized, lean-premixed, laboratory-scale combustor was used to perform the experimental study. Tests were conducted at a nominal combustor pressure of 1 atm over a range of operating conditions, as listed in Table~\ref{tab:par}.

{\bf e}gin{table}[!hbp]
	\centering
	\caption{Operating conditions}
	{\bf e}gin{tabular}{c|c}
		\hline
		\textbf{Parameters} & \textbf{Value} \\
		\hline
		Equivalence Ratio & 0.525, 0.55, 0.60, 0.65\\
		\hline
		Inlet Velocity & 25-50 m/s i m/s increments \\
		\hline
		Combustor Length & 25-59 inch in 1 inch increments\\
		\hline
	\end{tabular}
	\label{tab:par}
\end{table}

 In each test, the combustion chamber dynamic pressure and the global OH and CH chemiluminescence intensity were measured to study the mechanisms of combustion instability. The measurements were made simultaneously at a sampling rate of 8192 Hz~(per channel), and data were collected for 8 seconds, for a total of 65536 measurements~(per channel). A total of $780$ samples of data were collected from all the tests where in every test the combustion process was driven from stable to unstable by changing either the equivalence ratio, ${\partial h}i$. However, as the accurate model of the process is not available, an accurate label of transition of the process to unstable phase is not available. It is noted that the data consists the behavior of the process over a large number of operating condition and thus provides a rich set of data to test the efficacy of the algorithm in detecting classes irrespective of the underlying operating conditions.
 \subsection{Bearing Prognostic Data} This test data has been picked from NASA's prognostics data repository~\cite{NASAPHM, TMZT12}. A detailed description of the experiments could be found in~\cite{QLLY06}. The bearing test rig hosts four test bearings on one shaft which is driven by an AC motor at a constant speed. A constant force is applied on each of the bearings and accelerometer data is collected at every bearing at a sampling rate of \SI{20}{\kilo\hertz} for about \SI{1}{\second}. The tests are carried for $35$ days until a significant amount of debris was found in the magnetic plug of the test bearing. A defect in at least one of the bearings is found at the end of every test. In this paper, we will use the data from a bearing which shows anomalous behavior in the later parts of test. In particular, out of the three data sets, we use set one where an inner race fault occurred on Bearing $3$. In the analysis, we use data from Bearing $3$.
\section{Markov Modeling}

{\bf e}gin{figure}
	\centering 

	\includegraphics[width=0.75\textwidth]{figures/autocorr.eps}
	\caption{Autocorrelation function of time-series data during the unstable phase of combustion. The time-series data is down-sampled by the lag marked in red square. It is noted that the individual time-series have their own down-sampling lags.}
	\label{fig:autocorr}
\end{figure}

{\bf e}gin{figure*} 
	\centering
	\subfloat[Probability density function for the pressure time series data]{\includegraphics[width=0.75\textwidth]{figures/fig1.eps}\label{fig:datadistribution}}\quad
	\subfloat[Spectral decomposition of the stochastic matrix for 1-step Markov model]{\includegraphics[width=0.75\textwidth]{figures/fig3.eps}\label{fig:spectralprop}}\\
	\caption{The first plate in the above Figure shows the change in the empirical density calculated for the pressure time-series data as the process deviates from the stable operating condition to unstable operating condition. The second plate shows the spectral decomposition of the 1-step stochastic matrix for the data under stable and unstable operating conditions.}
	\label{fig:databehavior}

\end{figure*}

In this section, we present results for modeling and analysis of the time-series data which are presented in this paper.
\subsection{Combustion}
Time-series data is first normalized by subtracting the mean and dividing by the
standard deviation of its elements; this step corresponds to bias removal and
variance normalization. Data from engineering systems is typically oversampled
to ensure that the underlying dynamics can be captured (in the current experiments, it was $\SI{8192}{\hertz}$ ). Due to coarse-graining from the symbolization process, an over-sampled
time-series may mask the true nature of the system dynamics in the symbolic
domain (e.g., occurrence of self loops and irrelevant spurious transitions in
the Markov chain). Time-series is first down-sampled to find the next crucial
observation. The first minimum of auto-correlation function generated from the
observed time-series is obtained to find the uncorrelated samples in time. The
data sets are then down-sampled by this lag. The autocorrelation function for the time-series data during unstable case is shown in Figure~\ref{fig:autocorr} where the data are downsampled by the lag marked in \mbox{col}or{red} red \mbox{col}or{black} rectangle in Figure~\ref{fig:autocorr}. To avoid discarding significant amount of data due to downsampling, down-sampled data using different initial
conditions is concatenated. Further details of this preprocessing can be found
in~\cite{Srivastav2014}.

The continuous time-series data set is then partitioned using maximum entropy
partitioning (MEP), where the information rich regions of the data
set are partitioned finer and those with sparse information are partitioned
coarser. In essence, each cell in the partitioned data set contains
(approximately) equal number of data points under MEP. A ternary alphabet with
$\alphaphabetSet=\{0,1,2\}$ has been used to symbolize the continuous combustion
instability data. As discussed in section~\ref{sec:experiment}, we analyze data sets from different phases, as the process goes from stable through the transient to the unstable region (the ground truth is decided using the RMS-values of pressure).

{\bf e}gin{figure} 
	\centering
	\subfloat[Hierarchical cluster tree of stable states]{\includegraphics[width=0.75\textwidth]{figures/fig4.eps}\label{fig:clusterstable}}\quad
	\subfloat[Hierarchical cluster tree of unstable states]{\includegraphics[width=0.75\textwidth]{figures/fig5.eps}\label{fig:culsterunstable}}\\
	\caption{State clustering under stable and unstable conditions.}
	\label{fig:clusterbehavior}

\end{figure}

In Figure~\ref{fig:datadistribution}, we show the observed changes in the behavior of the data as the combustion operating condition changes from stable to unstable. A change in the empirical distribution of data from unimodal to bi-modal is observed as the system moves from stable to unstable. We selected $150$ samples of pressure data from the stable and unstable phases each to analyze and compare. First, we compare the expected size of temporal memory during the two stages of operation. There are changes in the eigenvalue decomposition rate for the 1-step stochastic matrix calculated from the data during the stable and unstable behavior, irrespective of the combustor length and inlet velocity. During stable conditions, the eigenvalues very quickly go to zero as compared to the unstable operating condition (see Figure~\ref{fig:spectralprop}). This suggests that the size of temporal memory of the discretized data increases as we move to the unstable operating condition. This indicates that under the stable operating condition, the discretized data behaves as symbolic noise as the predictive power of Markov models remain unaffected even if we increase the order of the Markov model. On the other hand, the predictive power of the Markov models can be increased by increasing the order of the Markov model during unstable operating condition, indicating more deterministic behavior. An $\epsilonsilon=0.05$ is chosen to estimate the depth of the Markov models for both the stable and unstable phases. Correspondingly, the depth was calculated as $2$ and $3$ for the stable and unstable conditions (see Figure~\ref{fig:databehavior}). The corresponding $D(\epsilonsilon)$ is used to construct the Markov models next. First a PFSA whose states are words over $\alphaphabetSet$ of length $D(\epsilonsilon)$ is created and the corresponding maximum-likely parameters ($\bm{M}$ and $\bm{\Pi}$) are estimated. Then, the hierarchical clustering algorithm using K-L distance is used to cluster and aggregate the states. It is noted that we create individual models for every sample, i.e., every sample is partitioned individually so that the symbols will have different meaning (i.e., they represent different regions in the measurement space of the signals) for every sample. Consequently, each sample will have a different state-space when viewed in the continuous domain. Thus, we do not show the mean behavior of the samples during any operating regime as the state-space would be inconsistent (even though the cardinality could be the same).

In Figure~\ref{fig:clusterbehavior}, we show the hierarchical cluster tree which details the structure of the state-space for the PFSA with depth $D(\epsilonsilon)$ for a typical sample during stable and unstable behavior. The cluster tree also suggests the symbolic noise behavior of the data during the stable regime (the states are very close to each other based on the K-L distance). However, clearly a coarse clustering of states in the model during the unstable behavior would lead to significant information loss (as the states are statistically different). However, to compare the two Markov models, we keep the cardinality of the final models the same. For example, the algorithm is terminated with $3$ states in the final Markov model during the stable as well as the unstable regime. and the final aggregated states are the three clusters depicted in the Figure~\ref{fig:clusterbehavior}. Once the final aggregated states are obtained, we estimate the parameters of the model using the Bayesian inference discussed in section~\ref{subsec:DBN}.

Next, we present some results for model selection using the information-theoretic criteria discussed earlier in section~\ref{subsec:MDL}. BIC and AIC are used to select the model which achieves the minimum score. As seen in the Figures~\ref{fig:MDLstable} through~\ref{fig:MDLunstable}, the model with $5$ states is selected for stable as well as for the unstable case (note that the original model for the stable class had $9$ states for depth $2$ and the unstable model had $27$ states for a depth of $3$). In contrast to cross-validation, the two criteria provide an unsupervised way for model selection. Thus we see that we need much smaller state-space to preserve the temporal statistics of the data and AIC and BIC provide us with a technique to select the compact model.

{\bf e}gin{figure} 
	\centering
	\subfloat[Model scores using the BIC and AIC criterion during a typical stable condition]{\includegraphics[width=0.75\textwidth]{figures/fig11a.eps}\label{fig:MDLstable}}\quad
	\subfloat[Model scores using the BIC and AIC criterion during a typical unstable condition]{\includegraphics[width=0.75\textwidth]{figures/fig11b.eps}\label{fig:MDLunstable}}\\
	\caption{Unsupervised model selection under stable and unstable conditions.}
	\label{fig:MDLscores}

\end{figure}
In Figure~\ref{fig:HammingCombustion}, we show the Hamming distance between the sequences generated by the original model and the reduced models for a typical sample each from stable and unstable combustion. The box-plots are generated by simulating the original model and the reduced-order model to generate symbol sequences of length $1000$ from $100$ different initial states (i.e., a total of $100$ strings are generated) and the Hamming distance between them is calculated. A bound on the Hamming distance between the sequences generated by the original model and final model is also calculated using the inequality~\eqref{eqn:bound}. The results are shown in Figure~\ref{fig:HammingCombustion}. It is possible to use the proposed Hamming distance metric to select a final model; however, this measures the distance between the distributions induced by the Markov models, and model selection using it is left as a future work. It is noted that the bounds on Hamming distance can provide a computationally convenient way to select model scores as it can be found from the symbol emission probabilities of the model instead of explicitly looking at the predictive capability by looking at the likelihoods of the symbol sequences.

{\bf e}gin{figure} 
	\centering
	\subfloat[Hamming distance between the original and final models for a typical stable combustion sample]{\includegraphics[width=0.75\textwidth]{figures/HammingDistanceStable.eps}\label{fig:Hammingstable}}\quad
	\subfloat[Hamming distance between the original and final models for a typical unstable combustion sample]{\includegraphics[width=0.75\textwidth]{figures/HammingDistanceUnStable.eps}\label{fig:Hammingunstable}}\\
	\caption{Box plot for Hamming distance between the original and reduced-order models obtained after merging based on the results in Section~\ref{sec:analysis}}
	\label{fig:HammingCombustion}

\end{figure}
\subsection{Bearing}
The same procedure of downsampling and depth estimation is followed for analysis of bearing data as was described in the previous section for combustion. A ternary alphabet is again chosen to discretize the continuous data after downsampling and the maximum entropy partitioning is used to find the partitions. Using the spectral method, a depth of $2$ (i.e., a total of $9$ states) is estimated for an $\epsilonsilon=0.02$ (we skip the plot of spectral decomposition plot for brevity). The BIC and AIC score for the different models is shown in Figure~\ref{fig:ModelBearing} and the model with five states is selected using the obtained scores (marked in black rectangle). In Figure~\ref{fig:HammingDistBearing}, we show the Hamming distance between the sequences generated by the original model (with $9$ states) and the reduced models and the corresponding bounds obtained by inequality~\eqref{eqn:bound}.

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/MDL_Bearing.eps}
	\caption{Model scores using the BIC and AIC criteria; selected models are depicted by black rectangles.}
	\label{fig:ModelBearing}
\end{figure}

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/HDDistribution.eps}
	\caption{Box plot of the Hamming distance between the original and reduced-order models along with the analytical bound presented in Section~\ref{sec:analysis}.}
	\label{fig:HammingDistBearing}
\end{figure}

\section{Classification and Anomaly Detection Results}
In this section, we present some results for anomaly detection and classification using the pressure time-series data to infer the underlying reduced-order Markov model. As we discussed earlier in section~\ref{subsec:combustionexperiment}, the exact transition point of the system from stable to unstable is unknown, we first present results on anomaly detection and clustering of the data into different clusters which can be then associated with the stable and unstable class. We will present two different metrics for anomaly detection that allows models of different state-space and structure to be compared. It is noted that the word metric is used here in a loose sense; it is meant to be a distance that could be used to compare two different Markov models.

\subsection{Anomaly Detection}
As individual time-series have different state-space, we define some metrics to compare them. These metrics reflect changes in the information complexity of Markov models and reveal different behavior of combustion process based on the changes in the inferred data model. In particular, the following two metrics are defined.
{\bf e}gin{enumerate}
\item Cluster Divergence: This measure is defined for individual Markov models based on the cluster structure of the state-space of the model. Physically, it represents the maximum statistical difference between the states of the Markov model measures using K-L distance. It is calculated for a particular model $\bm{M}hcal{M}$ as follows
{\bf e}gin{equation}\label{eqn:metric}
\Delta_{\bm{M}hcal{M}}=\max\limits_{q_i,q_j \in qSet} d(q_i,q_j)
\end{equation}
where $d$ is defined by equation~\eqref{eq:kldistance}.
\item Discrepancy Statistics: We measure the discrepancy between the i.i.d. statistics and the Markov statistics for the discretized data. This could be also interpreted as the information gain for Markov models. This measure also represents the information complexity of the data. If the i.i.d. statistics and the Markov statistics are very close, then the data has no temporal statistics; however, an increase in this measure would indicate the information gain by creating a temporal Markov model for the data. This is measured by the following equation.
{\bf e}gin{equation}
H_\bm{M}hcal{M}=\sum_{q \in qSet} {\partial r}ob(q) D_{KL}({\partial r}ob(\alphaphabetSet\mid q)\| {\partial r}ob(\alphaphabetSet))
\end{equation}
where ${\partial r}ob(\alphaphabetSet\mid q)$ represents the symbol emission probability conditioned on a state $q$ of the Markov model and ${\partial r}ob(\alphaphabetSet)$ represents the marginal symbol emission probability. The term $D_{KL}$ represents the symmetric K-L distance between the two distributions.
\end{enumerate}
In Figure~\ref{fig:divergence}, we present some results to show the behavior of $\Delta_{\bm{M}hcal{M}}$ with increasing pressure fluctuations. It is noted that every model has been created in an unsupervised fashion by first discretizing and then, estimating the memory of the discrete sequence. As seen in Figure~\ref{fig:complexorig}, there are three distinct behaviors that can be associated with $\Delta_\bm{M}hcal{M}$. With low pressure fluctuations, the metric is very close to $0$, indicating that the states of the model are very similar statistically. This is seen until data number $200$ with corresponding $P_{rms}\sim 0.065 $ psig, which leads to a gradual change to a point where the measure saturates with $P_{rms}\sim 0.12$ psig (when the process becomes unstable). Thus, with this gradual trend with increasing pressure fluctuations, we associate different behaviors with the process. However, as is seen in the Figure~\ref{fig:complexorig}, the transition from stable to unstable behavior is not clearly defined and is very difficult to label during the experiments as the process is very fast. We show the pressure signals from the three different clusters in Figure~\ref{fig:pressure} where it could be seen that the sample number $250$ could be seen to approach an approximate limit cyclic behavior (and thus, could be loosely classified as transient stage).  An important point to note at this point is that this measure is independent of any operating conditions and only depends on stability (or instability) of the process. This metric is thus used for anomaly detection. In Figure~\ref{fig:complexfinal}, we show the statistics of $\Delta_{\bm{M}hcal M}$ with four states. We see that there is some loss of information up on merging states in the unstable class; the stable cluster remains unchanged implying that the states are statistically similar and the model distortion up on merging of states is insignificant. Thus, $\Delta_{\bm{M}hcal M}$ can be reliably used to detect departure from stable behavior.

The statistics for the discrepancy measure for the full state models is shown in Figure~\ref{fig:InfoGain}. The plot in Figure~\ref{fig:InfoGain} also agrees qualitatively with the earlier results on $\Delta_{\bm{M}hcal M}$. From these plots, we can infer that the Markov statistics for the stable cluster is very similar to the i.i.d. statistics and thus the data is very much independently distributed and conditioning on the inferred states of the Markov models doesn't improve predictability (or information complexity) of the temporal model. Thus, these two measures help infer the changes in the behavior of the data during the combustion process and are useful for anomaly detection.

{\bf e}gin{figure} 
	\centering
	\subfloat[$\Delta_{\bm{M}hcal{M}}$ for the full state model for the time-series data with increasing pressure root mean square]{\includegraphics[width=0.75\textwidth]{figures/fig9a.eps}\label{fig:complexorig}}\\
	\subfloat[Typical pressure signals from the three clusters seen in Figure~\ref{fig:complexorig}]{\includegraphics[width=0.75\textwidth]{figures/PressureSignal_fig9b.eps}\label{fig:pressure}}
\end{figure}
{\bf e}gin{figure}\ContinuedFloat 
  	\centering
	\subfloat[$\Delta_{\bm{M}hcal{M}}$ for models with $4$ states for the time-series data with increasing pressure root mean square]{\includegraphics[width=0.75\textwidth]{figures/fig9c.eps}\label{fig:complexfinal}}

	\caption{Anomalous behavior of data in the combustion process}
	\label{fig:divergence}

\end{figure}

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/fig10.eps}
	\caption{Variation of discrepancy statistics $H_{\bm{M}hcal M}$ with increasing pressure fluctuations. This also shows an anomaly around the point $200$ and qualitatively agrees to the behavior of $\Delta_{\bm{M}hcal M}$.}
	\label{fig:InfoGain}
\end{figure}

To see more explicitly the changes in the underlying models, the models during stable and unstable phases are visualized in the information space. To do this, we reduce the state space of the models to just $2$ states and estimate the corresponding emission parameters. As the models have three symbols, the emission matrix has $2$ rows and each row corresponds to the symbol emission probabilities conditioned on the two states. Each of these rows for $100$ cases from stable and $100$ cases from unstable are plotted on a single simplex plane which is shown in Figure~\ref{fig:Simplex}. The Figure shows the clusters of stable and unstable cases in the information space and that the model with even $2$ states are clustered separately. This shows that there is a structured change in the temporal dynamics of the data at the two phases and that the inferred Markov models are able to capture this change. Furthermore, the distinctive features of the models are sufficiently retained even after significant reduction in the state-space of the models.

{\bf e}gin{figure} 
	\centering
	\includegraphics[width=0.75\textwidth]{figures/simplex_2states_1.eps}
	\caption{Cluster of stable and unstable phase in information space. Each point is a row of the emission matrix for the reduced Markov model with $2$ states. The plot shows the change in the Markov model as the process moves from stable and unstable. \mbox{col}or{red} Red diamonds \mbox{col}or{black} represent the unstable phase while \mbox{col}or{green} green diamonds \mbox{col}or{black} represent the stable phase.}
	\label{fig:Simplex}
\end{figure}
\subsection{Classification}
These models are then used to train classifiers using support vector machines (SVM) and decision trees (DT)~\cite{B06}. The rationale behind using multiple classifier is to show that the performance of the Markov models is independent of the classification technique (i.e., it works equally well with maximum margin classifiers or decision tree classifiers). The SVM classifier is trained using a radial basis function kernel while the decision tree is trained using the standard euclidean distance. The classifiers are trained with $100$ data points from each class and are tested on the remaining data (around $80$ and $380$ for stable and unstable respectively). The tests are repeated for $100$ different train and test data sets from the total data. The results of classification accuracy are listed in Table~\ref{table:classification}. The SVM classifier is able to achieve around $1.67\%$ error using models with $2$ states while the decision tree classifier is able to achieve around $4.70\%$ error using models with $4$ states. This provides another way of selecting the final model for state merging in a supervised learning setting. It is noted that the original models contain $9$ states for stable and $27$ for unstable class.
{\bf e}gin{table}
  \caption{Performance of classifiers with different number of states. Mean Error= Lower is better. }\label{table:classification}
  \centering
  {\bf e}gin{tabular}{cccc}
    \hline
    Number of States & Classifier & Classification Error ($\%$)\\
    \hline
    \multirow{2}{*}{$9$} & SVM & $3.48 \pm 0.74$ \\
    & DT & $9.83\pm 3.24$ \\
    \hline
    \multirow{2}{*}{$8$} & SVM & $3.62 \pm 0.71$ \\
    & DT & $9.38 \pm 3.11$ & \\
    \hline
    \multirow{2}{*}{$7$} & SVM & $2.87 \pm 0.68$ \\
    & DT & $7.70 \pm 2.61$ & \\
    \hline
    \multirow{2}{*}{$6$} & SVM & $2.48 \pm 0.61$ \\
    & DT & $7.00 \pm 2.55$ & \\
    \hline
    \multirow{2}{*}{$5$} & SVM & $2.05 \pm 0.54$ \\
    & DT & $6.10 \pm 2.17$ & \\
    \hline
    \multirow{2}{*}{$4$} & SVM & $1.86 \pm 0.43$ \\
    & DT & $4.72 \pm 2.29$ & \\
    \hline
    \multirow{2}{*}{$3$} & SVM & $1.69 \pm 0.45$ \\
    & DT & $5.56 \pm 1.90$ & \\
    \hline
    \multirow{2}{*}{$2$} & SVM & $1.67 \pm 0.43$ \\
    & DT & $4.83 \pm 1.80$ & \\
    \hline
  \end{tabular}
\end{table}

\section{Summary, Conclusions and Future Work}
In recent times the idea of representation learning has become very popular in the machine learning literature as it allows decoupling of data for model learning  from the end-objectives like classification or clustering. In this paper, we presented a technique for Markov modeling of time-series data using concepts of symbolic dynamics which allows inference of model structure as well as parameters for compact data representation. In the proposed technique we first estimate the memory size of the discretized time-series data. The size of memory is estimated using spectral decomposition properties of the one-step Markov model created from the symbol sequence. Then, a second pass of data is made to infer the model with the right memory and the corresponding symbol emission matrix is estimated. Then the equivalence class of states based on K-L distance between the states are estimated using hierarchical clustering of the corresponding states of the Markov model. The proposed concepts were validated using two different datasets-- combustion instability and bearing. Modeling of combustion instability still remains a puzzle in the combustion community. The Markov modeling technique was used to analyze the problem of combustion instability. The proposed ideas were tested on experimental data from a swirl-stabilized combustor used to study unstable thermo-acoustic phenomenon during combustion process. The proposed approach allows us to infer the complexity of the time-series data based on the inferred Markov model. Two different metrics were proposed for anomaly detection and classification of the stable and unstable classes. The results presented in this paper are encouraging as the inferred models are able to identify the stable and unstable phases independent of any other operating condition.

Simultaneous optimization of discretization and memory estimation for model inference is a topic of future research. While the results obtained with Markov modeling for the combustion instability problem are inspiring, further investigation with transient data is required for better characterization of the process. More thorough comparison of the proposed models with HMM models of similar state-space size is also an important topic of future work.
 
\section*{Acknowledgments}
The authors would like to thank Professor Domenic Santavicca and Mr. Jihang Li of Center for Propulsion, Penn State for kindly providing the experimental data for combustion used in this work.

\end{document}