Higher-Order Markov Chains: Theory & Applications

• Higher-order Markov chains are stochastic processes whose transitions depend on a finite history of states, extending the memory of classical models.
• They utilize high-order transition tensors and algebraic geometric techniques to model complex dependencies and reduce computational challenges.
• Efficient computation methods, including tensor power iterations and skeleton reduction, enable accurate estimation and dimensionality reduction in practical applications.

A higher-order Markov chain is a discrete- or continuous-time stochastic process in which the probability distribution of the next state depends on a finite, fixed number $m-1$ of previous states, generalizing the memoryless (first-order) Markov property to arbitrary depth. The transition dynamics are encoded by a high-order transition tensor, and such chains arise naturally in time series modeling, natural language processing, computational biology, and the analysis of systems with extended memory. Research on arXiv has contributed rigorous geometric, algebraic, statistical, and computational frameworks for understanding, classifying, and estimating higher-order Markov processes and their variants.

1. Formal Definition, Transition Tensors, and State Space

Let $\mathcal{S} = \{1,2,\dots,n\}$ be a finite state space. A time-homogeneous Markov chain of order $m-1$ is a stochastic process $\{X_t\}$ such that

$$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$

where $\mathcal{P} = [p_{i_1 i_2 \ldots i_m}]$ is an order-$m$ column-stochastic tensor: $$p_{i_1 i_2 \ldots i_m} \geq 0, \qquad \sum_{i_1=1}^{n} p_{i_1 i_2 \ldots i_m} = 1 \quad \forall (i_2,\ldots,i_m)\in \mathcal{S}^{m-1}.$$ When $m=2$, this reduces to the transition matrix for a first-order Markov chain. Every $m$-th order chain on $\mathcal{S} = \{1,2,\dots,n\}$0 can be embedded as a first-order Markov chain on $\mathcal{S} = \{1,2,\dots,n\}$1, but this lifted chain has a highly sparse and structured transition kernel (Xu, 3 Oct 2025).

The table below summarizes key mathematical objects:

Order	Transition Object	Dimension
1	Matrix $\mathcal{S} = \{1,2,\dots,n\}$2	$\mathcal{S} = \{1,2,\dots,n\}$3
$\mathcal{S} = \{1,2,\dots,n\}$4	Tensor $\mathcal{S} = \{1,2,\dots,n\}$5	$\mathcal{S} = \{1,2,\dots,n\}$6 ($\mathcal{S} = \{1,2,\dots,n\}$7 modes)

2. Geometry and Algebraic Structure

Higher-order Markov models, particularly in the Mixture Transition Distribution (MTD) class, exhibit rich algebraic geometry. For the $\mathcal{S} = \{1,2,\dots,n\}$8-th order MTD, transitions are parameterized as mixtures of first-order processes: $\mathcal{S} = \{1,2,\dots,n\}$9 with $m-1$0 and stochastic $m-1$1 (Sturmfels, 2012). The Zariski closure of the model in projective space forms an explicit algebraic variety. In the binary case ($m-1$2), the model is a linear space sliced into simplices, while for $m-1$3, it forms a cone over a Segre variety, cut out by a Gröbner basis of linear forms and $m-1$4 minors.

This algebraic structure has the following implications:

• Parameters are generically algebraically identifiable.
• All defining relations among joint probabilities are generated in degrees $m-1$5.
• The model geometry directly explains the presence of multiple local maxima in the likelihood for $m-1$6 and higher (Sturmfels, 2012).

3. Stationary Distributions, Regularity, and Long-Term Behavior

For a higher-order chain, a probability vector $m-1$7 is stationary if, for all $m-1$8,

$m-1$9

where $\{X_t\}$0 is the $\{X_t\}$1-step transition tensor defined by

$\{X_t\}$2

(Han et al., 10 Jun 2025). A key result is that if the transition tensor is regular (i.e., $\{X_t\}$3 for some $\{X_t\}$4), then there is a unique stationary distribution $\{X_t\}$5 attracting all initial distributions (Han et al., 10 Jun 2025).

Notably, the stationary set can have surprisingly rich structure—there exist $\{X_t\}$6-th order chains whose stationary vector set is any prescribed face of the simplex, or even the entire simplex, with explicit characterizations for second and higher orders (Li et al., 2013). Computational methods such as tensor power iteration and momentum-accelerated variants provide practical means for finding stationary distributions, with convergence guarantees under contraction conditions and regularity assumptions (Yu et al., 2020).

4. Classification, Skeleton Reduction, and Recurrence

Structural properties such as recurrence, periods, and irreducibility are determined by the transition kernel's "skeleton," a minimal directed graph derived from the tensor (Gallesco et al., 10 Jan 2026). The skeleton captures the essential context length and transition constraints, providing:

• Direct classification of recurrent classes and their periods via corresponding classes in a reduced binary matrix.
• Algorithmic reduction in complexity by focusing on the true minimal memory (skeleton order $\{X_t\}$7), often yielding dramatic computational savings.
• Efficient criteria for essential irreducibility, based on the existence of unique closed classes in the skeleton.

For example, a $\{X_t\}$8th-order chain may have skeleton order $\{X_t\}$9, reducing a $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$0-dimensional problem to $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$1 (Gallesco et al., 10 Jan 2026).

5. Estimation, Dimensionality Reduction, and Parsimonious Modeling

The curse of dimensionality—i.e., $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$2 parameters—has driven the development of parsimonious higher-order models:

• MTD Models: Parameterize higher-order chains as mixtures of first-order transitions—yielding linear rather than exponential parameter growth and enhancing interpretability (Sturmfels, 2012, Ost et al., 2022, Heiner et al., 2019).
• RHOMP: Retrospective chain models depend on a single past state per transition, yielding parameter count $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$3, efficient convex estimation, and robust predictive performance in user trajectory modeling (Wu et al., 2017).
• Max Markov Chain (MMC): Assumes a sparse, winner-take-all mechanism among lags, leading to closed-form MLEs after permutation search, and tractable inference via greedy or hill-climbing procedures (Zhang et al., 2022).
• Bayesian Nonparametrics: Conditional tensor factorization (CTF) places shrinkage or Dirichlet process priors on the tensor, enabling data-driven lag selection and hierarchical modeling of dependencies (Sarkar et al., 2015, Heiner et al., 2019).

The following table summarizes parameter counts for full and reduced models:

Model	Parameter Count
Full order-$$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$4	$$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$5
Mixture Transition Dist.	$$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$6
RHOMP/MMC	$$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$7
Bayesian CTF (sparse)	variable, << $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$8

6. Computational Approaches and Software

Efficient algorithms for computations with higher-order chains exploit the tensor structure:

• Tensor Box Product: Generalizes matrix multiplication for $$\Pr(X_{t+1}=i_1 \mid X_t=i_2,\,X_{t-1}=i_3,\,\ldots,\,X_{t-m+2}=i_m) = p_{i_1 i_2 \ldots i_m}$$9-way transitions; underlies $\mathcal{P} = [p_{i_1 i_2 \ldots i_m}]$0-step tensor powers and recursion for transition probabilities (Xu, 3 Oct 2025).
• Power Methods: Higher-order power iteration, momentum-based acceleration, and nonparametric quadratic extrapolation yield convergence to stationary distributions under weak contraction hypotheses (Yu et al., 2020).
• Reduction to First-Order: Algorithms form the augmented first-order chain ($\mathcal{P} = [p_{i_1 i_2 \ldots i_m}]$1 matrix) using the tensor's matricization and Khatri–Rao products (Xu, 3 Oct 2025).
• Skeleton Extraction: Prunes context trees and computes minimal order for efficient recurrence/irreducibility classification (Gallesco et al., 10 Jan 2026).
• Software: HOMC (MATLAB package) implements these operations, provides diagnostic routines, and enables experimentation with higher-order chains (Xu, 3 Oct 2025).

7. Extensions: Continuous-Time, Higher-Order SDEs, and Applications

Continuous-time generalizations embed the memory of a process in a function-valued state, leading to stochastic differential equations (SDEs) on path space. The transition law depends on the recent history via functionals, and parameter estimation is tackled via projected maximum likelihood and differential-equation solvers (Nag, 2021). Applications include high-memory stochastic volatility models and analogs of classical SDEs with non-Markovian drift/diffusion terms.

In combinatorial random structures, higher-order chains are essential for modeling walks with memory on hypergraphs, where the next move depends on group structure and the last $\mathcal{P} = [p_{i_1 i_2 \ldots i_m}]$2 visited nodes. Tensor-based frameworks succinctly encode transition probabilities and provide convergence and steady-state results via higher-order Perron–Frobenius theory (Cui et al., 8 Apr 2026). Such constructions illuminate phenomena unobservable by classical first-order Markovian or graph-based approaches.

---

The theory and methodology around higher-order Markov chains, as codified in the recent arXiv literature, incorporate algebraic geometry, tensor methods, Bayesian shrinkage, structural graph reduction, and dynamical systems. This body of work unifies classical probability, computational algebra, and applied statistics, enabling the scalable, interpretable, and theoretically grounded analysis of processes with extended memory.