Factorial Hidden Markov Models (fHMM)

• fHMM are sequential latent-variable models where multiple independent hidden chains jointly determine observations through a coupled emission model.
• They offer greater expressive capacity than single-chain HMMs by representing states as a Cartesian product, enabling robust modeling in domains like energy disaggregation and audio processing.
• Advanced inference methods such as tensorized filtering and variational techniques address scalability challenges inherent in the exponential state space of fHMM.

A factorial hidden Markov model (fHMM) is a sequential latent-variable model in which multiple independent hidden Markov chains evolve in parallel and jointly determine the observation distribution at each time step. In contrast to a standard HMM, which uses a single latent chain, an fHMM represents the hidden state as a Cartesian product of chain-specific states, preserving conditional independence in the dynamics while coupling the chains through the emission model. This structure yields a distributed latent representation with substantially greater expressive capacity than a single-chain HMM, and it has been used for additive Gaussian, multinomial, non-negative dictionary, copula-based, and physics-informed emissions in domains including energy disaggregation, audio source separation, genomics, natural language, traffic modeling, systems biology, and environmental event classification (Barrios et al., 8 Jul 2026, Li et al., 2015).

1. Core formulation

For $M$ latent chains, with chain $m$ taking values $z_t^{(m)} \in \{1,\ldots,K_m\}$, the canonical factorial transition structure is

$$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$

With independent initial distributions, the joint model can be written as

$$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$

This is the defining distinction of the fHMM: the prior dynamics factorize across chains, but the emission couples them at each time slice (Barrios et al., 8 Jul 2026).

The same model can be viewed as a dynamic Bayesian network with $M$ hidden nodes per time slice and one observed node, with edges $z_{t-1}^{(m)} \to z_t^{(m)}$ for each chain and $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$ for the emission. In that sense, the fHMM is not a separate probabilistic formalism so much as a particular structured HMM family whose graphical semantics are especially convenient when several latent causes jointly explain one observation sequence (Wang, 2012).

A standard construction flattens the joint latent configuration into a single HMM state. If chain $m$ has $K_m$ states, the equivalent single-chain HMM has

$m$0

states, with flattened initial and transition operators expressed as Kronecker products of the chain-specific quantities. This equivalence is exact, but it obscures the factorial structure that motivates the model in the first place (Barrios et al., 8 Jul 2026).

Emission models vary by application. In additive Gaussian fHMMs, widely used in the classical formulation and in later work on model selection, observations satisfy

$m$1

or, in appliance-level energy disaggregation,

$m$2

with one chain per appliance and chain states corresponding to characteristic power levels (Li et al., 2015, Zhong et al., 2014). In discrete-observation LLMs, emissions can instead be log-linear softmax functions over the vocabulary conditioned on all chains, which turns the latent configuration into a distributed representation for each token (Nepal et al., 2013).

2. Inference and computational structure

The central computational difficulty of the fHMM is that exact filtering, smoothing, and decoding over the flattened state space scale with $m$3, not with the individual $m$4. Standard forward filtering after flattening has complexity

$m$5

and storing the full flattened transition matrix requires a dense $m$6 object. This is the basic source of the state-space blow-up that defines the computational literature on fHMMs (Barrios et al., 8 Jul 2026).

Several inference families have been developed to mitigate this. Structured mean-field approximations decouple the posterior across chains while preserving per-chain temporal dependence; in non-negative audio models this yields inference complexity linear in the number of sources, rather than exponential, while remaining close to exact inference in SDR and much faster in practice (Mysore et al., 2012). Stochastic variational inference without message passing replaces explicit forward–backward recursions with amortized copula-based variational distributions and subchain sampling, thereby removing per-iteration dependence on the full sequence length and enabling training on sequences with $m$7 on the order of $m$8 (Ng et al., 2016). Factored expectation propagation has also been used in input–output fHMMs, where free-form structured variational inference is combined with expectation constraints and EP site updates for transition and emission parameters (1305.4153).

Exactness can also be preserved while reorganizing the computation. Tensorized filtering keeps messages, transitions, and emissions in their natural multidimensional form rather than flattening them. In that formulation, forward messages are tensors $m$9, the recursion is implemented by successive mode-wise contractions with the per-chain transition matrices, and the cost becomes

$z_t^{(m)} \in \{1,\ldots,K_m\}$0

rather than quadratic in $z_t^{(m)} \in \{1,\ldots,K_m\}$1. In the equal-size case $z_t^{(m)} \in \{1,\ldots,K_m\}$2, this is $z_t^{(m)} \in \{1,\ldots,K_m\}$3 instead of $z_t^{(m)} \in \{1,\ldots,K_m\}$4. On $z_t^{(m)} \in \{1,\ldots,K_m\}$5 randomized trials, the tensorized and flattened filters were numerically equivalent up to floating-point error; the reported speedups were approximately $z_t^{(m)} \in \{1,\ldots,K_m\}$6 when varying $z_t^{(m)} \in \{1,\ldots,K_m\}$7 and $z_t^{(m)} \in \{1,\ldots,K_m\}$8 when varying $z_t^{(m)} \in \{1,\ldots,K_m\}$9, and the tensorized implementation scaled to $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$0 on the test hardware while the flattened implementation failed around $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$1 due to memory exhaustion (Barrios et al., 8 Jul 2026).

A different line of work exploits sparsity in the emission factor graph rather than tensor structure. The Graph Filter and Graph Smoother localize the Bayes update by discarding distant likelihood factors according to graph distance, and their error is controlled in a local total variation norm. The resulting bounds are described as “dimension-free” in the sense that they do not necessarily deteriorate with overall model dimension, provided the local graph constants remain bounded (Rimella et al., 2019).

For multimodal posteriors, especially in Bayesian settings, posterior exploration rather than forward recursion can become the main bottleneck. Augmented ensemble MCMC addresses this by combining parallel tempering with an auxiliary-variable crossover move inspired by genetic algorithms, producing always-accepted exchange moves between tempered chains and materially improving mode jumping in simulation and cancer-genomics applications (Märtens et al., 2017).

3. Learning, identifiability, and model selection

Parameter estimation in fHMMs spans supervised maximum likelihood, EM-type latent-variable learning, variational Bayes, and fully Bayesian sampling. In supervised energy disaggregation, the model parameters are estimated directly from disaggregated appliance-level training data: empirical state means for power levels, empirical frequencies for initial states, empirical transition probabilities for homogeneous or time-stratified transitions, and residual variance for emission noise (Zhong et al., 2014). In Bayesian traffic and systems-biology variants, latent states, transition parameters, and regime-specific emission parameters are updated jointly by MCMC or expectation-propagation-based variational procedures (Zhang et al., 17 Jun 2025, 1305.4153).

A central theoretical issue is identifiability. For factorial Gaussian models, including FHMMs with Gaussian emissions, the emission matrix is not identifiable even if the true assignment matrix is known. The core algebraic reason is that the combinatorial assignment matrix $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$2 has rank

$$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$3

which implies a $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$4 dimensional null space and hence nontrivial perturbations of the dictionary that leave the likelihood unchanged. The paper giving this result shows that identifiability can be restored under a one component sharing assumption, where each factor shares a common emission component $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$5, reducing the number of free columns to $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$6 and matching the rank bound (Subakan et al., 2015).

That same work reframes emission estimation as dictionary learning. Under the shared-component factorial model and an incoherence condition on the shared atom, the emission matrix can be recovered from cluster centers $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$7 without alternating between the estimation of $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$8 and $$P(\mathbf{z}_t \mid \mathbf{z}_{t-1}) = \prod_{m=1}^M P\big(z_t^{(m)} \mid z_{t-1}^{(m)}\big), \qquad \mathbf{z}_t=(z_t^{(1)},\dots,z_t^{(M)}).$$9. This is significant because it isolates an identifiability pathology that is intrinsic to the factorial parameterization, rather than to a particular optimization algorithm (Subakan et al., 2015).

Model selection is another nontrivial aspect of fHMM learning because one must choose both the number of chains and the number of states per chain. Factorized Asymptotic Bayesian inference for FHMMs addresses this by integrating out transition probabilities while applying a Laplace approximation only to emission parameters. The resulting shrinkage factors penalize low-occupancy states and are shown to eliminate redundant states almost surely under the paper’s similarity assumptions. Empirically, refined FAB outperformed iFHMM and variational FHMM baselines in model-selection accuracy while maintaining competitive held-out perplexity (Li et al., 2015).

4. Major structured variants

A substantial part of the fHMM literature consists of structured variants that preserve the factorial prior while adding domain-specific constraints to the transitions or emissions.

In energy disaggregation, the factorial non-homogeneous HMM (FNHMM) introduces time-varying transition matrices $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$0 so that appliance usage can depend on time of day. The interleaved factorial non-homogeneous HMM (IFNHMM) adds a second constraint: at most one chain is allowed to change state per time step,

$$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$1

On 100 households from the Household Electricity Survey with 2-minute sampling, mean reconstruction errors $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$2 were reported as $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$3 for FHMM, $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$4 for FNHMM, $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$5 for IFHMM, and $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$6 for IFNHMM, with large household-to-household variability; the coefficient of variation of IFNHMM error was $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$7 (Zhong et al., 2014).

In single-channel speech separation, gain-adapted FHMM (GFHMM) augments a two-chain speech FHMM with a hidden gain parameter $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$8 representing target-to-interference ratio. The emission model uses a MIXMAX approximation in the log-spectrum, while $$p(\mathbf{z}_{1:T}, y_{1:T}) = \prod_{m=1}^M p(z_1^{(m)}) \prod_{t=2}^T p(z_t^{(m)} \mid z_{t-1}^{(m)}) \prod_{t=1}^T p(y_t \mid \mathbf{z}_t).$$9 is estimated by iterative quadratic optimization combined with factorial Viterbi decoding. Experiments on 180 mixtures with gain differences from $M$0 to $M$1 dB showed that GFHMM significantly outperformed both standard FHMM and its memoryless vector-quantization counterpart, while converging in fewer than three iterations on average (Radfar et al., 2019).

In audio source separation with non-negative spectrograms, the non-negative FHMM replaces Gaussian emissions with a PLCA/LDA-style discrete emission process over dictionary atoms and frequency quanta. A Bayesian version places a Dirichlet prior on framewise mixture weights $M$2, conditioned on the current HMM states through a mask over available atoms. Its variational inference algorithm was reported to achieve around a $M$3 increase in speed relative to exact inference in typical configurations, while maintaining comparable SDR (Mysore et al., 2012).

Input–output FHMMs incorporate exogenous covariates directly into the transition model. In systems biology, binary transcription-factor chains are driven by metabolite inputs through sigmoid or continuous-time transition functions, while gene-expression outputs are modeled by a linear-Gaussian emission with a sparse TF–gene matrix $M$4. The inference procedure combines EP for $M$5 and $M$6 with path mean-field updates for the latent state trajectories (1305.4153).

A different kind of structural decomposition appears in FHMM-IDM for car-following behavior. There, two latent chains are used: one for intrinsic driver action and one for external traffic scenario. The acceleration emission is Gaussian around the Intelligent Driver Model mean, with regime-specific IDM parameters chosen by one chain and a Gaussian scenario model over kinematic covariates chosen by the other. On the HighD dataset, the model uncovered interpretable regimes such as “Cautious Following,” “Aggressive Following,” “Congested Cruising,” “Steady-State Following,” and “High-Speed Seeking,” together with five traffic scenarios (Zhang et al., 17 Jun 2025).

5. Applications across domains

The breadth of fHMM applications follows directly from the model’s semantics: each chain can represent one latent source, subsystem, regime, clone, or contextual factor, while the emission accounts for their joint observational consequences.

Domain	Factorial interpretation	Representative outcome
Energy disaggregation	One chain per appliance	IFNHMM reduced mean reconstruction error from $M$7 to $M$8 on 100 households (Zhong et al., 2014)
Single-channel speech separation	One chain per speaker, plus gain adaptation	GFHMM significantly outperformed FHMM on 180 mixtures with gain differences from $M$9 to $z_{t-1}^{(m)} \to z_t^{(m)}$0 dB (Radfar et al., 2019)
Audio source separation	One chain per source, with state-gated dictionary atoms	Bayesian N-FHMM achieved around a $z_{t-1}^{(m)} \to z_t^{(m)}$1 speed increase in typical configurations (Mysore et al., 2012)
Tumor heterogeneity	One chain per clone, plus a normal-cell chain	HetFHMM outperformed PyClone, PhyloSub, and Rec-BTP on simulated data (Haffari et al., 2015)
Natural-language representation learning	Multiple latent chains per token	FHMM posterior features reduced PoS error from $z_{t-1}^{(m)} \to z_t^{(m)}$2 to $z_{t-1}^{(m)} \to z_t^{(m)}$3 in the reported setup (Nepal et al., 2013)
Traffic behavior modeling	Separate chains for driver action and traffic scenario	FHMM-IDM uncovered interpretable regime–scenario co-occurrence patterns (Zhang et al., 17 Jun 2025)
Haze/dust event classification	Separate chains for haze and dust presence	Micro-F1 $z_{t-1}^{(m)} \to z_t^{(m)}$4; Dust F1 improved from $z_{t-1}^{(m)} \to z_t^{(m)}$5 to $z_{t-1}^{(m)} \to z_t^{(m)}$6, Haze F1 from $z_{t-1}^{(m)} \to z_t^{(m)}$7 to $z_{t-1}^{(m)} \to z_t^{(m)}$8 (Zhang et al., 21 Aug 2025)

In tumor heterogeneity, the factorial interpretation is especially literal: each chain represents a clone-specific genotype trajectory along the genome, and the observed read counts arise from the mixture of clone-specific genotypes weighted by cellular prevalence. HetFHMM uses $z_{t-1}^{(m)} \to z_t^{(m)}$9 chains, with chain $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$0 representing normal cells and the remaining chains representing tumor clones, and each genotype state belongs to a 21-state copy-number/LOH space. On simulated data, it achieved better clustering quality and lower $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$1 for cellular prevalence than the reported baselines (Haffari et al., 2015).

In natural language processing, the same mathematical structure serves a different purpose. Instead of interpreting the chains as physical sources, the model treats them as distributed latent features for words in sequence. Posterior marginals over the chains become global-context-sensitive token representations, and with $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$2 layers of $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$3 states each the reported part-of-speech error fell to $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$4 from a $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$5 baseline, while noun-phrase chunking F1 reached $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$6 from a baseline $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$7 (Nepal et al., 2013).

Environmental classification offers a further example of emission flexibility. In the haze–dust model, the factorial prior keeps the Haze and Dust processes statistically independent in their transitions, while the emission uses a Gaussian copula to capture nonlinear dependence among PM10, wind speed, visibility, and relative humidity. Mutual-information weighting is then injected into Viterbi decoding to emphasize informative variables, which materially improved rare-class performance under strong class imbalance (Zhang et al., 21 Aug 2025).

6. Assumptions, limitations, and open directions

Across the literature, several assumptions recur. The latent chains are typically independent in their transitions; the observational coupling is confined to the emission model; dwell times are geometric unless explicit duration modeling is added; and the emission family is application-specific but usually chosen for tractability, such as additive Gaussian, multinomial, log-normal, or copula-based forms (Barrios et al., 8 Jul 2026, Zhong et al., 2014). These assumptions are often reasonable first approximations, but they are also the main source of model mismatch.

The most persistent limitation is computational. Even when flattening is avoided, storage still scales with $\{z_t^{(1)},\ldots,z_t^{(M)}\} \to y_t$8 in tensorized exact filtering, and approximate methods often require either restrictive factorization assumptions or sophisticated engineering. This suggests that scalability is not a single problem but a collection of related ones: time complexity, memory complexity, posterior multimodality, and optimization instability are all distinct failure modes (Barrios et al., 8 Jul 2026, Ng et al., 2016, Märtens et al., 2017).

Identifiability is a second structural limitation. The nonidentifiability of the emission matrix in standard factorial Gaussian models is not an optimization artifact; it is an intrinsic symmetry of the model class. Shared-component constructions, shrinkage, informative priors, or downstream supervision can mitigate this, but they do not remove the general lesson that richer latent factorizations can introduce new ambiguities as well as new expressiveness (Subakan et al., 2015, Li et al., 2015).

Domain-specific variants add further caveats. The interleaving constraint in IFNHMM can be violated when multiple appliances change state within the same low-frequency observation window. Gain-adapted speech models assume quasi-stationary gains over a mega-frame. FHMM-IDM inherits the regime-label switching issues of mixture models and HMMs. Graph-localized approximations depend on bounded graph-degree and correlation-decay conditions. These are not incidental implementation details; they are substantive modeling assumptions that delimit the operating regime of each variant (Zhong et al., 2014, Radfar et al., 2019, Zhang et al., 17 Jun 2025, Rimella et al., 2019).

Current research directions are correspondingly diverse. Tensorized backward recursions, full EM schemes, and exploitation of sparsity or low rank are explicit next steps in scalable exact inference (Barrios et al., 8 Jul 2026). In disaggregation and traffic modeling, richer covariates and stronger duration models such as semi-Markov extensions are natural continuations (Zhong et al., 2014, Zhang et al., 17 Jun 2025). In factorial Gaussian models, extending identifiability theory beyond one shared component remains open (Subakan et al., 2015). A plausible implication is that the most successful future fHMMs will not be generic replacements for HMMs, but carefully structured models that use factorization only where it aligns with substantive latent mechanisms and where the computational advantages can actually be realized.