Fitting Sparse Markov Models to Categorical Time Series Using Regularization

Let ${\mathbb{N}}=\{1,2,\ldots\}$ be the set of all positive integers. Let ${\cal X}_{n}=(X_{1},\ldots,X_{n})$ and $\tilde{X}_{t}^{(m)}=(X_{t},X_{t-1},\ldots,X_{t-m+1})$, for $m\geq 1$, $t\in{\mathbb{N}}$. Write $w$ for an ordered (finite) sequence of $\Sigma$-elements of length $|w|$. Let $wu$ denote the (ordered) concatenation of $w$ and $u$. Write $|\Sigma|=d$ and $\Sigma^{m}=\{\sigma_{1},\ldots,\sigma_{p}\}$ so that $p=|\Sigma|^{m}$. Let $N_{w}=\sum_{t=|w|}^{n}{1\!\!1}(\tilde{X}_{t}^{(|w|)}=w)$ where ${1\!\!1}(\cdot)$ denotes the indicator function. For any $S\subset\Sigma^{m}$ and $a\in\Sigma$, define $N_{S}=\sum_{t=m}^{n-1}{1\!\!1}(\tilde{X}_{t}^{(m)}\in S)$, $N_{S,a}=\sum_{t=m}^{n-1}{1\!\!1}(\tilde{X}_{t}^{(m)}\in S,X_{t+1}=a)$. In particular, $N_{\sigma_{j}}$ denotes the number of times the chain $\tilde{X}_{t}^{(m)}$ hits the $m$-tuple $\sigma_{j}$, and $N_{\sigma_{j},a}$ is the number of transitions from $\sigma_{j}$ to $a$.

Next we define the probabilities associated with the SMM (1.2). For $j=1,\ldots,p$ and $a\in\Sigma$, let

and $\bm{\pi}_{j}$ be the corresponding transition probability vector. Note that by the SMM property, for any $a\in\Sigma$, the transition probability $\bm{\pi}_{j,a}$ is a constant over all $j$ such that $\sigma_{j}\in{\cal C}_{i}$. However, we do not know the sets ${\cal C}_{i}$ and one of the challenges of fitting an SMM to a dataset is to be able to idenitify the sets ${\cal C}_{1},\ldots,{\cal C}_{k_{0}}$. To that end, define nonparametric estimators of $\bm{\pi}_{j,a}$ using their empirical versions:

where $N=n-m+1$ denotes the total number of $m$th order histrory in the observed variables ${\cal X}_{n}$.

Here we propose a new approach to fitting the SMM based on regularization.

Consider the penalized criterion function

over $\mathbf{b}_{j}=(b_{ja}:a\in\Sigma)\in\Pi_{d}$ for $j=1,\ldots,p$, where $\lambda>0$ is a penalty parameter, $\Pi_{d}$ is the $d$-dimensional simplex $\Pi_{d}=\{(u_{1},\ldots,u_{d})\in[0,1]^{d}:u_{1}+\ldots+u_{d}=1\}$ and where $\rho(\cdot,\cdot)$ is a distance measure between two $d$-dimensional probability vectors. Thus, (2.3) treats the estimators $\hat{\pi}_{ja}$ as (correlated) “observations” and penalizes the distance between all distinct pairs of the probability vectors in order to identify the identical probability vectors. In particular, the number of parameters grow approximately quadratically in the number of observations and with a suitable choice of the penalization term, one can identify the identical probability vectors. When $\rho(\mathbf{b}_{j_{1}},\mathbf{b}_{j_{2}})^{2}=\sum_{\ell=1}^{d}(b_{j_{1}\ell}-b_{j_{2}\ell})^{2}$, (2.3) gives a version of the Group LASSO of Yuan and Lin (2006)[15] that is designed for selecting pairs of full vectors that are close, and we have a convex optimization problem that can be solved for large $p$s. On the other hand, if we use the $\ell_{1}$ distance $\rho(\mathbf{b}_{j_{1}},\mathbf{b}_{j_{2}})=\sum_{\ell=1}^{d}|b_{j_{1}\ell}-b_{j_{2}\ell}|$, then only componentwise zero differences can be identified.

Once we minimize the criterion function in (2.3), it is a relatively easy task to find estimates of $k_{0}$ and the sets ${\cal C}_{i}$. Specifically, we start with a pair with the smallest $j_{1}$ and seek all $j_{2}>j_{1}$ that the distance between the solutions $b^{*}_{j_{1}}$ and $b^{*}_{j_{2}}$ is zero. Then, we set $\hat{{\cal C}}_{1}$ to be the set consisting of $j_{1}$ and all such $j_{2}$. In the next step, we consider all pairs that are not in $\hat{{\cal C}}_{1}$ and repeat the procedure until all pairs with estimated zero distances have been grouped. In case there are indices $j$ for which none of the estimated paired distances are zero, we keep them as a singletons, that is groups consisting of single elements. This gives the estimated groups $\hat{{\cal C}}_{i}:i=1,\ldots,\hat{k}$, with $\hat{k}$ giving an estimate of $k_{0}$.

The traditional clustering methodologies like $K$-means have many limitations. In most cases, we have to pre-specify the number of clusters, along with the possibility that we end up with a local minima instead of the global one. The advantage of clustering by solving the equation (2.3) for a range of $\lambda$ is that we get a solution path from at most $p$ many singleton clusters to only one cluster consisting of all the elements. Subsequently, we can fix some criterion function which will enable us to find the optimum cluster assignment among the all possible models in the solution path. Hence, not only we won’t need to fix the number of clusters beforehand, also we avoid the problems of being stuck at the local minima. Several efficient algorithms have been developed in recent years to solve the equation (2.3) when the penalty function $\rho$ is convex; e.g $\rho(b_{j1},b_{j2})=\lVert b_{j_{1}}-b_{j_{2}}\rVert_{p}$ for some $p\geq 1$. Pelckmans et al. (2005)[9], Lindsten et al. (2011)[7], Hocking et al. (2011)[5] and others recently proposed this method and established this to be more robust and scalable in comparison to the traditional ones. Lindsten et al. (2011)[7] used an off-the-shelf convex solver CVX to solve the convex clustering problem, which suffers from scalability issues. The theoretical perfect cluster recovery conditions have been derived by Zhu et al. (2014)[16] only for two clusters, while Panahi et al. (2017)[8] derived perfect recovery conditions for general $k$ clusters, but under the assumption of uniform weights. Sun et al. (2021)[13] provides sufficient conditions for theoretical recovery conditions under more general weight choices. They have also developed a faster algorithm called semismooth Newton based augmented Lagrangian method (SS-NAL), and derived the convergence criteria for their algorithm.

Chi and Lange (2015)[3] have developed an elegant method of solving the convex clustering problem by augmented lagrangian methods. For $\rho(x)=\lVert x\rVert_{2}$, we first view solving the equation (2.3) as the following constrained optimization problem

where $\mathcal{E}$ is the set of all distinct edges $\{l:l=(l_{1},l_{2}),l_{1}<l_{2},w_{l}>0\}$. Here, a new splitting variable $\mathbf{v}_{l}$ has been introduced to capture the difference between the group centroids, which helps the optimization procedure much easier. Chi and Lange (2015)[3] have used two algorithms for solving this constrained optimization problem, namely ADMM and AMA. For both these algorithms, one first have to incorporate augmented lagrangian as follows:

where $\mathbf{B},\mathbf{V}$ and $\bm{\Gamma}$ are the matrices with $\mathbf{b}_{j},\mathbf{v}_{l}$ and $\bm{\gamma}_{l}$ for $j=1,...,p$ and $l\in\mathcal{E}$ in their columns respectively. Splitting the variables in such fashion would allow us to update $\mathbf{B}$, $\mathbf{V}$ and $\bm{\Gamma}$ sequentially, given the other variables. The convergence of ADMM does not depend on the choice of $\nu$, it will converge for any $\nu>0$. On the other hand, AMA converges for any $0<\nu<2/p$. The performance of both these algorithms have been compared, and it is established that AMA is much faster than ADMM, especially when the weights are sparse. The major advantage of AMA lies in the inherent structure. After some basic algebra, Chi and Lange (2015)[3] have proved that we only need to update $\mathbf{B}$ and $\bm{\Gamma}$ in every step, and we can bypass updating $\mathbf{V}$ applying some linear relationship among the variables. Since, AMA provides much faster result, we use this method in our analysis. Suppose, $\mathbf{B}^{(t)}$ and $\bm{\Gamma}^{(t)}$ be the parameter values in the $t^{th}$ step. The updates in the next step are computed using the following relations:

where $\mathbf{g}_{l}^{(t+1)}=\mathbf{b}_{l_{1}}^{(t+1)}-\mathbf{b}_{l_{2}}^{(t+1)}$, $C_{l}=\{\bm{\gamma}_{l}:\lVert\bm{\gamma}_{l}\rVert_{2}\leq\lambda w_{l}\}$, and $\mathcal{P}_{A}(\mathbf{x})$ is the projection of $\mathbf{x}$ onto the set $A$. We continue till convergence. The convergence criterion has been discussed in Chi and Lange (2015)[3] in details, using the dual problem and duality gap.

So far, we have discussed the numerical methods to solve (2.3) for a given $\lambda$. But it is important to choose an optimum value of $\lambda$ for the optimization problem. In this section, we propose a data driven method to select this tuning parameter using BIC criterion. For a given $\lambda$, denote the obtained clusters as $\hat{\mathcal{C}}_{1}(\lambda),\ldots,\hat{\mathcal{C}}_{k_{\lambda}}(\lambda)$, where $k_{\lambda}$ is the number of clusters. Define the common transition probability for the $m$-tuples in the estimated group $\hat{\mathcal{C}}_{\alpha}(\lambda)$ as

The log-likelihood of the observations under the obtained cluster assignment for a particular $\lambda$ is given by

Hence, the BIC score corresponding to the obtained model is

By a grid search over a range of possible $\lambda$ values, we select the $\lambda$ for which BIC is minimized. The solution of the equation (2.3) corresponding to that $\lambda$ is considered as the estimated cluster assignment. In the next section, we provide theoretical results which will demonstrate that for a range of $\lambda$ values, we will be able to perfectly recover the true clusters for large $n$.

Consider solving the equation (2.3) with $\rho(b_{j_{1}},b_{j_{2}})=\lVert b_{j_{1}}-b_{j_{2}}\rVert_{2}$, and the optimum solution is denoted by $b_{i}^{*}(\lambda)$, for $i=1,2,...,p$. Let the true partition of the state space $\Sigma^{m}$ is $\{{\cal C}_{1},\ldots,{\cal C}_{k_{0}}\}$; with the corresponding transition probability vectors being $R_{1},...,R_{k_{0}}$ such that $R_{\ell,a}=P(X_{t+1}=a\big{|}Y_{t}=\sigma_{\ell})$. Denote $p_{\ell}=\big{|}{\cal C}_{\ell}\big{|}$, the size of the $\alpha^{th}$ partition. Following the notations of Yuan et al. (2021), define

Suppose, the following conditions hold.
(A1)  $w_{j_{1},j_{2}}=w_{j_{2},j_{1}}$ and $w_{j_{1},j_{2}}>0$ for any $j_{1},j_{2}\in{\cal C}_{\ell}$, $\ell=1,2,...,k_{0}$.
(A2)  $p_{\alpha}w_{i,j}>\mu_{i,j}^{(\alpha)}$, $\forall i,j\in{\cal C}_{\alpha}$ and $\forall\alpha=1,2,...,k_{0}$.

In other words, the condition (A1) implies the symmetry in the weight choices, and ideally the weight should be positive between two $m$-tuples belonging to the same partition. On the other hand, (A2) determines a lower bound for the weight between two $m$-tuple in a particular group. We will use these conditions to prove our results. Before going into the main results, let us recall a few other results as follows.
Proposition 1 (Consistency of Estimated Transition Probabilities):  Suppose $X_{1},\ldots,X_{n}$ be an SMM of order $m$, with the partitions $\{\mathcal{C}_{1},\ldots,\mathcal{C}_{k_{0}}\}$. Then, as $n\to\infty$,
(a) $\hat{\bm{\pi}}_{j}\xrightarrow[]{p}\bm{R}_{\alpha}$ for $j\in\mathcal{C}_{\alpha}$;
(b) $\dfrac{N_{\sigma_{j}}}{N}\xrightarrow[]{p}q_{j}$, where $q_{j}$ is the stationary probability of the state $\sigma_{j}$;
(c)

where $\Sigma_{\alpha}=diag(\bm{R}_{\alpha})-\bm{R}_{\alpha}\bm{R}_{\alpha}^{(T)}$. Since $\Sigma_{\alpha}$ is of rank $|\Sigma|-1$, the asymptotic Normal distribution is singular.
Proposition 2 (Sun et al. (2021)):  Suppose the above conditions hold, and $\lambda_{\text{min}}^{(n)}<\lambda_{\text{max}}^{(n)}$. Then for any $\lambda\in(\lambda_{\text{min}}^{(n)},\lambda_{\text{max}}^{(n)})$, $\bm{b}_{i}^{*}(\lambda)=\bm{b}_{j}^{*}(\lambda)$ for $i,j\in\mathcal{C}_{\alpha}$; $\alpha=1,..,k_{0}$ and $\bm{b}_{i}^{*}(\lambda)\neq\bm{b}_{j}^{*}(\lambda)$ for any $i\in\mathcal{C}_{\alpha},j\in\mathcal{C}_{\beta},\alpha\neq\beta$. In other words, for any $\lambda\in(\lambda_{\text{min}}^{(n)},\lambda_{\text{max}}^{(n)})$, we recover the true partition of the state space.
Thus, Proposition 1 tells us about the weak consistency and the CLT for the transition probability vectors. Proposition 2 deals with perfect recovery conditions under general weight choices, derived in Sun et al. (2021). These propositions will be the key tools in proving our results. In the next section, we provide our major theoretical findings.

The very first result will ensure that the resulting solution of the objective function in (2.4) will produce a valid probability distribution over $\Sigma$. Note that, we don’t consider this solution as our estimated transition probability for the group. However, if under some extra-ordinary circumstances the solution exhibits oracle properties, one can use the solutions as the common transition probabilities for the partitions.
Theorem 1:  For any $\lambda>0$, the optimal solution $\bm{b}_{i}^{*}(\lambda)$ is a valid probability distribution for $i=1,2,,,,p$; i.e.
(a) $b^{*}_{i,a}(\lambda)\geq 0$ for $a=1,...,d$,
(b) $\sum_{a=1}^{d}b^{*}_{i,a}(\lambda)=1$.
Next, we would like to derive the probability of true cluster recovery. There are two steps involved in this process. First, we need the true model in the solution path over varying $\lambda$. This implies the conditions of the proposition 2 must be satisfied, i.e $\lambda_{\text{min}}^{(n)}<\lambda_{\text{max}}^{(n)}$. From Theorem 2, we get a lower bound of the probability of the true model being present in the solution path. This will tell us that with increasing sequence length, the probability of not perfect recovery converges to $0$ in an exponential rate.
Theorem 2:  Define

Then, under the above conditions (A1) and (A2), as $n\to\infty$,

for $0<\epsilon<\dfrac{\delta\delta_{1}}{\delta_{1}+\delta_{2}}$, and for some constants $C_{1}^{(\alpha)},C_{2,j}>0$.
Once we have the true model in the solution path, the next step is to establish that the probability of selecting that model through BIC criterion converges to $1$ as $n\to\infty$. The next theorem will prove this statement.
Theorem 3:  Suppose the conditions of the Theorem (2) holds, and $\lambda_{\text{min}}^{(n)}<\lambda_{\text{max}}^{(n)}$. For any $\lambda$, denote the clustering assignment obtained by minimizing the equation (2.3) as $M_{\lambda}=\{\hat{{\cal C}}_{1}(\lambda),\ldots,\hat{{\cal C}}_{k_{\lambda}}(\lambda)\}$; where $k_{\lambda}$ is the number of clusters. Suppose $\ell_{n}(\lambda)$ be the log-likelihood of the observations corresponding to the cluster assignment $M_{\lambda}$, and the corresponding BIC score is $BIC_{n}(\lambda)=-2\ell_{n}(\lambda)+k_{\lambda}(d-1)\log n$. Choose some $\lambda_{0}\in(\lambda_{\text{min}}^{(n)},\lambda_{\text{max}}^{(n)})$. Then, for any $\lambda$ such that $M_{\lambda}\neq M_{\lambda_{0}}$,

as $n\to\infty$.
Next, we will provide some sufficient conditions for perfect cluster recovery under a particular weight choice using Gaussian kernels. These weights have been used in the previous analysis of convex clustering, and proven to produce good clustering results.
Theorem 4:  Define $p_{min}=\min_{\alpha}p_{\alpha},p_{max}=\max_{\alpha}p_{\alpha}$. Suppose the following assumptions hold.
(a) $w_{i,j}=e^{-\phi\lVert\hat{\bm{\pi}}_{i}-\hat{\bm{\pi}}_{j}\rVert_{2}^{2}}l^{k}_{i,j}$, where $l^{k}_{i,j}$ is the indicator function that $\hat{\bm{\pi}}_{i}$ belongs to $k$ nearest neighbour of $\hat{\bm{\pi}}_{j}$ or vice versa, for some $\phi>0$.
(b) $k\geq p_{max}-1$.
(c) For some $\epsilon<\epsilon_{max}=\dfrac{\delta}{2}-\dfrac{1}{2\phi\delta}\log\Big{(}2\Big{(}\dfrac{k+1}{p_{min}}-1\Big{)}\Big{)}$, $\lVert\hat{\bm{\pi}}_{j}-\bm{R}_{\alpha}\rVert_{2}<\dfrac{\epsilon}{2};$ $\forall j\in\mathcal{C}_{\alpha}$, $\forall\alpha=1,...,k_{0}$.
Under the above assumptions, the conditions (A1) and (A2) are satisfied. Moreover, $\delta_{1}\geq p_{min}e^{-\phi\epsilon_{max}^{2}}-2(k+1-p_{min})e^{-\phi(\delta-\epsilon_{max})^{2}}=\delta^{(min)}_{1}$, $\delta_{2}\leq 2(k+1-p_{min})e^{-\phi(\delta-\epsilon_{max})^{2}}=\delta_{2}^{(max)}.$

Here, we take $|\Sigma|=4$, the usual scenario when analyzing the DNA sequences. The order of the chain $m$ is taken to be both $2$ or $3$. For $m=2$, we equally divide the all $16$ tuples into $4$ groups of $4$ elements; and for $m=3$, we divide the $64$ triplets into $8$ groups of equal size $8$. For a particular group $g$, we generate $Z_{g,\ell}$ independently from $Unif(0,1)$, for $\ell=1,2,3,4$. The transition probability of that group is generated from Dirichlet distribution with parameter $(e^{Z_{g,1}},...,e^{Z_{g,4}})$. As of weights, we first take $w_{i,j}=1$ for all $i,j=1,2...,p,i<j$. Next we choose some sparse weights depending on the distance between the estimated transition probabilities $\hat{\bm{\pi}}_{i}$ and $\hat{\bm{\pi}}_{j}$. We have used the $k$-nearest neighbour based weights as proposed in Chi and Lange (2015), such as $w_{i,j}=\exp^{-\phi\lVert\hat{\bm{\pi}}_{i}-\hat{\bm{\pi}}_{j}\rVert_{2}^{2}}l^{k}_{i,j}$ where $l^{k}_{i,j}$ is the indicator function that $\hat{\bm{\pi}}_{i}$ belongs to $k$ nearest neighbour of $\hat{\bm{\pi}}_{j}$ or vice versa, for some $\phi>0$. In this example we use $\phi=100$. We also incorporate a third choice of weight, namely $w_{i,j}=\exp^{-\phi\lVert\hat{\bm{\pi}}_{i}-\hat{\bm{\pi}}_{j}\rVert_{\infty}^{2}}l^{k(\infty)}_{i,j}$ where $l^{k(\infty)}_{i,j}$ is the similar indicator function, but the nearest neighbour is computed w.r.t $l_{\infty}$ distance. We use two different values of the nearest neighbour, $k=5$ and $k=3$. The results are provided in the following tables.

From this study, it is clear that for both $m=2$ and $m=3$, choice of uniform weight results in really poor performance in terms of recovering the true cluster. Although the $ARI$ increases with increasing $n$, we might need really large sample size to get good result. On the other hand, choosing the weights using the number of nearest neighbour $k=3$ perform much better than that with $k=5$ in terms of both $ARI$ and perfect recovery, especially for lower sample size. Note that the model is balanced in this example, and the optimum choice of $k$ is $3$ (by Corollary 1). Hence, we can justify that fact using our simulation study. On the other hand, the optimum choice of $k$ is $7$ for $m=3$, but the choice $k=5$ is reliable as well. $k=3$ makes the weights too sparse in this case, which results in degradation of the clustering accuracy a little bit. For both $m=2$ and $m=3$, the probability of true recovery increases with increasing $n$.

This experiment ensures pretty good result for high $n$, mostly for $n\geq 10000$. It is thus worthy of investigation under what circumstances we will be able to get very good recovery for lower $n$, such as $n=1000$ or $n=2000$. The theoretical results suggest that if the cluster centroids are well separated, clustering performance gets better even for lower sample size. In this experiment, we have $4$ groups for $m=2$, with the minimum centroid difference $?$ in terms of $l_{2}$ distance and $?$ in terms of $l_{\infty}$ distance; while these values are $?$ and $?$ for $m=3$. In the next simulation study, we will demonstrate how the clustering accuracy improves for well separated centroids.

In this study as well, we take $|\Sigma|=4$ and $m=3$. We divide this $64$ triplets into four groups of size $18,18,15$ and $13$. For the $\alpha^{th}$ group, $R_{\alpha,\alpha}=0.7,R_{\alpha,\beta}=0.1$, $\alpha=1,2,3,4$, $\beta=1,2,3,4$, $\alpha\neq\beta$. As the choice of weight, we have first used the $k=15$ nearest neighbour weights w.r.t the $l_{2}$ distance in the Gaussian kernel, and used $\phi=100$. The second choice of weight is $w_{i,j}=\exp^{-\phi\lVert\hat{\bm{\pi}}_{i}-\hat{\bm{\pi}}_{j}\rVert_{\infty}}l^{k(\infty)}_{i,j}$, with $\phi=10$. Note that, here we have used the $l_{\infty}$ distance to find the nearest neighbour, but instead of incorporating the Gaussian kernel, we have used the natural exponential decay. The third weight is similar to the second one, only the distance is $l_{1}$ instead of $l_{\infty}$. Here we have only used relatively smaller samples, $n=1000$ and $n=2000$ respectively. The results are portrayed in the following Table.

From the experiment, we can infer the weight involving $l_{2}$ distance in the Gaussian kernel performs poorly compared to $l_{\infty}$ or $l_{1}$ distance. Especially, usage of $l_{\infty}$ distance is really effective in such scenarios, as it measures the maximum element-wise distance between two estimated transition probability vectors. In this way, we are able to separate out two vectors which are not likely to be in the same cluster. From this study, we can fairly conclude that when the centroids are well-separated in such fashion so that their difference is significantly high in a small number of co-ordinates, use of $l_{\infty}$ distance can provide the best possible result.