High-Performance Tensor Formulation of the Viterbi Algorithm for Hidden Semi-Markov Models

# High-Performance Tensor Formulation of the Viterbi Algorithm for Hidden Semi-Markov Models
Source: [https://arxiv.org/html/2609.16500](https://arxiv.org/html/2609.16500)
Lorenzo PiarulliElia BelliAffiliation:Department of Computer Science Sapienza University of Rome belli\.2006305@studenti\.uniroma1\.itDaniele De SensiAffiliation:Department of Computer Science Sapienza University of Rome desensi@di\.uniroma1\.it

###### Abstract

Hidden Semi\-Markov Models \(HSMMs\) are fundamental probabilistic models widely adopted across diverse domains, from computational biology to finance and signal processing\. The Viterbi algorithm decodes the most likely state sequence given an HSMM and can be applied iteratively for ab initio model learning\. However, existing Viterbi implementations remain sequential, and GPU\-accelerated solutions are entirely absent, making HSMM decoding impractical for large\-scale workloads\. We present a tensor\-based formulation of the Viterbi algorithm for HSMMs, restructuring the inner loops into tensor operations that naturally map onto SIMD units and massively parallel architectures\. Building on this formulation, we provide optimized implementations spanning single\- and multi\-core CPUs, and, for the first time, GPU\. Experimental evaluation demonstrates speedups of up to 14×\\timeson a single core, over 200×\\timeswith multi\-core, and over 570×\\timeson GPU over the state\-of\-the\-art sequential baseline, establishing a new performance baseline for large\-scale HSMM decoding\.

###### Index Terms:

Hidden Markov Models, Tensors, GPU

## IIntroduction

High\-Performance Computing \(HPC\) architectures are evolving at an extraordinary rate\. Over the past decades, we have transitioned from CPU\-only computation to massive multicore processors, and then to GPUs\. Now, driven by the rise of artificial intelligence, entirely new accelerator architectures are emerging, including dataflow engines, systolic arrays, and SIMD\-centric designs, that promise unprecedented throughput for structured, regular computations\. Yet, while hardware evolves rapidly, the algorithms that run on it do not always keep up\. Application scientists tend to be conservative: their software frameworks are enormously complex, and restructuring a working codebase to exploit a new architecture is a daunting, error\-prone endeavor\. In many cases, the cost and complexity of migration simply outweigh the perceived benefit, and teams understandably choose to keep a working pipeline rather than risk breaking it for uncertain gains\. As a result, many fundamental algorithms, including those that today underpin astonishing scientific discoveries, remain anchored to decades\-old sequential formulations that hide parallelization possibilities and, consequently, remain confined to single\-threaded CPU execution\. Worse still, practitioners often resort to simplified or truncated versions of their models simply because the full, general formulation would be computationally infeasible on the sequential hardware\.

Hidden Semi\-Markov Models \(HSMMs\) are versatile probabilistic frameworks with applications spanning diverse fields, from genome annotation and segmentation in computational biology\[[1](https://arxiv.org/html/2609.16500#bib.bib3)\]to finance\[[2](https://arxiv.org/html/2609.16500#bib.bib4)\], speech recognition\[[3](https://arxiv.org/html/2609.16500#bib.bib5)\], and signal processing\[[4](https://arxiv.org/html/2609.16500#bib.bib26)\]\. As a generalization of the classicalHidden Markov Model \(HMM\), an HSMM describes a system transitioning through a finite set of hidden states at discrete time intervals\. In this paradigm, the system’s internal state remains unobservable, manifesting only through visible emissions governed by state\-specific probabilities, while state progression is regulated by a fixed transition matrix\.

A practical illustration of this framework is found in sleep\-cycle prediction from physiological data: heartbeat measurements serve as the observations, while the underlying sleep stages represent the hidden states to be inferred\. Similarly, in HPC applications like genome annotation, a DNA sequence is modeled as a series of nucleotides \(observations\); the objective is to classify each nucleotide as belonging to either a coding or non\-coding region \(hidden states\)\. Crucially, these systems often exhibit temporal persistence, remaining in a specific state for numerous consecutive time steps before transitioning, a characteristic that motivates the use of state\-duration modeling\.

Standard HMMs, however, are fundamentally limited in capturing these temporal dynamics, as state occupancy is inherently restricted to a geometric, memoryless distribution\.HSMMsgeneralize this framework by permitting each hidden state to persist for a variable interval, explicitly defined by a duration distribution\[[4](https://arxiv.org/html/2609.16500#bib.bib26)\]\. Consequently, the model incorporates three core probabilistic elements: state transitions, which govern the likelihood of moving between states; state\-specific durations, which model the time spent within a state; and observation emissions, which define the probability of an observation given the current state\.

HSMMs are particularly indispensable in computational biology for tasks such as genome annotation and segmentation\[[1](https://arxiv.org/html/2609.16500#bib.bib3),[5](https://arxiv.org/html/2609.16500#bib.bib23)\]\. Once a genome has been sequenced, a central challenge lies in identifying which nucleotide sequences correspond to protein\-coding exons, non\-coding introns, regulatory elements, or intergenic regions\. Interpreting the nucleotide sequence as an observation sequence, genome annotation amounts to inferring the hidden functional category of each region\. The key difficulty is that genomic features can span from a few dozen to many thousands of nucleotides; modeling such extended segments requires explicit duration distributions that bypass the constant transition probability inherent to HMMs, making HSMMs a natural fit\. This utility extends to chromatin state annotation,CpGisland detection, and other problems where segment lengths carry vital biological meaning\[[6](https://arxiv.org/html/2609.16500#bib.bib6),[7](https://arxiv.org/html/2609.16500#bib.bib35)\]\.

Given a sequence of observations, two central tasks arise:*learning*, which estimates the model parameters from observed data, and*decoding*, which identifies the most likely hidden state sequence\. Three fundamental algorithms solve these tasks: the*Forward–Backward*algorithm\[[8](https://arxiv.org/html/2609.16500#bib.bib2),[4](https://arxiv.org/html/2609.16500#bib.bib26)\]computes state probabilities at each time step, the*Baum–Welch*algorithm\[[9](https://arxiv.org/html/2609.16500#bib.bib1),[8](https://arxiv.org/html/2609.16500#bib.bib2),[4](https://arxiv.org/html/2609.16500#bib.bib26)\], an instance of Expectation–Maximization\[[10](https://arxiv.org/html/2609.16500#bib.bib33)\], estimates model parameters, and the*Viterbi*algorithm performs decoding\. The Viterbi algorithm directly solves the genome annotation problem; moreover, by iteratively applying Viterbi decoding and re\-estimating parameters from the decoded sequences, one can implement a Viterbi training loop, a technique widely used for*ab initio*gene prediction\[[11](https://arxiv.org/html/2609.16500#bib.bib24)\]when no prior data is available\. Since the Viterbi algorithm serves both decoding and*ab initio*learning, it is the most widely adopted of the three, making it a primary candidate for acceleration\.

However, while the Viterbi algorithm for a standard HMM has a complexity ofO⁡\(T​N2\)O\(TN^\{2\}\), the HSMM formulation introduces an additional loop over all possible durations, raising the complexity toO⁡\(T​N2​D\)O\(TN^\{2\}D\), whereDDis the maximum admissible duration\. In practical applications whereDDranges from hundreds to thousands\[[5](https://arxiv.org/html/2609.16500#bib.bib23)\], this extra dimension poses a major computational bottleneck\. Despite the critical importance of these problems, the algorithmic formulations used for HSMM inference have remained largely unchanged since their original proposal\.

Historically, the classical Viterbi algorithm for HSMMs has been implemented via four nested loops\[[4](https://arxiv.org/html/2609.16500#bib.bib26),[8](https://arxiv.org/html/2609.16500#bib.bib2),[12](https://arxiv.org/html/2609.16500#bib.bib36)\]\. Despite the proliferation of HPC resources, existing implementations remain confined to scalar, single\-threaded CPU execution\[[13](https://arxiv.org/html/2609.16500#bib.bib16),[14](https://arxiv.org/html/2609.16500#bib.bib17),[15](https://arxiv.org/html/2609.16500#bib.bib32),[16](https://arxiv.org/html/2609.16500#bib.bib31)\], in stark contrast to standard HMMs, which benefit from an extensive ecosystem of GPU\-accelerated and SIMD\-optimized frameworks\[[17](https://arxiv.org/html/2609.16500#bib.bib12),[18](https://arxiv.org/html/2609.16500#bib.bib7),[19](https://arxiv.org/html/2609.16500#bib.bib10),[20](https://arxiv.org/html/2609.16500#bib.bib38),[21](https://arxiv.org/html/2609.16500#bib.bib37)\]\.

Consequently, researchers face a restrictive trade\-off: either endure prohibitively long execution times or artificially truncate their datasets, preventing the full expressive potential of HSMMs from being realized in large\-scale scientific workflows\.

Given the importance of HSMMs and the widening gap between the computational demands of real\-world applications and the capabilities of existing sequential implementations, this work proposes the following contributions:

1. 1\.We introduce a novelTensor\-Based formulation of the Viterbi algorithm for Hidden Semi\-Markov Models\. By restructuring the traditional three\-nested inner loops into dense tensor operations, our approach naturally maps onto the SIMD and SPMD execution models of modern HPC architectures\. This reformulation exposes significant optimization opportunities that were previously inaccessible in sequential scalar implementations\.
2. 2\.We leverage this formulation to deliver optimized implementations across CPUs utilizing SIMD vectorization and threading, and, for the first time, GPUs\. To facilitate broader adoption, all implementations are released as a high\-performance open\-source library designed for seamless integration into existing scientific workflows\.
3. 3\.We conduct an extensive performance evaluation across three CPU and five GPU architectures, demonstrating speedups of up to570×570\\timesover state\-of\-the\-art HSMM frameworks\.

## IISequential HSMM Viterbi Algorithm

We now formalize the elements that define a Hidden Semi\-Markov Model\. An HSMM is represented by the parameter tuple\[[4](https://arxiv.org/html/2609.16500#bib.bib26),[8](https://arxiv.org/html/2609.16500#bib.bib2),[12](https://arxiv.org/html/2609.16500#bib.bib36)\]:

λ=\(𝒮,𝒪,𝐀,𝐏,𝐁,𝝅\),\\lambda=\\bigl\(\\mathcal\{S\},\\;\\mathcal\{O\},\\;\\mathbf\{A\},\\;\\mathbf\{P\},\\;\\mathbf\{B\},\\;\\boldsymbol\{\\pi\}\\bigr\)\\,,\(1\)where𝒮=\{s1,…,sN\}\\mathcal\{S\}=\\\{s\_\{1\},\\dots,s\_\{N\}\\\}is a finite set ofNNhidden states,𝒪=\{o1,…,oM\}\\mathcal\{O\}=\\\{o\_\{1\},\\dots,o\_\{M\}\\\}is a finite set ofMMobservation symbols, and𝝅N​\[j\]=πj\\boldsymbol\{\\pi\}^\{N\}\[j\]=\\pi\_\{j\}is the initial probability of statesjs\_\{j\}\. The collection of all transition probabilities forms the*transition probability matrix*𝐀N×N\\mathbf\{A\}^\{N\\times N\}, where each entry𝐀⁡\[i,j\]=ai​j\\mathbf\{A\}\[i,\\,j\]=a\_\{ij\}represents the probability of transitioning from statesis\_\{i\}to statesjs\_\{j\}\. Similarly, the collection of all duration probabilities forms the*duration probability matrix*𝐏N×D\\mathbf\{P\}^\{N\\times D\}, where each entry𝐏​\[j,d\]=pj​\(d\)\\mathbf\{P\}\[j,\\,d\]=p\_\{j\}\(d\)represents the probability that statesjs\_\{j\}persists for exactlyddconsecutive time\-steps, withd∈\{1,…,D\}d\\in\\\{1,\\dots,D\\\}\. Finally, the emission probabilities define the*emission probability matrix*𝐁N×M\\mathbf\{B\}^\{N\\times M\}, where𝐁​\[j,o\]=bj​\(o\)\\mathbf\{B\}\[j,\\,o\]=b\_\{j\}\(o\)is the probability that statesjs\_\{j\}emits observationoo\. These matrix definitions of the model parameters will be central to the tensor formulation presented in Sec\.[III](https://arxiv.org/html/2609.16500#S3)\. Table[I](https://arxiv.org/html/2609.16500#S2.T1)summarizes all components of the model\.

TABLE I:Components of an HSMM and Viterbi AlgorithmThe Viterbi algorithm is one of the fundamental algorithms for HSMMs\. It addresses the following problem: given a modelλ\\lambdaand a sequence of observationso0,…,oT−1o\_\{0\},\\dots,o\_\{T\-1\}, the goal is to find the most likely hidden state sequence𝐬∗=\(q0∗,…,qT−1∗\)\\mathbf\{s\}^\{\*\}=\(q^\{\*\}\_\{0\},\\dots,q^\{\*\}\_\{T\-1\}\)\. To achieve this, for every time stepttand each statesj∈𝒮s\_\{j\}\\in\\mathcal\{S\}, the algorithm computesδt​\(j\)\\delta\_\{t\}\(j\), which represents the likelihood of the most probable state sequence ending in statesjs\_\{j\}at timett\. Both the sequential and the tensor formulations \(Sec\.[III](https://arxiv.org/html/2609.16500#S3)\) maintain these values in aΔN×T\{\\Delta\}^\{N\\times T\}matrix:

𝚫N×T=\[δ1​\(0\)δ1​\(1\)⋯δ1​\(T−1\)δ2​\(0\)δ2​\(1\)⋯δ2​\(T−1\)⋱δN​\(0\)δN​\(1\)⋯δN​\(T−1\)\]\.\\small\\boldsymbol\{\\Delta\}^\{N\\times T\}\\;=\\;\\begin\{bmatrix\}\\delta\_\{1\}\(0\)&\\delta\_\{1\}\(1\)&\\cdots&\\delta\_\{1\}\(T\\\!\-\\\!1\)\\\\\[2\.0pt\] \\delta\_\{2\}\(0\)&\\delta\_\{2\}\(1\)&\\cdots&\\delta\_\{2\}\(T\\\!\-\\\!1\)\\\\\[2\.0pt\] \\vdots&\\vdots&\\ddots&\\vdots\\\\\[2\.0pt\] \\delta\_\{N\}\(0\)&\\delta\_\{N\}\(1\)&\\cdots&\\delta\_\{N\}\(T\\\!\-\\\!1\)\\end\{bmatrix\}\.
Following a dynamic programming approach, the algorithm proceeds in three distinct stages:initialization,induction, andbacktracking\.

### II\-AInitialization Phase

The initialization phase covers the firstDDtime steps, whereDDis the maximum admissible state duration\. During this interval, we must account for the possibility that the system has occupied statesjs\_\{j\}sincet=1t=1with no prior transition\. For each statesjs\_\{j\}and each1≤t≤D1\\leq t\\leq D, the initialization value combines\(a\)the initial state probabilityπj\\pi\_\{j\},\(b\)the probability that statesjs\_\{j\}persists for exactlytttime steps, and\(c\)the joint emission probability of observationso0,…,ot−1o\_\{0\},\\dots,o\_\{t\-1\}under statesjs\_\{j\}\. Fort\>Dt\>D, no state can have persisted since the beginning, so this contribution is no longer considered\.

δt​\(j\)=πj⏟\(a\)⋅pj​\(t\)⏟\(b\)⋅∏τ=0t−1bj​\(ot−τ\)⏟\(c\),1≤t≤D\\delta\_\{t\}\(j\)=\\underbrace\{\\hbox\{\\pagecolor\{initbg\}$\\displaystyle\\pi\_\{j\}$\}\}\_\{\\text\{\(a\)\}\}\\;\\cdot\\;\\underbrace\{\\hbox\{\\pagecolor\{durbg\}$\\displaystyle p\_\{j\}\(t\)$\}\}\_\{\\text\{\(b\)\}\}\\;\\cdot\\;\\underbrace\{\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\prod\_\{\\tau=0\}^\{t\-1\}b\_\{j\}\(o\_\{t\-\\tau\}\)$\}\}\_\{\\text\{\(c\)\}\}\\;,\\quad 1\\leq t\\leq D\(2\)

### II\-BInduction phase

After initializing the firstDDtime steps, we proceed with the most computationally intensive phase: the induction\. For each time stepttand each current statesjs\_\{j\}, the goal is to find the previous statesis\_\{i\}and the durationddthat together maximize the likelihood of reachingsjs\_\{j\}at timett\. The formulation is given by Equation \([3](https://arxiv.org/html/2609.16500#S2.E3)\), and the corresponding four\-nested\-loop pseudocode is shown in Algorithm[1](https://arxiv.org/html/2609.16500#alg1)\.

The computation combines four factors\. Term\(a\)is the previously computed valueδt−d​\(i\)\\delta\_\{t\-d\}\(i\): it encodes the likelihood of the best path ending in statesis\_\{i\}at timet−dt\\\!\-\\\!d, under the assumption that a transition tosjs\_\{j\}occurred there\. Term\(b\)is the transition probabilityai​ja\_\{ij\}from statesis\_\{i\}to statesjs\_\{j\}\.

Together,\(a\)and\(b\)form the inner maximization: for a fixed durationdd, we evaluate all possible source statessis\_\{i\}and select the one that yields the highest likelihood\.

The result is then multiplied by term\(c\), the probabilitypj​\(d\)p\_\{j\}\(d\)that statesjs\_\{j\}persists for exactlyddconsecutive time steps, and by term\(d\), the cumulative emission probability of all observations fromt−d\+1t\-d\+1tottunder statesjs\_\{j\}\. The maximization repeats this for all durationsd∈\{1,…,min⁡\(t,D\)\}d\\in\\\{1,\\dots,\\min\(t,D\)\\\}and all source statessis\_\{i\}and selects the best combination\. The resulting optimal combination\(d∗,i∗\)\(d^\{\*\},\\,i^\{\*\}\)for each statesjs\_\{j\}at time stepttis stored inψt​\(j\)\\psi\_\{t\}\(j\), while the corresponding likelihood is stored inδt​\(j\)\\delta\_\{t\}\(j\)\.

δt​\(j\)=max∀i,∀d⁡\[δt−d​\(i\)⏟\(a\)⋅ai​j⏟\(b\)⋅pj​\(d\)⏟\(c\)⋅∏k=0d−1bj​\(ot−k\)⏟\(d\)\],\\delta\_\{t\}\(j\)=\\max\_\{\\forall i,\\forall d\}\\Bigl\[\\underbrace\{\\hbox\{\\pagecolor\{deltabg\}$\\displaystyle\\delta\_\{t\-d\}\(i\)$\}\}\_\{\(a\)\}\\cdot\\underbrace\{\\hbox\{\\pagecolor\{transbg\}$\\displaystyle a\_\{ij\}$\}\}\_\{\(b\)\}\\cdot\\underbrace\{\\hbox\{\\pagecolor\{durbg\}$\\displaystyle p\_\{j\}\(d\)$\}\}\_\{\(c\)\}\\cdot\\underbrace\{\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\prod\_\{k=0\}^\{d\-1\}b\_\{j\}\(o\_\{t\-k\}\)$\}\}\_\{\(d\)\}\\Bigr\],\(3\)
Algorithm 1Sequential HSMM Viterbi Induction\.1:for

t=2t=2to

TTdo

2:for

j=1j=1to

NNdo

3:

δt​\(j\)←−∞\\delta\_\{t\}\(j\)\\leftarrow\-\\infty
4:for

d=1d=1to

min⁡\(t,D\)\\min\(t,\\,D\)do

5:for

i=1i=1to

NNdo

6:val

←\\leftarrowδt−d​\(i\)\\delta\_\{t\-d\}\(i\)

⋅\\cdotai​ja\_\{ij\}

⋅\\cdotpj​\(d\)p\_\{j\}\(d\)

⋅\\cdot∏k=0d−1bj​\(ot−k\)\\prod\_\{k=0\}^\{d\-1\}b\_\{j\}\(o\_\{t\-k\}\)

7:ifval

\>δt​\(j\)\>\\delta\_\{t\}\(j\)then

8:

δt​\(j\)←\\delta\_\{t\}\(j\)\\leftarrowval

9:

ψt​\(j\)←\(d,i\)\\psi\_\{t\}\(j\)\\leftarrow\(d,\\,i\)

Note that induction starts fromt=2t=2; fort≤Dt\\leq D, the valueδt​\(j\)\\delta\_\{t\}\(j\)is determined by the maximum of two cases: the system has remained in statejjsincet=1t=1\(initialization\), or it transitioned from a previous stateiiat some timeτ<t\\tau<t\(induction\)\. The algorithm takes the greater of these two probabilities\.

### II\-CBacktracking Phase

Onceδt​\(j\)\\delta\_\{t\}\(j\)has been computed∀t,∀j\\forall\\,t,\\;\\forall\\,j, we can recover the optimal state sequence\. Starting from the last time step, we select the state with the highest delta value\. The backtracking then proceeds backwards, untilt=0t=0, returning the most likely chain of states\. Since we operate in a Semi\-Markov regime, the recovered path will typically exhibit states persisting across multiple consecutive time steps, reflecting the explicit duration modeling that distinguishes the HSMM from an HMM\.

## IIITensor\-Based Viterbi Algorithm

The tensor formulation of the Viterbi algorithm follows the same subdivision as the standard one: an*initialization phase*, an*induction phase*, and a*backtracking phase*\. However in this work, we express both initialization and induction using a new tensor formulation\. The backtracking phase remains sequential as it does not represent a computational bottleneck\. Instead, optimization efforts focus on the induction phase, where tensor reformulation yields the greatest benefit\.

### From Loops to Tensors

As described in the sequential Algorithm[1](https://arxiv.org/html/2609.16500#alg1), for every time stepttwe need to evaluate all possible combinations of a previous statesis\_\{i\}and a durationddfor each current statesjs\_\{j\}\. Formally, for every statesj∈𝒮s\_\{j\}\\in\\mathcal\{S\}we must consider every pair\(si,d\)\(s\_\{i\},\\,d\)withsi∈𝒮∖\{sj\}s\_\{i\}\\in\\mathcal\{S\}\\setminus\\\{s\_\{j\}\\\}andd∈\{1,…,D\}d\\in\\\{1,\\dots,D\\\}\. Once all possible combinations of\(si,d\)\(s\_\{i\},d\)have been computed for each statesj∈𝒮s\_\{j\}\\in\\mathcal\{S\}, our aim is to computeδt​\(j\)\\delta\_\{t\}\(j\)for every time stept≤Tt\\leq T\.

To compute these combinations and thenδt​\(j\),∀j\\delta\_\{t\}\(j\),\\forall jsequentially, one must iterate over the current states, then over the possible durations, and then again over the previous states, yieldingthree nested loopsalready contained inside the main time\-step loop\. The sequential algorithm therefore requires four nested loops\.

### III\-AThe Brick Representation

![Refer to caption](https://arxiv.org/html/2609.16500v1/brick_representation.png)Fig\. 1:Visualization of theBrick3D tensor and subdivision intos​jsj\-slices\.Our key idea is to replace the three nested loops with structured tensor operations that make data reuse explicit and expose independent computations along each axis\. The loops over\(sj,si,d\)\(s\_\{j\},s\_\{i\},d\)can be naturally mapped onto a three\-dimensional tensor, as shown in Figure[1](https://arxiv.org/html/2609.16500#S3.F1)a\. We choose the following layout:

- •yy\-axis⟶\\;\\longrightarrow\\;target statessjs\_\{j\}, withj∈\{1,…,N\}j\\in\\\{1,\\dots,N\\\};
- •xx\-axis⟶\\;\\longrightarrow\\;source statessis\_\{i\}, withi∈\{1,…,N\}i\\in\\\{1,\\dots,N\\\};
- •zz\-axis⟶\\;\\longrightarrow\\;durationsdd, withd∈\{1,…,D\}d\\in\\\{1,\\dots,D\\\}\.

With this convention the set of all combinations\(sj,si,d\)\(s\_\{j\}\\,,s\_\{i\}\\,,d\)can be visualized as a 3\-dimensional tensor of sizeN×N×DN\\times N\\times Dthat we will callBrick\(ℬ\\mathcal\{B\}\), shown in Figure[1](https://arxiv.org/html/2609.16500#S3.F1)a\. Each*slice*of theBrickalong theyy\-axis corresponds to the complete set of\(si,d\)\(s\_\{i\},\\,d\)combinations for a single target statesjs\_\{j\}, which we callsjs\_\{j\}\-slices \(Figure[1](https://arxiv.org/html/2609.16500#S3.F1)b\)\. Using this representation, we reformulate the sequential induction of Algorithm[1](https://arxiv.org/html/2609.16500#alg1)into three key stages: \(i\) theBrick Construction, which occurs once \(Sec\.[III\-B](https://arxiv.org/html/2609.16500#S3.SS2)\); \(ii\) theBrick Update, where theBrickis modified at each time stepttto incorporate temporal factors \(as described in Sec\.[III\-C](https://arxiv.org/html/2609.16500#S3.SS3)\); and \(iii\) theMaximum Extraction, which identifies the maximum value and the corresponding coordinates\(si,d\)\(s\_\{i\},d\)within eachsjs\_\{j\}\-slice\. While phases \(ii\) and \(iii\) are executed iteratively at each time step, phase \(i\) is a pre\-computation step performed outside the time steps loop\. This pipeline, including initialization, takes shape within the Algorithm[2](https://arxiv.org/html/2609.16500#alg2)that will be described line\-by\-line below\.

Algorithm 2Tensor\-based HSMM Viterbi\.0:Initialization\(1≤t≤D\)\(1\\leq t\\leq D\)

1:

𝐄N×D←Emission Product Computation\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\mathbf\{E\}^\{N\\times D\}$\}\\leftarrow\\textit\{Emission Product Computation\}
2:

𝚫N×D←𝝅N×↑⊙𝐏N×D⊙𝐄N×D\\boldsymbol\{\\Delta\}^\{N\\times D\}\\leftarrow\\hbox\{\\pagecolor\{initbg\}$\\displaystyle\\boldsymbol\{\\pi\}^\{N\\times\\uparrow\}$\}\\odot\\hbox\{\\pagecolor\{durbg\}$\\displaystyle\\mathbf\{P\}^\{N\\times D\}$\}\\odot\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\mathbf\{E\}^\{N\\times D\}$\}
2:Induction\(2≤t≤T\)\(2\\leq t\\leq T\)

3:

ℬf​i​r​s​tN×N×D←𝐀N×N×↑⊙𝐏N×↑×D\\mathcal\{B\}\_\{first\}^\{N\\times N\\times D\}\\leftarrow\\hbox\{\\pagecolor\{transbg\}$\\displaystyle\\mathbf\{A\}^\{N\\times N\\times\\uparrow\}$\}\\odot\\hbox\{\\pagecolor\{durbg\}$\\displaystyle\\mathbf\{P\}^\{N\\times\\uparrow\\times D\}$\}// Brick Construction

4:for

t=2t=2to

TTdo

5:

𝐄N×D←Emission Product Computation\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\mathbf\{E\}^\{N\\times D\}$\}\\leftarrow\\textit\{Emission Product Computation\}
6:

𝚫p​a​s​tN×D←𝚫​\(t−D:t−1\)N×D\\hbox\{\\pagecolor\{deltabg\}$\\displaystyle\\boldsymbol\{\\Delta\}\_\{past\}^\{N\\times D\}$\}\\leftarrow\\boldsymbol\{\\Delta\}\(t\-D:t\-1\)^\{N\\times D\}// Past Delta Extraction

7:

ℬ\(t\)N×N×D←ℬf​i​r​s​t⊙𝚫p​a​s​t↑×N×D⊙𝐄N×↑×D\\hbox\{\\pagecolor\{brickbg\}$\\displaystyle\\mathcal\{B\}\_\{\(t\)\}^\{N\\times N\\times D\}$\}\\leftarrow\\mathcal\{B\}\_\{first\}\\odot\\hbox\{\\pagecolor\{deltabg\}$\\displaystyle\\boldsymbol\{\\Delta\}\_\{past\}^\{\\uparrow\\times N\\times D\}$\}\\odot\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\mathbf\{E\}^\{N\\times\\uparrow\\times D\}$\}// Brick Update

8:for

j=1j=1to

NNdo

9:

Ψj​\(t\)←arg⁡maxd,i⁡ℬ\(t\)​\[sj​\-slice\]N×D\\Psi\_\{j\}\(t\)\\leftarrow\\arg\\max\_\{d,\\,i\}\\;\\hbox\{\\pagecolor\{brickbg\}$\\displaystyle\\mathcal\{B\}\_\{\(t\)\}\[s\_\{j\}\\text\{\-slice\}\]^\{N\\times D\}$\}// Max\. Extraction

9:Backtracking

10:

𝐪∗=Backtracking​\(ΨN×T,ΔN×T\)\\mathbf\{q\}^\{\*\}=\\textit\{Backtracking\}\(\\Psi^\{N\\times T\},\\;\\Delta^\{N\\times T\}\)

### III\-BBrick Construction

The key observation is the following: thetransition probability matrix𝐀N×N\\mathbf\{A\}^\{N\\times N\}and theduration probability matrix𝐏N×D\\mathbf\{P\}^\{N\\times D\}do not depend on the time steptt\. Their contribution to theBrickcan therefore be precomputed once, before entering the time\-step loop\. This corresponds to line 3 of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\. As shown in Figure[2](https://arxiv.org/html/2609.16500#S3.F2)a, the transition matrix𝐀\\mathbf\{A\}is a two\-dimensionalN×NN\\times Nmatrix\. Within ourBrick, it occupies the*front face*, i\.e\. it is aligned with theyy\- andxx\-axes\. The duration probability matrix𝐏\\mathbf\{P\}is a two\-dimensionalN×DN\\times Dmatrix positioned on the*side face*, i\.e\. aligned with theyy\- andzz\-axes \(Figure[2](https://arxiv.org/html/2609.16500#S3.F2)a\)\.

![Refer to caption](https://arxiv.org/html/2609.16500v1/broadcasted_product.png)Fig\. 2:Alignment of𝐀\\mathbf\{A\}and𝐏\\mathbf\{P\}within theBrick3D tensor and the corresponding representation of the broadcasted product\.To combine these two matrices into a single three\-dimensional tensor, we perform the following product:

ℬf​i​r​s​tN×N×D=𝐀~N×N×↑⏟\(a\)⊙𝐏~N×↑×D⏟\(b\)\\mathcal\{B\}\_\{first\}^\{N\\times N\\times D\}=\\underbrace\{\\hbox\{\\pagecolor\{transbg\}$\\displaystyle\\widetilde\{\\mathbf\{A\}\}^\{N\\times N\\times\\uparrow\}$\}\}\_\{\(a\)\}\\;\\odot\\;\\underbrace\{\\hbox\{\\pagecolor\{durbg\}$\\displaystyle\\widetilde\{\\mathbf\{P\}\}^\{N\\times\\uparrow\\times D\}$\}\}\_\{\(b\)\}\(4\)where the↑\\uparrowsymbol denotes the axis along which broadcasting occurs, expanding the tensor to match the dimensions of the corresponding operand\. Concretely, as shown in Figure[2](https://arxiv.org/html/2609.16500#S3.F2)b\-c, the operation can be understood in two steps:

\(a\)Broadcast𝐀\\mathbf\{A\}: replicate theN×NN\\times NmatrixDDtimes along thezz\-axis, obtaining a tensor𝐀~N×N×D\\widetilde\{\\mathbf\{A\}\}^\{N\\times N\\times D\}\(Figure[2](https://arxiv.org/html/2609.16500#S3.F2)b\-\(a\)\)\.

\(b\)Broadcast𝐏\\mathbf\{P\}: replicate theN×DN\\times DmatrixNNtimes along thexx\-axis, obtaining a tensor𝐏~N×N×D\\widetilde\{\\mathbf\{P\}\}^\{N\\times N\\times D\}\(Figure[2](https://arxiv.org/html/2609.16500#S3.F2)b\-\(b\)\)\.

Then, theBrickis given by the element\-wise product of the two broadcast matrices:ℬf​i​r​s​tN×N×D=𝐀~N×N×D⊙𝐏~N×N×D\\mathcal\{B\}\_\{first\}^\{N\\times N\\times D\}=\\widetilde\{\\mathbf\{A\}\}^\{N\\times N\\times D\}\\odot\\widetilde\{\\mathbf\{P\}\}^\{N\\times N\\times D\}\(Figure[2](https://arxiv.org/html/2609.16500#S3.F2)c\)\. We refer to this as a broadcasted product: a fundamental operation of our tensor\-based Viterbi algorithm\.

### III\-CBrick Update

After constructing the staticBrickℬf​i​r​s​t\\mathcal\{B\}\_\{first\}, we enter the main loop over time\-steps\. At each steptt, two additional matrices must be incorporated: the*past delta values*𝚫p​a​s​t\\boldsymbol\{\\Delta\}\_\{past\}and the*emission probability product*𝐄\\mathbf\{E\}, both of which depend strictly ontt\. This phase corresponds to lines 5–7 of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\.

#### III\-C1Extracting Past Delta Matrix

The past delta values are the simpler of the two time\-dependent matrices\. We collect from the stored𝚫\\boldsymbol\{\\Delta\}matrix a window of sizeN×DN\\times Dobtaining𝚫p​a​s​tN×D\\boldsymbol\{\\Delta\}\_\{past\}^\{N\\times D\}\.

𝚫p​a​s​tN×D=\[δ1​\(t−D\)⋯δ1​\(t−1\)δ2​\(t−D\)⋯δ2​\(t−1\)δN​\(t−D\)⋯δN​\(t−1\)\]\.\\boldsymbol\{\\Delta\}\_\{past\}^\{N\\times D\}\\;=\\;\\begin\{bmatrix\}\\delta\_\{1\}\(t\-D\)&\\cdots&\\delta\_\{1\}\(t\-1\)\\\\\[2\.0pt\] \\delta\_\{2\}\(t\-D\)&\\cdots&\\delta\_\{2\}\(t\-1\)\\\\\[2\.0pt\] \\vdots&\\vdots&\\vdots\\\\\[2\.0pt\] \\delta\_\{N\}\(t\-D\)&\\cdots&\\delta\_\{N\}\(t\-1\)\\end\{bmatrix\}\.As shown in Figure[3](https://arxiv.org/html/2609.16500#S3.F3), this window starts att−Dt\-D\(or11ift<Dt<D\) and ends att−1t\-1, reversed so that the first column corresponds to the most recent past step \(durationd=1d=1\)\. The resulting matrix𝚫p​a​s​tN×D\\boldsymbol\{\\Delta\}\_\{past\}^\{N\\times D\}is aligned on thexx\- andzz\-axes of theBrick\(lower face\)\.

![Refer to caption](https://arxiv.org/html/2609.16500v1/extracting_pastdeltas.png)Fig\. 3:Extracting from𝚫\\boldsymbol\{\\Delta\}past\[t−D,t\)\[t\\\!\-\\\!D,t\)values to obtain𝚫p​a​s​t\\boldsymbol\{\\Delta\}\_\{past\}\.
#### III\-C2Computing Emission Probability Matrix

The emission probability matrix is the more complex of the two factors\. Recall from the term \(d\) of Equation[3](https://arxiv.org/html/2609.16500#S2.E3)that, for a given statesjs\_\{j\}and durationdd, we need the product of the emission probabilities over theddmost recent observations\. In the sequential algorithm this partial product is trivially computed inside the duration loop\. In the tensor formulation, however, the duration loop has been eliminated\. We therefore need to compute the entireN×DN\\times Dmatrix of cumulative emission at each time step\.

𝐄N×D=\[b1​\(ot\)∏k=01b1​\(ot−k\)…∏k=0d−1b1​\(ot−k\)b2​\(ot\)∏k=01b2​\(ot−k\)…∏k=0d−1b2​\(ot−k\)⋱bN​\(ot\)∏k=01bN​\(ot−k\)…∏k=0d−1bN​\(ot−k\)\]\.\\mathbf\{E\}^\{N\\times D\}=\\begin\{bmatrix\}b\_\{1\}\(o\_\{t\}\)&\\prod\_\{k=0\}^\{1\}b\_\{1\}\(o\_\{t\-k\}\)&\\dots&\\prod\_\{k=0\}^\{d\-1\}b\_\{1\}\(o\_\{t\-k\}\)\\\\\[1\.0pt\] b\_\{2\}\(o\_\{t\}\)&\\prod\_\{k=0\}^\{1\}b\_\{2\}\(o\_\{t\-k\}\)&\\dots&\\prod\_\{k=0\}^\{d\-1\}b\_\{2\}\(o\_\{t\-k\}\)\\\\\[1\.0pt\] \\vdots&\\vdots&\\ddots&\\vdots&\\\\\[1\.0pt\] b\_\{N\}\(o\_\{t\}\)&\\prod\_\{k=0\}^\{1\}b\_\{N\}\(o\_\{t\-k\}\)&\\dots&\\prod\_\{k=0\}^\{d\-1\}b\_\{N\}\(o\_\{t\-k\}\)\\end\{bmatrix\}\.
As shown in Figure[4](https://arxiv.org/html/2609.16500#S3.F4), we extract theDDmost recent observation indices, look up the corresponding emission prob for all states inside𝐁N×M\\mathbf\{B\}^\{N\\times M\}, reverse the order \(so that the first column corresponds to durationd=1d=1\)\. Then, we apply a cumulative product along the duration axis obtaining𝐄N×D\\mathbf\{E\}^\{N\\times D\}\(Figure[4](https://arxiv.org/html/2609.16500#S3.F4)\)\. Our tensor formulation highlights significant redundancies, leading us to propose an optimized caching strategy as detailed in Sec\.[IV](https://arxiv.org/html/2609.16500#S4)\.

Fig\. 4:Computing Emission Product to obtain the Emission Probability Matrix\.Once the past delta matrix𝚫p​a​s​t\\boldsymbol\{\\Delta\}\_\{past\}and the emission probability matrix𝐄\\mathbf\{E\}have been computed for a given time steptt, we combine them with the precomputedBrickℬ\\mathcal\{B\}through two successive broadcasted products:

ℬ\(t\)N×N×D=ℬf​i​r​s​tN×N×D⊙𝚫~p​a​s​t​\(t\)↑×N×D⏟\(c\)⊙𝐄~\(t\)N×↑×D⏟\(d\)\\hbox\{\\pagecolor\{brickbg\}$\\displaystyle\\mathcal\{B\}\_\{\(t\)\}^\{N\\times N\\times D\}$\}\\;=\\;\\mathcal\{B\}\_\{first\}^\{N\\times N\\times D\}\\;\\odot\\;\\underbrace\{\\hbox\{\\pagecolor\{deltabg\}$\\displaystyle\\widetilde\{\\boldsymbol\{\\Delta\}\}\_\{past~\(t\)\}^\{\\uparrow\\times N\\times D\}$\}\}\_\{\(c\)\}\\;\\odot\\;\\underbrace\{\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\widetilde\{\\mathbf\{E\}\}\_\{\(t\)\}^\{N\\times\\uparrow\\times D\}$\}\}\_\{\(d\)\}\(5\)
\(c\)𝚫p​a​s​t\\boldsymbol\{\\Delta\}\_\{past\}lies on thexx–zzplane, and broadcast theN×DN\\times DmatrixNNtimes along theyy\-axis, obtaining a tensor𝚫~p​a​s​tN×N×D\\widetilde\{\\boldsymbol\{\\Delta\}\}\_\{past\}^\{N\\times N\\times D\}\.

\(d\)𝐄\\mathbf\{E\}lies on theyy–zzaxis, and broadcast theN×DN\\times DmatrixNNtimes along thexx\-axis, obtaining a tensor𝐄~N×N×D\\widetilde\{\\mathbf\{E\}\}^\{N\\times N\\times D\}\.

The resulting tensorℬN×N×D\\mathcal\{B\}^\{N\\times N\\times D\}is the fully populatedBrickfor time steptt: each entryℬ⁡\[j,i,d\]\\mathcal\{B\}\[j,\\,i,\\,d\]encodes the likelihood of transitioning from statesis\_\{i\}to statesjs\_\{j\}with durationdd, given the observations up to timett\. Having obtained all combinations of\(si,d\)\(s\_\{i\},d\)for eachsjs\_\{j\}, we must now identify the one that yields the maximum likelihood for eachsjs\_\{j\}\.

### III\-DMaximum Extraction

We now need to extract, for each target statesjs\_\{j\}, the combination\(si∗,d∗\)\(s\_\{i\}^\{\*\},\\,d^\{\*\}\)that maximizes the tensor entry\. Recall thatsjs\_\{j\}is indexed along theyy\-axis; the corresponding slice is therefore a two\-dimensionalN×DN\\times Dmatrix\. For eachsjs\_\{j\}\-slice we seek the maximum value and its associated coordinates\(si,d\)\(s\_\{i\},d\)as shown in Figure[5](https://arxiv.org/html/2609.16500#S3.F5)\. This yields the resulting maximum values vector𝚫​\(t\)N\\mathbf\{\\Delta\}\(t\)^\{N\}\. This phase corresponds to lines 8–9 of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\.

Fig\. 5:Argmax computation in eachsjs\_\{j\}\-slice\.Seeking the maximum value is a fundamental problem in parallel computing and several strategies can be employed to improve its performance; we discuss our strategy in Sec\.[IV](https://arxiv.org/html/2609.16500#S4)\.

### III\-EInitialization Phase

As described in Sec\.[II\-A](https://arxiv.org/html/2609.16500#S2.SS1), the initialization phase covers the firstDDtime steps \(1≤t≤D1\\leq t\\leq D\), accounting for the possibility that statesjs\_\{j\}has persisted sincet=0t=0with no prior transition\. In the tensor formulation, the entire initialization is expressed as a singlebroadcasted productof three matrices \(lines 1–2 of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\):

𝚫\[∀j,:D\]N×D=𝝅N×↑⊙𝐏N×D⊙𝐄N×D\\boldsymbol\{\\Delta\}\[\\forall j,\\,1\\\!:\\\!D\]^\{N\\times D\}=\\hbox\{\\pagecolor\{initbg\}$\\displaystyle\\boldsymbol\{\\pi\}^\{N\\times\\uparrow\}$\}\\odot\\hbox\{\\pagecolor\{durbg\}$\\displaystyle\\mathbf\{P\}^\{N\\times D\}$\}\\odot\\hbox\{\\pagecolor\{embg\}$\\displaystyle\\mathbf\{E\}^\{N\\times D\}$\}\(6\)Here,\(a\)𝝅N\\boldsymbol\{\\pi\}^\{N\}is the initial state probability vector, broadcasted along the duration axis;\(b\)𝐏N×D\\mathbf\{P\}^\{N\\times D\}is the duration probability matrix; and\(d\)𝐄N×D\\mathbf\{E\}^\{N\\times D\}is the cumulative emission matrix, where each entry𝐄⁡\[j,d\]\\mathbf\{E\}\[j,d\]accumulates the emission probabilities of the firstddobservations under statesjs\_\{j\}\. Since no entry depends on any other, allN×DN\\times Dvalues are computed with no loop\-carried dependencies\. During these firstDDsteps, both initialization and induction contribute; the finalδt​\(j\)\\delta\_\{t\}\(j\)is taken as the maximum of the two \(lines 8–9 of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\)\.

## IVImplementations and Optimizations

This section details the optimization strategies for the Tensor\-Based Viterbi algorithm from Sec\.[III](https://arxiv.org/html/2609.16500#S3)\. The tensor formulation naturally enables several optimizations: it exposes data reuse patterns for cache\-friendly access and structures computation along independent axes, mapping efficiently onto SIMD and multi\-threaded execution models\.

We developed four versions of the Tensor Viterbi algorithm:Tens\-Py,Tens\-1c,Tens\-mcandTens\-gpu\.Tens\-Pyis a direct transcription of Algorithm[2](https://arxiv.org/html/2609.16500#alg2)in NumPy, and it was the first implementation developed\. No low\-level optimization is attempted; the implementation serves as a readable reference and validation baseline\. It was necessary to analyze the algorithm, identify optimization opportunities, and test them before delving into the low\-level optimized implementations\.

We highlight that all operations are implemented in log\-space\. Since the quantities involved lie in the interval\[0,1\]\[0,\\,1\], repeated multiplications would quickly lead to numerical underflow\. Working in log\-space is a standard practice adopted by all major Markov model frameworks\[[8](https://arxiv.org/html/2609.16500#bib.bib2),[13](https://arxiv.org/html/2609.16500#bib.bib16)\]\. The practical consequence is straightforward: products become sums, cumulative products become cumulative sums, and broadcasted products become broadcasted sums\.

TABLE II:Summary of Implementations### IV\-ANumPy Version

Tens\-Pyfollows the tensor algorithm precisely and is implemented using NumPy, which allows the algorithm to be expressed directly in terms of tensor operations without requiring further low\-level optimizations\. The memory layout is the one used by NumPy \(row\-major, depth\-first\)\. The broadcasted sums are expressed in two phases: the 2D matrix is broadcast over the absent axis to obtain a 3D tensor, then the two 3D tensors are summed element\-wise\. This pattern is applied to both theBrick ComputationandBrick Updatephases\. ThePast Delta Extractionphase is performed using NumPy slicing operations, and the argmax is computed using NumPy’sargmaxfunction independently for each destination statesjs\_\{j\}\.

Once this first version was coded, analysis of the prototype revealed an optimization opportunity: the emission accumulation, originally recomputed inside the duration loop at every time\-step, can be decoupled from the main recurrence and maintained through a rolling cache\. After the warm\-up phase \(t\>Dt\>D\), the full emission buffer is obtained by a single element\-wise addition and a cache shift, reducing the per\-time\-step cost fromO⁡\(D​N\)O\(DN\)to amortizedO⁡\(N\)O\(N\)and eliminating the loop\-carried dependency on the duration axis\. This optimization has been critical in the C\+\+ versions\.

### IV\-BCPU Implementation

Tens\-1candTens\-mcare developed in C\+\+ without relying on any tensor library, since neither BLAS\[[22](https://arxiv.org/html/2609.16500#bib.bib29)\]nor frameworks such as xTensor\[[23](https://arxiv.org/html/2609.16500#bib.bib30)\]provide broadcasted sums natively\. Implementing the operations explicitly also enabled us to fuse theBrick UpdateandArgmaxphases, avoiding the need to store the fully populatedBrickin memory before computing the maximum\.

These versions allow us to exploit the optimization opportunities unveiled by the tensor formulation and produce tuned variants of the proposed algorithm\.

#### IV\-B1Memory Access Patterns

We introduced two complementary flat layouts designed for spatial locality\. TheΔp​a​s​t\\Delta\_\{past\}buffer uses a time\-major layout\(t⋅N\+j\)\(t\\cdot N\+j\)so that all states at a given time\-step are contiguous\. TheΔ\\DeltaandΨ\\Psiarrays use a state\-major layout\(j⋅T\+t\)\(j\\cdot T\+t\)as required by theBacktracking Phase\. TheBricktensor uses a\(j,d,i\)\(j,d,i\)layout where, for a fixed statejj, the entireD×ND\\times Nblock is contiguous, fitting in L2 across the duration loop if it is small enough\. This is a deliberate choice: the hot inner loop sweeps over statessis\_\{i\}within a fixed\(sj,d\)\(s\_\{j\},d\)slice, achieving contiguous access\.

#### IV\-B2Cached Emissions Computing

We build a 2D emission buffer𝐄\\mathbf\{E\}indexed by\(d⋅N\+j\)\(d\\cdot N\+j\), storing the cumulative sum of emission log\-probabilities over each observation window\. Fort≤Dt\\leq D, the cumulative sum is computed from scratch\. Fort\>Dt\>D, the buffer is updated incrementally using anemission cache𝐄c​a​c​h​e\\mathbf\{E\}\_\{cache\}: the new observation log\-probability at timettis added to a right\-shifted copy of𝐄c​a​c​h​e\\mathbf\{E\}\_\{cache\}, and the result is saved back to𝐄c​a​c​h​e\\mathbf\{E\}\_\{cache\}for the next time\-step\. This reduces the per\-state emission update fromO⁡\(D\)O\(D\)to amortizedO⁡\(1\)O\(1\)after the warm\-up phase, and makes the emission values for all\(d,j\)\(d,j\)pairs available independently before theArgmaxloop begins\.

#### IV\-B3Fused Brick Update and Maximum Extraction

The two most intensive phases of the algorithm can be fused in this implementation, performing a fusedBrick Update\+\\,\+\\,Argmaxwith no self\-transition exclusion \(the transition matrix encodes this structurally with zero\-values\)\. The inner loop is entirely branchless, using annotated ternary conditional assignments\. This branchless pattern allows the compiler to generate predicated instructions rather than conditional branches, eliminating branch misprediction penalties\. In the three\-level loop\(sj,d,si\)\(s\_\{j\},d,s\_\{i\}\), iterations are fully independent across all three axes, making each axis independently parallelizable or vectorizable\.

#### IV\-B4Multi\-Core Implementation

TheTens\-mcimplementation parallelizes the single\-core version using OpenMP\. A single\#pragma omp parallelregion spawns a persistent thread team, using implicit barriers betweenttsteps to avoid repeated fork/join overhead\.

Pre\-computation phases exploit full independence across their iteration spaces\. TheBrick Construction, for instance, distributesN⋅D⋅NN\\cdot D\\cdot Nindependent element\-wise sums across all three axes viacollapse\(3\)\. Similarly, theCached Emission ComputingexposesD⋅ND\\cdot Nindependent work units viacollapse\(2\)\. The key parallelization challenge lies in the fusedBrick Update and Maximum Extraction, which in the single\-core version offers onlyNNindependent tasks\. We decompose it into two phases: Phase A distributes work over\(sj,d\)\(s\_\{j\},d\)pairs, where each thread sweeps over all source statessis\_\{i\}to find the local maximum, exposingN⋅DN\\cdot Dindependent tasks\. Phase B then reduces overddper destination statesjs\_\{j\}\. This decomposition enables full thread utilization even whenNNalone is smaller than the available core count\.

### IV\-CGPU Implementation

The tensor\-based Viterbi algorithm is inherently suited for GPU architectures; therefore, building upon our initial CPU version, we developed what is the first GPU\-accelerated Viterbi implementation for Hidden Semi\-Markov Models to our knowledge\. TheTens\-gpuimplementation is developed in CUDA and ported to HIP viahipify\[[24](https://arxiv.org/html/2609.16500#bib.bib27)\], maintaining a single codebase for both NVIDIA and AMD architectures\.

The primary computational bottleneck lies in constructing \(ℬf​i​r​s​t\\mathcal\{B\}\_\{first\}\) and updating \(ℬ\\mathcal\{B\}\) theBricktensor\. To maximize hardware utilization, we employ a fine\-grained work distribution where each individual thread is responsible for computing a singleBrickelement\. Then, after obtaining the finalBrickversion within each time\-step iteration, the threads operating on elements of the samesjs\_\{j\}\-slice must cooperate to perform the maximum extraction\.

To implement this mapping, we decompose theN×N×DN\\times N\\times Dtensor intoN×NN\\times Nvectors of sizeDDaligned along thezz\-axis, as illustrated in Figure[6](https://arxiv.org/html/2609.16500#S4.F6)\. This data layout is mapped onto a 2D grid ofN×NN\\times Nthread blocks, where each block manages a single vector\. TheDDelements are evenly divided between the threads in the block\. This configuration ensures that each thread block handles an entire temporal slice of the tensor regardless of the durationDD\.

Fig\. 6:Mapping theBrickto GPU Thread Blocks\.#### IV\-C1Memory Access Patterns and Coalescing

The emission cache \(orderedj⋅D\+dj\\cdot D\+d\), the brickℬf​i​r​s​t\\mathcal\{B\}\_\{first\}\(orderedj⋅N⋅D\+i⋅D\+dj\\cdot N\\cdot D\+i\\cdot D\+d\), and𝚫\\boldsymbol\{\\Delta\}\(indexedi⋅T\+\(t−1−d\)i\\cdot T\+\(t\-1\-d\)\) all maintaindd\-contiguity for consecutive and coalesced thread access\. Meanwhile, the emission probability𝐄\\mathbf\{E\}for the current observationoto\_\{t\}and statesjs\_\{j\}is shared by all threads within a block, allowing the GPU to serve it via a single broadcast read\.

#### IV\-C2Cached Emissions Computing

Naively computing𝐄\\mathbf\{E\}at each time\-step requiresO⁡\(D\)O\(D\)work per thread\. To eliminate redundant operations, we employ the same cached strategy explored in Sec\.[IV\-B2](https://arxiv.org/html/2609.16500#S4.SS2.SSS2)\. This reduces complexity toO⁡\(1\)O\(1\)via a double\-buffering scheme\. Thread\(j,i,d\)\(j,i,d\)reads𝐄\[:,d−1\]\\mathbf\{E\}\[:,d\-1\]from𝐄c​a​c​h​e\\mathbf\{E\}\_\{cache\}, addsbj​\(ot\)b\_\{j\}\(o\_\{t\}\), and writes the result to𝐄f​i​n​a​l\\mathbf\{E\}\_\{final\}\. The value is kept in a register for the current iteration, and the two buffers alternate roles across time\-steps\.

#### IV\-C3SeparatedBrickUpdate and Maximum Extraction

As shown in Sec\.[III](https://arxiv.org/html/2609.16500#S3), the tensor formulation decomposes each time\-step into two phases:Brick UpdateandMaximum Extraction\. These phases are fused in the CPU implementation but split across two GPU kernels to better exploit its architecture\. Note, the time\-invariant componentℬf​i​r​s​t\\mathcal\{B\}\_\{first\}is precomputed once before the induction loop\.

##### First Kernel

each thread\(sj,si,d\)\(s\_\{j\},s\_\{i\},d\)computes one element of the finalBrickℬ\\mathcal\{B\}and stores it directly in shared memory, bypassing global memory\. An intra\-block parallel reduction then identifies the local maximum score and its\(si,d\)\(s\_\{i\},d\)coordinates inO⁡\(log⁡D\)O\(\\log D\)steps: cross\-warp steps use\_\_syncthreads\(\), while the finallog2⁡\(warpsize\)\\log\_\{2\}\(\\text\{warpsize\}\)steps switch to warp\-shuffle instructions \(\_\_shfl\_down\_sync\), keeping the running maximum entirely in registers\.

##### Second Kernel

aggregates the per\-block local maximums to determine the globalarg⁡max\\arg\\maxover the full\(si,d\)\(s\_\{i\},d\)plane for eachsjs\_\{j\}\.

Splitting the computation into two kernels avoids grid\-wide synchronization, which would otherwise require cooperative\-groups designs with additional occupancy constraints\.

## VExperimental Results

To assess the performance of our tensor\-based implementations, we present an in\-depth experimental analysis\. We first describe the evaluation environment \(Sec\.[V\-A](https://arxiv.org/html/2609.16500#S5.SS1)\) and validate our implementations \(Sec\.[V\-B](https://arxiv.org/html/2609.16500#S5.SS2)\)\. Then, we evaluate the performance of both CPU \(Sec\.[V\-C](https://arxiv.org/html/2609.16500#S5.SS3)\) and GPU \(Sec\.[V\-D](https://arxiv.org/html/2609.16500#S5.SS4)\) implementations, including an investigation into the impact of varying hardware architectures \(Sec\.[V\-E](https://arxiv.org/html/2609.16500#S5.SS5)\)\. Last, we analyze the energy consumption \(Sec\.[V\-F](https://arxiv.org/html/2609.16500#S5.SS6)\) and conduct a stress test using extreme\-scale inputs \(Sec\.[V\-G](https://arxiv.org/html/2609.16500#S5.SS7)\)\.

### V\-AEvaluation Environment

#### V\-A1Baseline

As a baseline for validation and performance comparison, we selectedhsmmlearn\[[16](https://arxiv.org/html/2609.16500#bib.bib31)\], a C\+\+ library \(with a Python API\) for HSMMs with explicit duration distributions originating from the Rhsmmpackage\[[15](https://arxiv.org/html/2609.16500#bib.bib32)\], making it one of the most established HSMM codebases\. We chose it for two reasons\. First, it implements the general HSMM formulation with explicit, non\-parametric duration modeling and categorical emissions, matching the problem addressed in this work\. Second, among the HSMM frameworks analyzed in a recent survey\[[13](https://arxiv.org/html/2609.16500#bib.bib16)\], it is one of the few combining a general\-purpose formulation with a C\+\+ backend, ensuring our comparison targets optimized compiled code\. We refer to this single\-core baseline asBase\-1c\. We also considerededhsmm\[[25](https://arxiv.org/html/2609.16500#bib.bib34)\], which is inspired byhsmmlearn, but its Viterbi algorithm is implemented in Cython, potentially limiting low\-level compiler optimization compared to a native C\+\+ implementation\.

Ashsmmlearnis strictly limited to a single\-core implementation, we developed a multi\-core variant,Base\-mc, to ensure a fair comparison\. This was achieved by parallelizing the C\+\+ Viterbi decoder with OpenMP, specifically targeting the loop over states, the only loop in the classical four\-nested\-loops formulation with fully independent iterations\.Base\-mcrepresents the maximum parallelism extractable from the sequential formulation without significant algorithmic restructuring, serving as our primary baseline for both multi\-core and GPU comparisons\. It exhibits near\-linear scaling with thread count, provided the number of threads does not exceedNN\. Beyond that point, the speedup saturates and degrades slightly due to synchronization overhead, confirming that the traditional formulation fundamentally limits the exploitable parallelism to a single axis\. All implementations, including both the baselines and our tensor\-based versions, utilize double\-precision \(FP64\) floating\-point arithmetic\.

#### V\-A2Problem Size

To evaluate our formulation under realistic conditions, we selected problem sizes guided by computational genomics\. We tested configurations with a number of states \(NN\) ranging from1010to7575, spanning prokaryotic gene finders\[[26](https://arxiv.org/html/2609.16500#bib.bib19)\]at the lower end, chromatin state annotation\[[27](https://arxiv.org/html/2609.16500#bib.bib20)\]in the mid\-range \(1515–2525\), and eukaryotic gene finders such as AUGUSTUS\[[28](https://arxiv.org/html/2609.16500#bib.bib21)\]at the upper end\.

For the sequence length \(TT\), we tested from10310^\{3\}to10710^\{7\}time steps, covering typical gene\-finding invocations\[[29](https://arxiv.org/html/2609.16500#bib.bib22),[5](https://arxiv.org/html/2609.16500#bib.bib23)\]up toT=106T=10^\{6\}\[[11](https://arxiv.org/html/2609.16500#bib.bib24)\]\. For the maximum duration\(D\)\(D\), we adopted values from100100to10,00010,000, ranging from the explicit intron duration cutoff of SNAP\[[5](https://arxiv.org/html/2609.16500#bib.bib23)\]to stress\-test scenarios capturing the longest gene structure features in the human genome\[[30](https://arxiv.org/html/2609.16500#bib.bib25)\]\.

#### V\-A3Architectures

We evaluated our implementations on a representative set of high\-performance CPU and GPU architectures, summarized in Table[III](https://arxiv.org/html/2609.16500#S5.T3)\.

TABLE III:Hardware specifications\. SM: Streaming Multiprocessor; CU: Compute Unit; GCD: Graphics Compute Die\.

### V\-BValidation Results

All implementations produce identical output tohsmmlearnacross every tested configuration, achieving100%100\\%decoding accuracy\. To ensure exact equivalence, we include the tail adjustment phase without further optimization, matchinghsmmlearn’s boundary handling\. Notably,Tens\-Py, our direct NumPy transcription of the tensor formulation \(Algorithm[2](https://arxiv.org/html/2609.16500#alg2)\), already achieves a4\.5×4\.5\\timesspeedup overBase\-1c\. Despite comparing interpreted Python against compiled C\+\+, this result demonstrates that the tensor reformulation alone yields substantial gains before any low\-level optimization is applied\.

### V\-CCPU Speedup Analysis

We begin by comparingTens\-1cagainstBase\-1con Intel Xeon 8480\+\. Figure[7\(a\)](https://arxiv.org/html/2609.16500#S5.F7.sf1)reports the speedup forT=105T=10^\{5\}acrossN∈\{10,15,25,50,75\}N\\in\\\{10,15,25,50,75\\\}andD∈\{100,250,500,1000\}D\\in\\\{100,250,500,1000\\\}\.Tens\-1cis consistently faster, with speedups ranging from8\.7×8\.7\\times\(N=75N\\\!=\\\!75,D=1000D\\\!=\\\!1000\) to14\.1×14\.1\\times\(N=10N\\\!=\\\!10,D=100D\\\!=\\\!100\)\. The speedup decreases as eitherNNorDDgrows, reflecting the point at whichD×ND\\times Nexceeds L2 cache capacity\.

![Refer to caption](https://arxiv.org/html/2609.16500v1/tens-1c-vs-base-1c_xeon8480_intel_t100000.png)\(\(a\)\)Tens\-1cspeedup overBase\-1candTens\-1cruntime \(in parentheses; s: seconds; m: minutes\)\. Xeon 8480\+,T=105T=10^\{5\}\.
\(\(b\)\)Baselines and Tensors profiling metrics on Xeon 8480\+ CPU with ICX compiler\.

Fig\. 7:Speedup and runtime ofTens\-1cand profiling report\.We profiled the execution using LIKWID\[[31](https://arxiv.org/html/2609.16500#bib.bib28)\], which reveals that the tensor reformulation reduces retired instructions by36×36\\times\(from 610B to 16\.9B\)\. Figure[7\(b\)](https://arxiv.org/html/2609.16500#S5.F7.sf2)highlights the resulting hardware efficiency: the vectorization ratio increases from<0\.001%<0\.001\\%inBase\-1cto30\.9%30\.9\\%inTens\-1c, while branch misprediction overhead drops from6\.9%6\.9\\%to0\.46%0\.46\\%\. The figure also shows a reduced unique DRAM footprint, falling from2\.632\.63GB to1\.421\.42GB\. This, combined with hardware counter evidence of a shift from L2 reuse to streaming L3 access \(not shown in the figure\), explains both the massive absolute gains and their gradual erosion at largerNNandDD\. Nevertheless, even at the largest configuration, the speedup remains significant:Tens\-1ccompletes in12\.212\.2minutes, whereasBase\-1crequires over1\.761\.76hours\.

We now compareTens\-mcagainstBase\-mcon Xeon 8480\+\. Figure[8](https://arxiv.org/html/2609.16500#S5.F8)reports the speedup forT=104T=10^\{4\}andT=105T=10^\{5\}across the same value ofDD\(duration\) andNN\(number of states\) considered before\. The trend is clear:Tens\-mcbecomes increasingly effective asNNandDDgrow\. AtT=105T=10^\{5\}the speedup overBase\-mcreaches11\.4×11\.4\\times\(N=50N\\\!=\\\!50,D=1000D\\\!=\\\!1000\), while atT=104T=10^\{4\}a peak of11\.5×11\.5\\timesis observed at the same configuration\. This behavior is the direct consequence of the parallelization strategy\.Base\-mccan only distribute work over theNNdestination statessjs\_\{j\}, since the sequential formulation carries loop dependencies across durations and source states\.

![Refer to caption](https://arxiv.org/html/2609.16500v1/tens-mc-vs-base-mc_xeon8480_intel_t10k_t100k.png)Fig\. 8:Speedup ofTens\-mcoverBase\-mcon Xeon 8480\+\. The cell shows the speedup and the runtime ofTens\-mc\.In contrast,Tens\-mcreformulates the computation as broadcasted sums with no loop\-carried dependencies, exposingN⋅DN\\cdot Dindependent tasks\. This exposes sufficient parallelism to fully saturate the 112 available threads, even whenNNalone would leave most cores idle\. LIKWID profiling confirms this efficiency:Tens\-mcretires only 9\.9B instructions compared to 21\.9B forBase\-mc, while driving L3 bandwidth to near\-peak utilization of the memory hierarchy\. Furthermore, Figure[7\(b\)](https://arxiv.org/html/2609.16500#S5.F7.sf2)shows thatTens\-mcachieves an80\.4%80\.4\\%vectorization ratio \(versus0\.003%0\.003\\%forBase\-mc\), reduces bad speculation from6\.25%6\.25\\%to1\.16%1\.16\\%, and lowers DRAM volume from4\.594\.59GB to1\.431\.43GB\.

At small problem sizes, the trend reverses: forT=104T=10^\{4\}withN=10N=10andD=100D=100, theTens\-mcspeedup falls to0\.1×0\.1\\times\. In this regime, pre\-computation and barrier overhead dominate because the per\-step workload is insufficient to amortize them\. This deficit shrinks to0\.5×0\.5\\timesatT=105T\\\!=\\\!10^\{5\}and vanishes, becoming a14×14\\timesimprovement, atT=106T\\\!=\\\!10^\{6\}\(as we will show in Fig\.[10\(a\)](https://arxiv.org/html/2609.16500#S5.F10.sf1)\), confirming that the overhead is successfully amortized over the longer sequences typical of genomic applications\. In Sec\.[V\-E](https://arxiv.org/html/2609.16500#S5.SS5), we will evaluate bothTens\-mcandTens\-1cacross additional CPU architectures, demonstrating that they consistently outperform their respective baselines\.

Moreover, it is worth remarking that when compiled with GCC, theTens\-mcruntime forD=100D=100,N=10N=10, andT=104T=10^\{4\}drops to0\.340\.34seconds\. This highlights specific inefficiencies in the Intel compiler’s OpenMP implementation for small input sizes\. Nevertheless, we report all data using the Intel compiler as it consistently outperformed GCC across all other configurations\.

![Refer to caption](https://arxiv.org/html/2609.16500v1/tens-gpu-vs-base-mc_h100_cuda__xeon8480_intel_t10k_t100k.png)Fig\. 9:Speedup ofTens\-gpuon H100 overBase\-mcon Xeon 8480\+\. Cells showTens\-gpuspeedup and runtime\.![Refer to caption](https://arxiv.org/html/2609.16500v1/sc_1000000t_annotated.png)\(\(a\)\)Speedup ofTens\-1coverBase\-1c\.
\(\(b\)\)Speedup ofTens\-mcandTens\-gpuoverBase\-mcfor N=50\.
Fig\. 10:Speedup ofTens\-1c,Tens\-mc, andTens\-gpuover different CPUs and GPUs, for T=1,000,000\.
### V\-DGPU Speedup Analysis

Since no GPU implementation of the HSMM Viterbi exists in the literature we useBase\-mc, the multi\-core variant we developed by integrating OpenMP parallelization into the sequential C\+\+ backend ofhsmmlearn, as our reference\. We rely onBase\-mcas a baseline for two reasons: in addition to the total lack of reference GPU implementations, the loop\-carried dependencies in the standard sequential formulation across durations and states \(sis\_\{i\}\) fundamentally preclude a direct GPU port of the traditional Viterbi algorithm for HSMMs\. Consequently, the algorithmic restructuring proposed in this work is the necessary prerequisite for GPU acceleration\. Figure[9](https://arxiv.org/html/2609.16500#S5.F9)reports the resulting speedup ofTens\-gpu\(H100\) overBase\-mc\(Xeon 8480\+\) forT∈\{104,105\}T\\in\\\{10^\{4\},10^\{5\}\\\}\.

Tens\-gpuachieves speedups across the entire parameter space, peaking at36\.8×36\.8\\times\(N=25N\\\!=\\\!25,D=1000D\\\!=\\\!1000,T=105T\\\!=\\\!10^\{5\}\)\. Performance scales with bothNNandDD: largerNNincreases theN×NN\\times Nthread\-block grid, better saturating the H100’s 132 SMs, while largerDDprovides longer reduction vectors per block, improving warp\-shuffle efficiency\. Even at the smallest configuration \(N=10N\\\!=\\\!10,D=100D\\\!=\\\!100\), the speedup remains3\.63\.6–4\.0×4\.0\\times, despite only 100 thread blocks being insufficient to fully occupy all SMs\.

For largeNN\(N=75N\\\!=\\\!75\), the speedup plateaus \(e\.g\.,17\.7×17\.7\\timesatD=1000D\\\!=\\\!1000,T=105T\\\!=\\\!10^\{5\}\): the per\-block shared memory footprint grows withDD, and the second reduction kernel oversis\_\{i\}states becomes a bottleneck asNNincreases\. AtT=106T=10^\{6\}the advantage grows further, reaching54\.3×54\.3\\timesatN=15N\\\!=\\\!15,D=500D\\\!=\\\!500:Base\-mcrequires over5\.75\.7minutes on 112 CPU cores, whileTens\-gpucompletes the same decoding in6\.46\.4seconds\. Notably, the absolute GPU runtimes remain sub\-second for most configurations and never exceed1010seconds even at the largest tested \(N=75N\\\!=\\\!75,D=1000D\\\!=\\\!1000,T=105T\\\!=\\\!10^\{5\}\), making interactive\-scale HSMM decoding on large inputs feasible for the first time\. These runtimes open the door to problem sizes that were previously intractable, as we will explore in Sec\.[V\-G](https://arxiv.org/html/2609.16500#S5.SS7)\.

### V\-EArchitecture Comparison

Figure[10\(a\)](https://arxiv.org/html/2609.16500#S5.F10.sf1)comparesTens\-1cagainstBase\-1catT=106T=\\\!10^\{6\}withN∈\{10,15,25\}N\\in\\\{10,15,25\\\}on the three CPUs described in Table[III](https://arxiv.org/html/2609.16500#S5.T3)\.Tens\-1cdelivers consistent speedups of1010–11×11\\timeson Grace,1111–14×14\\timeson Xeon, and1111–12×12\\timeson EPYC across all analyzed configurations, confirming that the gains are portable across architectures\.

In absolute terms, forN=25N\\\!=\\\!25andD=1000D\\\!=\\\!1000,Base\-1crequires2\.32\.3hours on Xeon, whichTens\-1creduces to12\.412\.4minutes\. ForN=50N\\\!=\\\!50andN=75N\\\!=\\\!75\(not shown\), speedups remain consistent\. AtN=75N\\\!=\\\!75andD=500D\\\!=\\\!500,Tens\-1creduces the runtime from9\.29\.2hours to5656minutes on Xeon; forD=1000D\\\!=\\\!1000,Base\-1ctimed out \(exceeding1515hours\), whereasTens\-1ccompleted in22hours\. Grace exhibits the lowest runtime, due to its \(almost2×2\\times\) higher memory bandwidth\.

Figure[10\(b\)](https://arxiv.org/html/2609.16500#S5.F10.sf2)extends the comparison toTens\-mcandTens\-gpuatT=106T\\\!=\\\!10^\{6\}withN=50N\\\!=\\\!50across all the CPUs and GPUs introduced in Table[III](https://arxiv.org/html/2609.16500#S5.T3)\. ForTens\-mc, speedups are reported overBase\-mcon the same CPU\. ForTens\-gpu, speedups are reported over the CPU whereBase\-mcis fastest \(Grace\)\.

On the CPU side, atD=1000D\\\!=\\\!1000Tens\-mcshows a speedup overBase\-mcof2×2\\timeson EPYC and8×8\\timeson Xeon\. On the GPU side, the speedup overBase\-mcgrows steadily withDD, as increasing the duration expands the per\-block workload and improves SM occupancy\. AtD=1000D\\\!=\\\!1000, the speedups overBase\-mcrange from5×5\\timeson MI250X to12×12\\timeson H200\. The lower performance for MI250X can be attributed to its lower memory bandwidth and fewer compute units per GCD\.

Across all GPUs, the absolute runtimes atD=500D\\\!=\\\!500remain below3030seconds for a million\-step sequence\. For comparison, the original unmodifiedhsmmlearnbaselineBase\-1crequires3\.83\.8hours on Grace for this configuration;Tens\-gpuon a H200 completes the same decoding in2424seconds, a reduction of570×570\\times\.

It is worth mentioning that, to maximize the breadth of our architectural comparison, we opted for a single CUDA/HIP GPU codebase, and OpenMP\-only multi\-core parallelization\. Further specialization is possible on both fronts, and the GPU kernels could exploit architecture\-specific features such as distinct memory hierarchies or generation\-specific instructions\.

### V\-FEnergy Consumption

Figure[11\(a\)](https://arxiv.org/html/2609.16500#S5.F11.sf1)reports the energy consumption of all implementations atN=50N\\\!=\\\!50,T=104T\\\!=\\\!10^\{4\}, normalized toBase\-1c\. For the sake of space, and because we observed a similar trend for the other values ofDD, we only report the data forD∈\{100,1000\}D\\in\\\{100,1000\\\}\. For the CPU versions, we measure the energy on the EPYC 7A53, whereas forTens\-gpuwe measure the energy on the MI250X\. In both cases, energy is monitored through the Cray Power Management \(PM\) counters\[[32](https://arxiv.org/html/2609.16500#bib.bib39)\]\. The dominant factor is execution time: since the instantaneous power draw remains comparable across CPU implementations, energy reductions closely track runtime reductions\.

\(\(a\)\)Energy consumption overBase\-1cforN=50N=50,T=104T=10^\{4\}on EPYC 7A53 and MI250X\.\(\(b\)\)Tens\-gpuperformance on a stress test case withT=107T=10^\{7\},N=100N=100,D=104D=10^\{4\}\.
Fig\. 11:Energy consumption and stress test\.Tens\-1creduces energy consumption by∼10×\{\\sim\}10\\timesoverBase\-1cacross all testedDDvalues, consistent with its single\-core speedup\. In the multi\-core regime,Tens\-mcconsumes525\.5525\.5J atD=1000D\\\!=\\\!1000, a3×3\\timesreduction compared toBase\-mc\(1\.51\.5kJ\), demonstrating that the tensor reformulation translates its runtime advantage into proportional energy savings\.Tens\-gpuachieves the lowest energy footprint, requiring only327\.7327\.7J atD=1000D\\\!=\\\!1000and85\.285\.2J atD=100D\\\!=\\\!100\(just2%2\\%ofBase\-1c\)\. Although the GPU exhibits higher instantaneous power draw, its shorter execution time more than compensates, rendering it the most energy\-efficient platform forD\>100D\\\!\>\\\!100\. These results confirm that the tensor formulation not only accelerates HSMM Viterbi decoding but also enables a significantly more energy\-efficient profile\.

### V\-GStress Test: Beyond Current Workloads

To demonstrate the practical impact of our formulation, we evaluateTens\-gpuon an extreme\-scale configuration:N=100N\\\!=\\\!100states,D=10,000D\\\!=\\\!10\{,\}000maximum duration, andT=107T\\\!=\\\!\\\!10^\{7\}time steps\. This scale is entirely inaccessible to the baseline;Base\-1ctriggers a memory allocation failure \(std::bad\_alloc\) for configurations exceedingT=106T\\\!=\\\!10^\{6\}andD=1,000D\\\!=\\\!1\{,\}000, precluding direct measurement\. By extrapolating from runtimes measured atN=75N\\\!=\\\!75,D=100D\\\!=\\\!100,T=106T\\\!=\\\!10^\{6\}\(whereBase\-1crequires22h andBase\-mc6\.76\.7min\) using theO⁡\(T⋅N2⋅D\)O\(T\\cdot N^\{2\}\\cdot D\)theoretical complexity, we estimate this extreme configuration would require approximately148148days forBase\-1cand2\.12\.1days forBase\-mcon 112 cores\. Such runtimes render not only individual decoding tasks impractical but also make Viterbi training, which requires dozens of such iterations, entirely infeasible on traditional architectures\.

Figure[11\(b\)](https://arxiv.org/html/2609.16500#S5.F11.sf2)reports the per\-iteration runtime ofTens\-gpuon five GPUs\. The H200 leads at53\.753\.7minutes, followed by the MI300X at57\.357\.3minutes, the H100 at1\.71\.7hours, the A100 at2\.02\.0hours, and the MI250X at3\.43\.4hours\. The H200 and MI300X’s advantage over the H100 is consistent with their higher memory bandwidth and larger number of compute units: at this scale, theN×N=10,000N\\times N=10\{,\}000thread\-block grid fully saturates both architectures, and performance becomes bandwidth\-bound, favoring the MI300X and H200\. The A100 trails the H100 due to its lower bandwidth and fewer SMs\. The MI250X, despite its110110CUs per GCD, is bottlenecked by its HBM2e bandwidth, the lowest among the five\.

These results demonstrate thatTens\-gpureduces a previously intractable workload, estimated to take over a month on a single core, to less than an hour on a single GPU, making whole\-genome\-scale HSMM decoding and iterative Viterbi training practically feasible even for larger sequences\.

## VIRelated Work

The Viterbi algorithm has been fundamental in high\-performance bioinformatics and signal processing for decades, yet existing acceleration efforts are almost exclusively devoted to standard HMMs\. Moving from HMMs to HSMMs introduces explicit state\-duration handling that substantially increases computational complexity: the Viterbi iteration must compute a maximum over all candidate durations, making the inner loops data\-dependent and inherently difficult to parallelize\. This combination of computational burden and parallelization difficulty helps explain why performant, hardware\-aware HSMM decoders remain absent from the literature\.

### VI\-AHidden Markov Models \(HMMs\)

HSMMs generalize standard HMMs by introducing explicit state\-duration distributions, raising computational complexity fromO⁡\(T​N2\)O\(TN^\{2\}\)toO⁡\(T​N2​D\)O\(TN^\{2\}D\)\. The Viterbi algorithm has been extensively accelerated for the simpler HMM formulation, including SIMD\-vectorized CPU frameworks\[[18](https://arxiv.org/html/2609.16500#bib.bib7),[20](https://arxiv.org/html/2609.16500#bib.bib38)\], CUDA\-based GPU implementations\[[33](https://arxiv.org/html/2609.16500#bib.bib9),[19](https://arxiv.org/html/2609.16500#bib.bib10),[34](https://arxiv.org/html/2609.16500#bib.bib11),[35](https://arxiv.org/html/2609.16500#bib.bib8)\], hardware–software co\-design and domain\-specific approaches\[[17](https://arxiv.org/html/2609.16500#bib.bib12),[36](https://arxiv.org/html/2609.16500#bib.bib14)\], and distributed computing\[[37](https://arxiv.org/html/2609.16500#bib.bib13)\]\. However, the additional duration dimension cannot simply be wrapped around existing HMM accelerators: it introduces a cumulative emission product over theddmost recent observations, requires accessing a variable\-depth window of past delta values, and turns the per\-state maximum into a joint maximization over both states and durations\. As a consequence, none of these efforts extend to HSMMs, and the HSMM formulation remains entirely unaddressed\.

### VI\-BHidden Semi\-Markov Models \(HSMMs\)

Several statistical frameworks implement major HSMM algorithms in R or Python, with performance\-critical routines in C/C\+\+\[[13](https://arxiv.org/html/2609.16500#bib.bib16)\]\. Domain\-specific solutions also exist, such as biomvRhsmm\[[14](https://arxiv.org/html/2609.16500#bib.bib17)\]for genomic segmentation, and Pertsinidou and Limnios\[[38](https://arxiv.org/html/2609.16500#bib.bib18)\]that propose Viterbi algorithms based on the backward recurrence Markov chain formulation\. All of these implementations, however, are sequential and single\-threaded, and none explicitly targets modern high\-performance CPUs or GPUs\. The most closely related work is Lu et al\.\[[39](https://arxiv.org/html/2609.16500#bib.bib15)\], who propose a Tensor\-based HSMM \(T\-HSMM\) for user activity analysis in Cyber\-Physical\-Social Systems \(CPSSs\)\. Their objective differs fundamentally from ours: their*tensor*refers to embedding multiple correlated entities in a unified higher\-dimensional space for activity modeling, rather than targeting computational acceleration, whereas ours reshapes the three inner loops of the Viterbi inductive phase into several 3D tensor operations that expose parallelism for high\-performance CPU and GPU execution\. Moreover, Lu et al\. collapse the duration dimension into a single scalar expected value per state, which alters the HSMM semantics and does not solve the exact Viterbi decoding problem\. Our formulation instead preserves the full duration dimension and performs exact decoding\. A direct head\-to\-head performance comparison is therefore not meaningful, since the two methods solve different problems\.

### VI\-CSummary

To the best of our knowledge, no prior work presents a high\-performance implementation of the Viterbi algorithm for HSMMs\. Our work fills this gap by proposing a tensor\-based reformulation of the HSMM Viterbi algorithm, opening the way to new optimization strategies, accelerator implementations, and application\-specific mappings for domains that require Hidden Semi\-Markov Model modeling\.

## VIIDiscussion

We now discuss the main design choices behind our formulation, the trade\-offs they involve, and the technical directions they leave open\.

### VII\-AAlternative Algorithmic Formulations

A lower\-complexity formulation is in principle available by factorizing the induction, reducing over source states before combining the duration and emission terms, which lowers the per\-step cost fromO⁡\(N2​D\)O\(N^\{2\}D\)toO⁡\(N​D\+N2\)O\(ND\+N^\{2\}\)\. We do not adopt it because the saving in arithmetic is offset by a loss of hardware efficiency\. On CPU, the factorization removes the time\-invariant*Brick*precomputation and with it the contiguous\(j,d,i\)\(j,d,i\)layout that keeps the working set resident across the duration sweep\. On GPU, it splits a single joint maximization over the\(si,d\)\(s\_\{i\},d\)plane into two reductions that must run one after the other, adding a second grid\-wide synchronization per time step and reducing occupancy, which is exactly the pattern our two\-kernel design avoids\. Evaluating this factorization under a different tensor formulation, built around its own data layout and reduction scheme, would nonetheless be an interesting direction\.

### VII\-BMapping onto Specialized Accelerators

Since emerging AI accelerators and dataflow architectures are designed precisely to execute dense tensor operations, mapping our formulation onto tensor cores, TPUs, systolic arrays, or FPGA dataflow designs is a natural direction to consider\. The obstacle is the kind of reduction involved\. Since all quantities are handled in log\-space, each step combines values with an addition and then selects a maximum\. These accelerators are instead built around multiply\-accumulate pipelines, so a maximum\-based reduction does not map directly onto their native primitives and would need a dedicated mapping strategy\. Studying how to support such operations on this class of hardware would therefore be valuable well beyond our setting, since it would open these units to dynamic programming algorithms in general\. A further consideration is that we use double precision to match the baseline exactly, whereas peak throughput on these units is available only at lower precision, so any port must first verify that the dynamic range of the problem allows a narrower format\.

### VII\-CSequence\-Level and Distributed Parallelism

Our implementations decode a single sequence at a time, from start to end\. A natural extension is to split a long sequence into chunks, decode them in parallel, and then reconcile the results at the chunk boundaries\. This would also enable multi\-node execution, where each node handles a portion of the sequence and the boundary values are exchanged through collective operations\. Decoding several independent sequences at once is another promising direction, since it would keep the device busy on small inputs, where a single sequence leaves many units idle\.

### VII\-DModel Assumptions and Algorithmic Scope

Our formulation uses one global maximum durationDDfor all states, which keeps the tensor dense and the work per thread uniform\. Giving each state its own boundDjD\_\{j\}would avoid computing durations that a state can never take, at the cost of an irregular*Brick*\. We also assume discrete emissions, which makes the emission term a simple table lookup and enables our cached update; continuous densities such as Gaussians would require a different caching strategy\. Finally, the Forward\-Backward and Baum\-Welch procedures iterate over the same\(sj,si,d\)\(s\_\{j\},s\_\{i\},d\)combinations and only replace the maximum with a sum, so they can reuse the same tensor operations and be accelerated in the same way\.

## VIIIConclusions

We presented a tensor\-based formulation of the Viterbi algorithm for Hidden Semi\-Markov Models that restructures the three inner loops of the sequential algorithm into dense tensor operations, exposing optimization opportunities that are inaccessible to the traditional scalar formulation\. Building on it, we delivered optimized single\-core CPU, multi\-core CPU, and, for the first time for HSMMs, GPU implementations, released as the open\-source librarytensor\-hsmm\.111https://github\.com/HLC\-Lab/tensor\-hsmm/

Across three CPU and five GPU architectures, our implementations achieve speedups of up to14×14\\timeson a single core, over200×200\\timeswith multi\-core, and over570×570\\timeson GPU with respect to the sequentialBase\-1cbaseline, while producing output identical tohsmmlearnin every tested configuration\. The gains are structural rather than platform\-specific: a direct NumPy transcription of the formulation already outperforms the compiled sequential baseline by4\.5×4\.5\\times, and profiling attributes the compiled speedups to a36×36\\timesreduction in retired instructions together with vectorization ratios rising from below0\.001%0\.001\\%to80\.4%80\.4\\%\. Because instantaneous power draw is comparable across implementations, these runtime reductions translate into proportional energy savings, with the GPU version consuming as little as2%2\\%of the baseline energy\. Most consequentially, a configuration estimated to require over a month of single\-core execution completes in under an hour on a single GPU, bringing whole\-genome\-scale HSMM decoding and iterative Viterbi training within practical reach and establishing a new performance baseline for large\-scale HSMM inference\.

## Acknowledgements

We acknowledge ISCRA for awarding this project access to the LEONARDO supercomputer, owned by the EuroHPC Joint Undertaking, hosted by CINECA \(Italy\)\. We acknowledge the EuroHPC Joint Undertaking, the LUMI consortium, and BSC for granting access to the LUMI and MareNostrum 5 supercomputers\. These resources, hosted by CSC \(Finland\) and the Barcelona Supercomputing Center \(Spain\), were provided through the EuroHPC Regular Access program\. The authors used Claude Opus 4\.6 and Gemini 3 for editing the paper; all ideas, content, and conclusions are their own\.

## References

- \[1\]L\. Gabriel, T\. Brůna, K\. J\. Hoff, M\. Ebel, A\. Lomsadze, M\. Borodovsky, and M\. Stanke\(2024\)BRAKER3: fully automated genome annotation using rna\-seq and protein evidence with genemark\-etp, augustus, and tsebra\.Genome research34\(5\),pp\. 769–777\.External Links:[Document](https://dx.doi.org/10.1101/gr.278090.123),[Link](https://doi.org/10.1101/gr.278090.123)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p2.1),[§I](https://arxiv.org/html/2609.16500#S1.p5.1)\.
- \[2\]S\. Qin, Z\. Tan, and Y\. Wu\(2024\)On robust estimation of hidden semi\-Markov regime\-switching models\.Annals of Operations Research\.External Links:[Document](https://dx.doi.org/10.1007/s10479-024-05989-4),[Link](https://doi.org/10.1007/s10479-024-05989-4)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p2.1)\.
- \[3\]H\. Zen, K\. Tokuda, T\. Masuuko, T\. Kobayasih, and T\. Kitamura\(2007\)A hidden semi\-markov model\-based speech synthesis system\.IEICE TRANSACTIONS on InformationE90\-D\(5\),pp\. 825–834\.External Links:[Document](https://dx.doi.org/10.1093/ietisy/e90-d.5.825),[Link](https://doi.org/10.1093/ietisy/e90-d.5.825),ISSN 1745\-1361Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p2.1)\.
- \[4\]S\. Yu\(2010\)Hidden semi\-markov models\.Artificial Intelligence174\(2\),pp\. 215–243\.Note:Special Review IssueExternal Links:ISSN 0004\-3702,[Document](https://dx.doi.org/https%3A//doi.org/10.1016/j.artint.2009.11.011),[Link](https://doi.org/10.1016/j.artint.2009.11.011)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p2.1),[§I](https://arxiv.org/html/2609.16500#S1.p4.1),[§I](https://arxiv.org/html/2609.16500#S1.p6.1),[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§II](https://arxiv.org/html/2609.16500#S2.p1.1)\.
- \[5\]I\. Korf\(2004\)Gene finding in novel genomes\.BMC bioinformatics5\(1\),pp\. 59\.External Links:[Link](https://doi.org/10.1186/1471-2105-5-59)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p5.1),[§I](https://arxiv.org/html/2609.16500#S1.p7.1),[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p2.1)\.
- \[6\]J\. Ernst and M\. Kellis\(2012\)ChromHMM: automating chromatin\-state discovery and characterization\.Nature Methods9\(3\),pp\. 215–216\.External Links:[Link](https://doi.org/10.1038/nmeth.1906),[Document](https://dx.doi.org/10.1038/nmeth.1906)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p5.1)\.
- \[7\]R\. Durbin, S\. R\. Eddy, A\. Krogh, and G\. Mitchison\(1998\)Biological sequence analysis: probabilistic models of proteins and nucleic acids\.Cambridge University Press\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1017/CBO9780511790492),[Link](https://doi.org/10.1017/CBO9780511790492)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p5.1)\.
- \[8\]L\.R\. Rabiner\(1989\)A tutorial on hidden markov models and selected applications in speech recognition\.Proceedings of the IEEE77\(2\),pp\. 257–286\.External Links:[Document](https://dx.doi.org/10.1109/5.18626),[Link](https://doi.org/10.1109/5.18626)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p6.1),[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§II](https://arxiv.org/html/2609.16500#S2.p1.1),[§IV](https://arxiv.org/html/2609.16500#S4.p3.1)\.
- \[9\]L\. E\. Baum, T\. Petrie, G\. Soules, and N\. Weiss\(1970\)A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains\.The Annals of Mathematical Statistics41\(1\),pp\. 164–171\.External Links:ISSN 00034851, 21688990,[Link](http://www.jstor.org/stable/2239727)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p6.1)\.
- \[10\]A\. P\. Dempster, N\. M\. Laird, and D\. B\. Rubin\(1977\)Maximum likelihood from incomplete data via the em algorithm\.Journal of the Royal Statistical Society: Series B \(Methodological\)39\(1\),pp\. 1–22\.External Links:[Document](https://dx.doi.org/https%3A//doi.org/10.1111/j.2517-6161.1977.tb01600.x),[Link](https://doi.org/10.1111/j.2517-6161.1977.tb01600.x)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p6.1)\.
- \[11\]A\. Lomsadze, V\. Ter\-Hovhannisyan, Y\. O\. Chernoff, and M\. Borodovsky\(2005\)Gene identification in novel eukaryotic genomes by self\-training algorithm\.Nucleic acids research33\(20\),pp\. 6494–6506\.External Links:[Link](https://doi.org/10.1093/nar/gki937)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p6.1),[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p2.1)\.
- \[12\]Y\. Guédon\(2003\)Estimating hidden semi\-markov chains from discrete sequences\.Journal of Computational and Graphical Statistics12\(3\),pp\. 604–639\.External Links:[Link](https://doi.org/10.1198/1061860032030),[Document](https://dx.doi.org/10.1198/1061860032030)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§II](https://arxiv.org/html/2609.16500#S2.p1.1)\.
- \[13\]C\. Bérard, M\. CROS, J\. DURAND, C\. Lothodé, S\. Plancade, R\. Trepos, and N\. Vergne\(2025\)Review of hsmm r and python softwares\.A Comprehensive Guide to HSMM: Theory, Software, and Advanced Extensions,pp\. 47–77\.External Links:[Link](https://doi.org/10.1002/9781394427581.ch2)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§IV](https://arxiv.org/html/2609.16500#S4.p3.1),[§V\-A1](https://arxiv.org/html/2609.16500#S5.SS1.SSS1.p1.1),[§VI\-B](https://arxiv.org/html/2609.16500#S6.SS2.p1.1)\.
- \[14\]Y\. Du, E\. Murani, S\. Ponsuksili, and K\. Wimmers\(2014\)BiomvRhsmm: genomic segmentation with hidden semi\-markov model\.BioMed Research International2014\(1\),pp\. 910390\.External Links:[Link](https://doi.org/10.1155/2014/910390)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§VI\-B](https://arxiv.org/html/2609.16500#S6.SS2.p1.1)\.
- \[15\]J\. Bulla and I\. Bulla\(2013\)Hsmm: hidden semi Markov models\.Note:R package, archived May 2022\.External Links:[Link](https://cran.r-project.org/package=hsmm)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§V\-A1](https://arxiv.org/html/2609.16500#S5.SS1.SSS1.p1.1)\.
- \[16\]Hsmmlearn: a library for hidden semi\-Markov models with explicit durationsNote:C\+\+/Cython with Python interface\. Wraps the C\+\+ code from the R hsmm package\. Archived January 2023\.External Links:[Link](https://github.com/jvkersch/hsmmlearn)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§V\-A1](https://arxiv.org/html/2609.16500#S5.SS1.SSS1.p1.1)\.
- \[17\]C\. Firtina, K\. Pillai, G\. S\. Kalsi, B\. Suresh, D\. S\. Cali, J\. S\. Kim, T\. Shahroodi, M\. B\. Cavlak, J\. Lindegger, M\. Alser,et al\.\(2024\)Aphmm: accelerating profile hidden markov models for fast and energy\-efficient genome analysis\.ACM Transactions on Architecture and Code Optimization21\(1\),pp\. 1–29\.External Links:[Link](https://doi.org/10.1145/3632950)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[18\]S\. R\. Eddy\(2011\)Accelerated profile hmm searches\.PLoS computational biology7\(10\),pp\. e1002195\.External Links:[Link](https://doi.org/10.1371/journal.pcbi.1002195)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[19\]L\. Yu, Y\. Ukidave, and D\. Kaeli\(2014\)GPU\-accelerated hmm for speech recognition\.In2014 43rd International Conference on Parallel Processing Workshops,pp\. 395–402\.External Links:[Link](https://doi.org/10.1109/ICPPW.2014.59)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[20\]H\. Jiang, N\. Ganesan, and Y\. Yao\(2018\)CUDAMPF\+\+: a proactive resource exhaustion scheme for accelerating homologous sequence search on cuda\-enabled gpu\.IEEE Transactions on Parallel and Distributed Systems29\(10\),pp\. 2206–2222\.External Links:ISSN 2161\-9883,[Link](http://dx.doi.org/10.1109/TPDS.2018.2830393)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1),[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[21\]S\. Hassan, S\. Sarkka, and A\. Garcia\-Fernandez\(2021\)Temporal parallelization of inference in hidden markov models\.IEEE Transactions on Signal Processing69,pp\. 4875–4887\.External Links:ISSN 1941\-0476,[Link](http://dx.doi.org/10.1109/TSP.2021.3103338),[Document](https://dx.doi.org/10.1109/tsp.2021.3103338)Cited by:[§I](https://arxiv.org/html/2609.16500#S1.p8.1)\.
- \[22\]L\. S\. Blackford, A\. Petitet, R\. Pozo, K\. Remington, R\. C\. Whaley, J\. Demmel, J\. Dongarra, I\. Duff, S\. Hammarling, G\. Henry,et al\.\(2002\)An updated set of basic linear algebra subprograms \(BLAS\)\.ACM Transactions on Mathematical Software28\(2\),pp\. 135–151\.External Links:[Link](https://doi.org/10.1145/567806.567807)Cited by:[§IV\-B](https://arxiv.org/html/2609.16500#S4.SS2.p1.1)\.
- \[23\]xtensor\-stack\(2025\)Xtensor: Multi\-dimensional arrays with broadcasting and lazy computing\.Note:C\+\+14 header\-only library\.https://github\.com/xtensor\-stack/xtensorCited by:[§IV\-B](https://arxiv.org/html/2609.16500#S4.SS2.p1.1)\.
- \[24\]AMD\(2024\)HIPIFY: convert CUDA to portable C\+\+ code\.Note:Accessed: 2026\-04\-06External Links:[Link](https://github.com/ROCm/HIPIFY)Cited by:[§IV\-C](https://arxiv.org/html/2609.16500#S4.SS3.p1.1)\.
- \[25\]Edhsmm: an\(other\) implementation of explicit duration hidden semi\-Markov models in Python 3Note:Python/Cython\. Archived May 2024\.External Links:[Link](https://github.com/poypoyan/edhsmm)Cited by:[§V\-A1](https://arxiv.org/html/2609.16500#S5.SS1.SSS1.p1.1)\.
- \[26\]A\. V\. Lukashin and M\. Borodovsky\(1998\)GeneMark\. hmm: new solutions for gene finding\.Nucleic Acids Research26\(4\),pp\. 1107–1115\.External Links:[Link](https://doi.org/10.1093/nar/26.4.1107)Cited by:[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p1.1)\.
- \[27\]A\. Kundaje, W\. Meuleman, J\. Ernst, M\. Bilenky, A\. Yen, P\. Kheradpour, Z\. Zhang, A\. Heravi\-Moussavi, Y\. Liu, V\. Amin,et al\.\(2015\)Integrative analysis of 111 reference human epigenomes\.Nature518\(7539\),pp\. 317\.External Links:[Link](https://doi.org/10.1038/nature14248)Cited by:[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p1.1)\.
- \[28\]M\. Stanke\(2003\)Gene prediction with a hidden markov model and a new intron submodel\.Bioinformatics\.External Links:[Link](https://doi.org/10.1093/bioinformatics/btg1080)Cited by:[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p1.1)\.
- \[29\]C\. Burge and S\. Karlin\(1997\)Prediction of complete gene structures in human genomic dna\.Journal of molecular biology268\(1\),pp\. 78–94\.External Links:[Link](https://doi.org/10.1006/jmbi.1997.0951)Cited by:[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p2.1)\.
- \[30\]M\. K\. Sakharkar, V\. T\. Chow, and P\. Kangueane\(2004\)Distributions of exons and introns in the human genome\.In silico biology4\(4\),pp\. 387–393\.External Links:[Link](https://doi.org/10.3233/ISB-00142)Cited by:[§V\-A2](https://arxiv.org/html/2609.16500#S5.SS1.SSS2.p2.1)\.
- \[31\]J\. Treibig, G\. Hager, and G\. Wellein\(2010\)LIKWID: a lightweight performance\-oriented tool suite for x86 multicore environments\.InProceedings of PSTI2010, the First International Workshop on Parallel Software Tools and Tool Infrastructures,pp\. 207–216\.External Links:[Link](https://doi.org/10.1109/ICPPW.2010.38),[Document](https://dx.doi.org/10.1109/ICPPW.2010.38)Cited by:[§V\-C](https://arxiv.org/html/2609.16500#S5.SS3.p2.1)\.
- \[32\]HPE Cray\(2026\)Cray Performance and Analysis Tools cray\_pm\.Note:Accessed: 2026\-04\-06External Links:[Link](https://cpe.ext.hpe.com/docs/24.03/performance-tools/index.html)Cited by:[§V\-F](https://arxiv.org/html/2609.16500#S5.SS6.p1.1)\.
- \[33\]M\. HoseinyFarahabady and A\. Y\. Zomaya\(2025\)GPU\-accelerated out\-of\-core hmm inference with concurrent cuda streams\.InInternational Conference on Computational Science,pp\. 369–376\.External Links:[Link](https://doi.org/10.1007/978-3-031-97635-3_44)Cited by:[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[34\]A\. Mohammadidoost and M\. Hashemi\(2020\)High\-throughput and memory\-efficient parallel viterbi decoder for convolutional codes on gpu\.arXiv preprint arXiv:2011\.09337\.External Links:[Link](https://doi.org/10.48550/arXiv.2011.09337)Cited by:[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[35\]V\. Roubtsova\(2023\)Parallel algorithm for a hidden markov model with an indefinite number of states and heterogeneous observation data\.\.InIWOCL,pp\. 31–1\.External Links:[Link](https://doi.org/10.1145/3585341.3587954)Cited by:[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[36\]L\. Hummelgren, V\. Palmkvist, L\. Stjerna, X\. Xu, J\. Jaldén, and D\. Broman\(2024\)Trellis: a domain\-specific language for hidden markov models with sparse transitions\.InProceedings of the 17th ACM SIGPLAN International Conference on Software Language Engineering,pp\. 196–209\.External Links:[Link](https://doi.org/10.1145/3687997.3695641)Cited by:[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[37\]I\. Sassi, S\. Anter, and A\. Bekkhoucha\(2021\)ParaDist\-hmm: a parallel distributed implementation of hidden markov model for big data analytics using spark\.International Journal of Advanced Computer Science and Applications12\(4\)\.External Links:[Document](https://dx.doi.org/10.14569/IJACSA.2021.0120438),[Link](http://dx.doi.org/10.14569/IJACSA.2021.0120438)Cited by:[§VI\-A](https://arxiv.org/html/2609.16500#S6.SS1.p1.1)\.
- \[38\]C\. Pertsinidou and N\. Limnios\(2015\)Viterbi algorithms for hidden semi\-markov models with application to dna analysis\.RAIRO\-Operations Research49\(3\),pp\. 511–526\.External Links:[Link](https://doi.org/10.1051/ro/2014053)Cited by:[§VI\-B](https://arxiv.org/html/2609.16500#S6.SS2.p1.1)\.
- \[39\]Z\. Lu, L\. T\. Yang, A\. Azman, F\. Zhou, S\. Zhang, and X\. Fu\(2025\)Tensor\-based hidden semi\-markov model for cpss user activity analysis and services\.IEEE Transactions on Services Computing\.External Links:[Link](https://doi.org/10.1109/TSC.2025.3618011)Cited by:[§VI\-B](https://arxiv.org/html/2609.16500#S6.SS2.p1.1)\.