A Hidden Markov Model Variant for Sequence Classification

1 Introduction

The need to operate on sequence data is prevalent in a variety of real world applications ranging from protein/DNA classification, speech recognition, intrusion detection and text classification. Sequence data can be distinguished from the more-typical vector representation in that the length of sequences within a dataset can vary and that the order of symbols within a sequence carries meaning.

For sequence classification, a variety of strategies, depending on the problem type, can be used to map sequences to a representation that can be handled by traditional classifiers. A simple technique involves selecting a fixed number of elements from the sequence and then using those elements as a fixed-length vector in the classification engine. In another technique, a small subsequence length, ℓ, is selected, and a sizeM ℓ vector is constructed containing the counts of all length ℓ subsequences from the original sequence. This vector can then be used for classification [Leslie et al., 2002]. A third method for classifying sequence data requires only a positive definite mapping defined over pairs of sequences rather than any direct mapping of sequences to vectors. This strategy, known as the kernel trick, is often used in conjunction with support vector machines (SVMs) and allows for a wide variety of sequence similarity measurements to be employed.

∗ Funding: NSF III-0905117.

Hidden Markov Models (HMM) [Rabiner and Juang, 1986; Eddy, 1998] have a rich history in sequence data modeling (in speech recognition and bioinformatics applications) for the purposes of classification, segmentation, and clustering. HMMs' success is based on the convenience of their simplifying assumptions. The space of probable sequences is constrained by assuming only pairwise dependencies over hidden states. Pairwise dependencies also allow for a class of efficient inference algorithms whose critical steps build on the ForwardBackward algorithm [Rabiner and Juang, 1986].

We present an HMM variant over a set of sequences, with one transition matrix per sequence, as a novel alternative for handling sequence data. After training, the per-sequence transition matrices of the HMM variant are used as fixed-length vector representations for each associated sequence. The HMM variant is also similar to a number of topic models, and we describe it in the context of Latent Dirichlet Allocation [Blei et al., 2003]. We then describe three methods to infer the parameters of our HMM variant, explore connections between these methods, and provide rationale for the classification behavior of the parameters derived through each.

We perform a comprehensive set of experiments, evaluating the performance of our method in conjunction with support vector machines, to classify sequences of amino acids into structural classes (fold recognition and remote homology detection problem [Rangwala and Karypis, 2006]). The combination of these methods, their interpretations, and their connections to prior work constitutes a new twist on classic ways of understanding sequence data that we believe is valuable to anyone approaching a sequence classification task.

2 Problem Statement

Given a set of N sequences, we would like to find a set of fixed-length vectors, A1...N , that, when used as input to a function f(A), maximize the probability of reconstructing the original set of sequences. Under our scheme, f(A) is a Hidden Markov Model variant with one transition matrix, An, assigned to each sequence, and a single emissions matrix, B, and start probability vector, a, for the entire set of sequences. By maximizing the likelihood of the set of sequences under the HMM variant model, we will also find the set of transition matrices that best represent our set of sequences. We further postulate that this maximum likelihood representation will achieve good classification results if each sequence is later associated with a meaningful label.

2.1 Model Description

We define a Hidden Markov Model variant that represents a set of sequences. Each sequence is associated with a separate transition matrix, while the emission matrix and initial state transition vector are shared across all sequences. We use the value of each transition matrix as a fixed-length representation of the sequence. We define the parameters and notation for the model in Table 1.

The joint probability of the model is shown below:

(1) p(x, z|a,A,B) =

This differs from the standard hidden Markov model only in the addition of a transition matrix, An (highlighted in bold in Equation 1), for each sequence, where the index n indicates a sequence in the training set. Under the standard HMM, a single transition matrix, A, would be used for all sequences.

To regularize the model, we further augment the basic HMM by placing Dirichlet priors on a, each row of A, and each row of B. The prior parameters are the uniform Dirichlet parameters γ, α, and β for a, A, and B respectively. The probability of the model with priors is shown below, where the prior probabilities are the first three terms in the product below and take the form Dir(x; a, K) = Γ(Ka) Γ(a)K ∏i xa-1 i : (2) p(x, z, a,A,B|α, β, γ ) =

One potential difficulty that could be expected in classifying simple HMMs by transition matrix is that the probability of a sequence under an HMM does not change under a permutation of the hidden states. This problem is avoided when we force each sequence to share an emissions matrix, which locks the meaning of each transition matrix row to a particular emission distribution. If the emission matrix were not shared, then two HMMs with permuted hidden states could have transition matrices that with large Euclidean distances. For instance, the following HMMs have different transition matrices, but the probability of an observed sequence is the same under each:

However, a Euclidean distance between their two transition matrices, A1 and A2 is large.

3 Background

3.1 Mixtures of HMMs

Smyth introduces a mixture of HMMs in [Smyth, 1997] and presents an initialization technique that is similar to our model in that an individual HMM is learned for each sequence, but differs from our model in that the emission matrices are not shared between HMMs. In [Smyth, 1997], these initial N models are used to compute the set of all pairwise distances between sequences, defined as the symmetrized log likelihood of each element of the pair under the other's respective model. Clusters are then computed from this distance matrix, which are used to initialize a set of K < N HMMs where each sequence is associated with one of K labels. Smyth notes that while the log probability of a sequence under an HMM is an intuitive distance measure between sequences, it is not intuitive how the parameters of the model are meaningful in terms of defining a distance between sequences. In this research, we demonstrate experimentally that the transition matrix of our model is useful for sequence classification when combined with standard distance metrics and tools.

3.2 Topic Models

Simpler precursors of LDA [Blei et al., 2003] and pLSI [Hofmann, 1999], which represent an entire corpus of documents with a single topic distribution vector, are very similar to the basic Hidden Markov Model, which assigns a single transition matrix to the entire set of sequences that are being modeled. To extend the HMM to a pLSI analogue, all that is needed is to split the single transition matrix into a per-sequence transition matrix. To extend this model to an LDA analogue, we must go a step further and attach Dirichlet priors to the transition matrices, as in our model.

Inference of the LDA model (Figure 1a) on a corpus of documents learns a matrix of document-topic probabilities. A row of this matrix, sometimes described as a mixed-membership vector, can be viewed as a measurement of how a given document is composed from the set of topics. In our HMM variant (Figure 1b), a single transition matrix, An, can be thought of as the analogue to a document-topic matrix row and can be viewed as a measurement of how a sequence is composed of pairs of adjacent symbols. The LDA model also includes a "topic-word" matrix, which indicates the probability of a word given a topic assignment. This matrix has the same meaning as the emissions matrix, B, in the HMM variant.

The Fisher kernel [Jaakkola and Haussler, 1999] and the Probabilistic Product Kernel [Jebara et al., 2004] (PPK), are principled methods that allow probabilistic models to be incorporated into SVM kernels. The HMM variant is similar to these methods in that it uses latent information from a generative model as input to a discriminative classifier. It differs from these methods, however, both in which portions of the generative model that are incorporated into the discriminative classifier and in the assumptions about how differences in generating distributions comparisons between training examples.

4 Learning the model parameters

4.1 Baum-Welch

A well-known method for learning HMM model parameters is the Baum-Welch algorithm. The Baum-Welch algorithm is an expectation maximization algorithm for the standard HMM model, and the basic algorithm is easily modified to learn the multiple transition matrices of our variant. The parameter updates shown below converges to a maximum a posteriori (MAP) estimate of p(z, a, A, B|x, γ, α, β) [Rabiner and Juang, 1986]:

where f and b are the forward and backward recursions defined below:

The complexity of the Baum-Welch-like algorithm for our variant is identical to the complexity of Baum-Welch for the standard HMM. The update for Aij in the original HMM involves summing over ∑n Tn terms, while the update for a single Anij is a sum over Tn terms, making the total number of terms over all the An's in our variant, ∑n Tn, the same number as the original algorithm.

4.2 Gibbs Sampling

Two Gibbs sampling schemes are commonly used to infer Hidden Markov Model parameters [Scott, 2002]. Unlike the Baum-Welch algorithm which returns a MAP estimate of the parameters, these sampling schemes allow the expectation of the parameters to be computed over the posterior distribution p(z, a, A, B|x, γ, α, β).

In the Direct Gibbs sampler (DG), hidden states and parameters are initially chosen at random, then new hidden states are sampled using the current set of parameters:

In the Forward Backward sampler (FB), the initial settings and parameter updates are the same as the DG scheme, but the hidden states are sampled in order from Tn down to 1 using values from the forward recursion. Specifically, each hidden state znt is sampled given znt+1 = j from a multinomial with parameters

In both algorithms, after the hidden states are sampled, parameters are sampled from Dirichlet conditional distributions, shown for A below, where I(ω) = 1 if ω is true and 0 otherwise:

The FB sampler has been shown to mix more quickly than the DG sampler, especially in cases where adjacent hidden states are highly correlated [Scott, 2002]. We therefore use the FB sampler in our implementation.

4.3 Variational Algorithm

Another approach for inference of the HMM variant parameters is through variational techniques. We employ a mean field variational algorithm that follows a similar pattern as EM. When the variational update steps are run until convergence, Kullback-Leibler divergence between the variational distribution, q(z, a, A, B), and the model's conditional probability distribution, p(z, a,A,B|x, γ, α, β), is minimized. The transition matrices returned by the variational algorithm are the expectations of those matrices under the variational distribution. Thus, like the Gibbs sampling algorithm, the parameters returned by the variational algorithm approximate the expectations of the parameters under the conditional distribution.

Our mean field variational approximation is shown below:

with variational parameters hnti, which approximate each znti, and ˜αnij, ˜ βim, and ˜ γi, which can be thought of as Dirichlet parameters approximating α, β, and γ.

When we maximize the variational free energy with respect to the variational parameters, we obtain the following update equations, where Ψ(x) = d log Γ(x) dx :

Notice that the update for hnti depends only on the adjacent h's, hnt-1i and hnt+1i as well as the expectations of the transition probabilities from the adjacent h's and the expectation of the emission probabilities from the current hnti. This mean field algorithm can therefore be understood as an equivalent of the Direct Gibbs sampling method except that at subsequent time steps interactions occur between variational parameters rather than through the sampled values of z. A complete derivation of the variational algorithm is included on the authors' website 1 .

1 http://www.cs.gmu.edu/~mlbio/ijcai11

Table 2: AUC results from all of the multi-class SVM experiments are displayed. The best performing algorithm, the best performing setting of K, and the best combination of K and algorithm is marked in bold. The Gibbs-Sampling-derived representation most frequently returned the best AUC score on the majority of the datasets.

5 Experimental Setup

5.1 Protocol

To evaluate our fixed-length representation scheme, for each dataset (described in Section 5.2), we created three sets of fixed-length representations per trial over ten trials by running each of the three inference algorithms: (i) Baum-Welch, (ii) Gibbs Sampling, and (iii) the mean field variational algorithm, on the entire set of input data. We varied the number of hidden states, K, from 5 to 20 in increments of 5. This procedure created a total of 120 (3×10×4) fixed-length representations for each dataset.

The fixed-length vector data was then used as input to a support vector machine (SVM) classifier 2 . We used the SVM to either perform either multiway classification on the dataset under the Crammer-Singer [Crammer and Singer, 2002] construction or the one-versus-rest approach, where a binary classifier was trained for each of the classes.

We compare classification results from our model with results from the Spectrum(2) kernel for all experiments. The Spectrum(ℓ) kernel is a string kernel whose vector representation is the set of counts of substrings of observed symbols length ℓ in a given string [Leslie et al., 2002]. For the one-versus rest experiments, we compare our results to more biologically sensitive kernels for protein classification, described in Rangwala et. al [Rangwala and Karypis, 2005].

5.2 Protein Datasets

The Structural Classification of Proteins (SCOP) [Murzin et al., 1995] database categorizes proteins into a multilevel hierarchy that captures commonalities between protein structure at different levels of detail. To evaluate our representation, we ran sets of protein classification experiments on the three top levels of the SCOP taxonomy: class, fold, and superfamily. Our datasets, which were obtained from previous studies [Rangwala and Karypis, 2006; Kuang et al., 2004], were derived from either the SCOP 1.67 or the SCOP 1.53 versions and filtered at 25% and 40% pairwise sequence identities. A protein sequence dataset filtered at 25% identity will have no two sequences with more than 25% sequence identity.

2 We used SVM-light and SVM-struct for classification (http://www.cs.cornell.edu/People/tj/svm_light/) [Joachims, 1999].

We partitioned the data into a single test and training set for each category. At the class level, the original dataset was split randomly in to training and test sets. To eliminate high levels of similarity between sequences that could lead to trivially good classification results, we imposed constraints on the training/test set partitioning for classification in the fold and superfamily experiments. For the fold level classification problem, the training sets were partitioned so that no examples that shared the fold and superfamily labels were included in both the training and test sets. Similarly, for the superfamily level classification problem (referred to as the remote homology detection problem [Leslie et al., 2002; Rangwala and Karypis, 2005]), no examples that shared the superfamily and family levels were included in both the training and test sets.

5.3 EvaluationMetrics

We evaluated each classification experiment by computing the area under the ROC curve (AUC), a plot of the true positive rate against the false positive rate, constructed by adjusting the SVM's intercept parameter. We also computed the AUC50 value, which is a normalized computation of the area under the ROC curve until the first 50 false positives have been detected. We were worried about variance over different Baum-Welch runs due to convergence of the algorithm to different local optima. To mitigate this concern, we ran both the Baum-Welch algorithm and the other inference algorithms, for consistency, 10 separate times on each dataset. The results presented for each inference method are averages over individual results of the 10 trials across the different classes.

6 Results and Discussion

6.1 Protein Sequence Classification

Table 2 shows a comparison of results (average AUC scores) across the inference algorithms in three taxonomic categories (class, fold, and superfamily) using the multiclass SVM. Although the AUC scores are close for each algorithm, in most cases, the Gibbs sampling algorithm outperforms the other algorithms.

Table 3 shows a comparison of results over the inference algorithms but only for the one-versus-rest superfamily classification experiment on the SCOP 1.53 dataset. Similar to the multiclass experiments using the linear kernel, the Gibbs sampling algorithm outperforms the other inference methods in the one-versus-rest experiments. Although the values of the best performing algorithm's AUC and AUC50 scores do not significantly change from the linear to the Gaussian kernel, the vari- Gaussian Kernel ational algorithm shows a large improvement, ranging from 6% to 30%.

6.2 Analysis of inference algorithms

The differences in AUC values resulting from the different training algorithms (Tables 2 and 3) can be explained, at least in part, by a high level overview of how each algorithm operates. While the Baum-Welch algorithm returns MAP parameters of the model, both the Gibbs sampling method and the variational algorithm return expectations of the parameters under an approximate of the posterior distribution. The MAP solution from the Baum-Welch algorithm is likely to reach a local maximum of the posterior, while the other algorithms should tend to average over posterior parameters.

The Gibbs sampling algorithm and the variational algorithm each compute expectations of the parameters under an approximate posterior distribution, but each uses a different method to construct this approximation. The variational algorithm will be less likely to converge to a good approximation of the marginal distribution because the mean field variational approximation necessarily does away with the direct coupling between adjacent hidden states characteristic of the HMM.

6.3 Comparative Performance

Tables 4 and 5 show a comparison between the HMM variant and common classification methods for the multiclass and one-versus rest experiments respectively. The AUC and AUC50 scores indicate that our scheme produces a representation that is roughly equivalent in power to the Spectrum kernel for protein classification. In defense of the HMM variant, the size of the vector representation produced by the spectrum kernel is significantly larger than the typical representations produced by our HMM variant. The Mismatch(5,1) kernel, used for SCOP 1.53 superfamily classification (Table 5), is similar to the Spectrum(5) kernel but also counts substrings of length 5 that differ by one amino acid residue from those found in an observed sequence. The size of the vector representation associated with this kernel can be up to 205 . This value is large compared to the largest vector representation in our experiments, which is 400 for the HMM variant with 20 hidden states. Nearly all of these high-performing kernel methods, unlike the HMM variant, employ domain specific knowledge, such as carefully tuned position-specific scoring matrices, to aid classification. In contrast, the only parameter that needs to be adjusted in the HMM variant is the number of hidden states.

7 Conclusions and FutureWork

Our HMM variant is an extension of the standard HMM that assigns individual transition matrices to each sequence in a dataset but keeps a single emissions matrix for the entire dataset. We describe three inference algorithms, two of which, a Baum-Welch-like algorithm and a Gibbs sampling algorithm, are similar to standard methods used to infer HMM parameters. A third, the variational inference algorithm, is related to algorithms used for inference on topic models and more complex HMM extensions. We demonstrate, by comparing results on protein sequence classification using our method in conjunction with SVMs, that each of these algorithms infers transition matrices that capture useful characteristics of individual sequences. Because our model fits within a large existing body of work on generative models, we are especially interested in related models that perform classification directly.