arXiv CSopen access

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

arXiv:2609.16500v1 [cs.LG] 15 Sep 2026

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

Elia Belli

Daniele De Sensi

Department of Computer Science Sapienza University of Rome [email protected]

Department of Computer Science Sapienza University of Rome [email protected]

Department of Computer Science Sapienza University of Rome [email protected]

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 largescale 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 multicore CPUs, and, for the first time, GPU. Experimental evaluation demonstrates speedups of up to 14× on a single core, over 200× with multi-core, and over 570× on 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

I. I NTRODUCTION 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.

SC26, November 15-20, 2026, Chicago, Illinois, USA 979-8-3195-4789-7/26/$31.00 ©2026 IEEE

Hidden Semi-Markov Models (HSMMs) are versatile probabilistic frameworks with applications spanning diverse fields, from genome annotation and segmentation in computational biology [1] to finance [2], speech recognition [3], and signal processing [4]. As a generalization of the classical Hidden 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 sleepcycle 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. HSMMs generalize this framework by permitting each hidden state to persist for a variable interval, explicitly defined by a duration distribution [4]. 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], [5]. 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, CpG island detection, and other problems where segment lengths carry vital biological meaning [6], [7]. 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 [4], [8] computes state probabilities at each time step, the Baum–Welch algorithm [4], [8], [9], an instance of Expectation–Maximization [10], 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] 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 of O(T N 2 ), the HSMM formulation introduces an additional loop over all possible durations, raising the complexity to O(T N 2 D), where D is the maximum admissible duration. In practical applications where D ranges from hundreds to thousands [5], 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], [8], [12]. Despite the proliferation of HPC resources, existing implementations remain confined to scalar, single-threaded CPU execution [13]–[16], in stark contrast to standard HMMs, which benefit from an extensive ecosystem of GPU-accelerated and SIMD-optimized frameworks [17]–[21]. 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) We introduce a novel Tensor-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) 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) We conduct an extensive performance evaluation across three CPU and five GPU architectures, demonstrating speedups of up to 570× over state-of-the-art HSMM frameworks. II. S EQUENTIAL HSMM V ITERBI A LGORITHM We now formalize the elements that define a Hidden SemiMarkov Model. An HSMM is represented by the parameter tuple [4], [8], [12]:
λ = S, O, A, P, B, π , (1) where S = {s1 , . . . , sN } is a finite set of N hidden states, O = {o1 , . . . , oM } is a finite set of M observation symbols, and π N [j] = πj is the initial probability of state sj . The collection of all transition probabilities forms the transition probability matrix AN ×N , where each entry A[i, j] = aij represents the probability of transitioning from state si to state sj . Similarly, the collection of all duration probabilities forms the duration probability matrix PN ×D , where each entry P[j, d] = pj (d) represents the probability that state sj persists for exactly d consecutive time-steps, with d ∈ {1, . . . , D}. Finally, the emission probabilities define the emission probability matrix BN ×M , where B[j, o] = bj (o) is the probability that state sj emits observation o. These matrix definitions of the model parameters will be central to the tensor formulation presented in Sec. III. Table I summarizes all components of the model. TABLE I: Components of an HSMM and Viterbi Algorithm Parameter

Matrix

Description

S = {s1 , . . . , sN } O = {o1 , . . . , oM } πj aij

SN OM πN AN ×N

pj (d)

PN ×D

bj (o) δt (j)

BN ×M ∆N ×T

ψt (j)

ΨN ×T

Set of N hidden states. Set of M observation symbols. Initial state probability for sj . Transition prob. from si to sj ; self-transitions governed by zero-duration. Prob. that sj persists for exactly d steps, d ∈ [1, D]. Emission prob. of observing o in state sj . Likelihood of the most probable sequence ending in sj at time t. Coordinates (si , d) of the best predecessor for sj at time t.

The Viterbi algorithm is one of the fundamental algorithms for HSMMs. It addresses the following problem: given a model λ and a sequence of observations o0 , . . . , oT −1 , the goal is to find the most likely hidden state sequence s∗ = (q0∗ , . . . , qT∗ −1 ). To achieve this, for every time step t and each state sj ∈ S, the algorithm computes δt (j), which represents the likelihood of the most probable state sequence ending in state sj at time t. Both the sequential and the

tensor formulations (Sec. III) maintain these values in a ∆N ×T matrix: 

∆N ×T

δ1 (0)  δ2 (0)  =   ..  . δN (0)

δ1 (1) δ2 (1) .. . δN (1)



··· ··· .. . ···

δ1 (T −1) δ2 (T −1)   . ..   . δN (T −1)

d−1 h i Y δt (j) = max δt−d (i) · aij · pj (d) · bj (ot−k ) , (3) ∀i,∀d | {z } |{z} | {z } k=0 (a) (b) (c) {z } | (d)

Algorithm 1 Sequential HSMM Viterbi Induction. Following a dynamic programming approach, the algorithm proceeds in three distinct stages: initialization, induction, and backtracking. A. Initialization Phase The initialization phase covers the first D time steps, where D is the maximum admissible state duration. During this interval, we must account for the possibility that the system has occupied state sj since t = 1 with no prior transition. For each state sj and each 1 ≤ t ≤ D, the initialization value combines (a) the initial state probability πj , (b) the probability that state sj persists for exactly t time steps, and (c) the joint emission probability of observations o0 , . . . , ot−1 under state sj . For t > D, no state can have persisted since the beginning, so this contribution is no longer considered.

1: for t = 2 to T do 2: for j = 1 to N do 3: δt (j) ← −∞ 4: for d = 1 to min(t, D) do 5: for i = 1 to N do Q 6: val ← δt−d (i) · aij · pj (d) · d−1 k=0 bj (ot−k ) 7: if val > δt (j) then 8: δt (j) ← val 9: ψt (j) ← (d, i)

Note that induction starts from t = 2; for t ≤ D, the value δt (j) is determined by the maximum of two cases: the system has remained in state j since t = 1 (initialization), or it transitioned from a previous state i at some time τ < t (induction). The algorithm takes the greater of these two probabilities. C. Backtracking Phase

δt (j) = πj · pj (t) · |{z} | {z } (a)

(b)

t−1 Y

bj (ot−τ ) ,

1≤t≤D

(2)

τ =0

|

{z (c)

}

B. Induction phase After initializing the first D time steps, we proceed with the most computationally intensive phase: the induction. For each time step t and each current state sj , the goal is to find the previous state si and the duration d that together maximize the likelihood of reaching sj at time t. The formulation is given by Equation (3), and the corresponding four-nested-loop pseudocode is shown in Algorithm 1. The computation combines four factors. Term (a) is the previously computed value δt−d (i): it encodes the likelihood of the best path ending in state si at time t − d, under the assumption that a transition to sj occurred there. Term (b) is the transition probability aij from state si to state sj . Together, (a) and (b) form the inner maximization: for a fixed duration d, we evaluate all possible source states si and select the one that yields the highest likelihood. The result is then multiplied by term (c) , the probability pj (d) that state sj persists for exactly d consecutive time steps, and by term (d) , the cumulative emission probability of all observations from t − d + 1 to t under state sj . The maximization repeats this for all durations d ∈ {1, . . . , min(t, D)} and all source states si and selects the best combination. The resulting optimal combination (d∗ , i∗ ) for each state sj at time step t is stored in ψt (j), while the corresponding likelihood is stored in δt (j).

Once δt (j) has been computed ∀ t, ∀ 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, until t = 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. III. T ENSOR -BASED V ITERBI A LGORITHM 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, for every time step t we need to evaluate all possible combinations of a previous state si and a duration d for each current state sj . Formally, for every state sj ∈ S we must consider every pair (si , d) with si ∈ S \ {sj } and d ∈ {1, . . . , D}. Once all possible combinations of (si , d) have been computed for each state sj ∈ S, our aim is to compute δt (j) for every time step t ≤ T . To compute these combinations and then δt (j), ∀j sequentially, one must iterate over the current states, then over the possible durations, and then again over the previous states,

yielding three nested loops already contained inside the main time-step loop. The sequential algorithm therefore requires four nested loops. A. The Brick Representation

Algorithm 2 Tensor-based HSMM Viterbi. I NITIALIZATION

(1 ≤ t ≤ D)

1: EN ×D ← Emission Product Computation 2: ∆N ×D ← π N ×↑ ⊙ PN ×D ⊙ EN ×D I NDUCTION (2 ≤ t ≤ T ) ×N ×D 3: BfNirst ← AN ×N ×↑ ⊙ PN ×↑×D

// Brick Construction

4: for t = 2 to T do

Fig. 1: Visualization of the Brick 3D tensor and subdivision into sj-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) can be naturally mapped onto a threedimensional tensor, as shown in Figure 1a. We choose the following layout: • y-axis −→ target states sj , with j ∈ {1, . . . , N }; • x-axis −→ source states si , with i ∈ {1, . . . , N }; • z-axis −→ durations d, with d ∈ {1, . . . , D}. With this convention the set of all combinations (sj , si , d) can be visualized as a 3-dimensional tensor of size N × N × D that we will call Brick (B), shown in Figure 1a. Each slice of the Brick along the y-axis corresponds to the complete set of (si , d) combinations for a single target state sj , which we call sj -slices (Figure 1b). Using this representation, we reformulate the sequential induction of Algorithm 1 into three key stages: (i) the Brick Construction, which occurs once (Sec. III-B); (ii) the Brick Update, where the Brick is modified at each time step t to incorporate temporal factors (as described in Sec. III-C); and (iii) the Maximum Extraction, which identifies the maximum value and the corresponding coordinates (si , d) within each sj -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 that will be described line-by-line below. B. Brick Construction The key observation is the following: the transition probability matrix AN ×N and the duration probability matrix PN ×D do not depend on the time step t. Their contribution to the Brick can therefore be precomputed once, before entering the time-step loop. This corresponds to line 3 of Algorithm 2. As shown in Figure 2a, the transition matrix A is a twodimensional N × N matrix. Within our Brick, it occupies the front face, i.e. it is aligned with the y- and x-axes. The duration probability matrix P is a two-dimensional N × D matrix positioned on the side face, i.e. aligned with the y- and z-axes (Figure 2a).

5:

EN ×D ← Emission Product Computation

6:

×D ∆N ← ∆(t − D : t − 1)N ×D past

7:

×D N ×N ×D B(t) ← Bf irst ⊙ ∆↑×N ⊙ EN ×↑×D // Brick Update past

8: 9:

for j = 1 to N do Ψj (t) ← arg maxd, i B(t) [sj -slice]N ×D

// Past Delta Extraction

// Max. Extraction

BACKTRACKING 10: q∗ = Backtracking(ΨN ×T , ∆N ×T )

Fig. 2: Alignment of A and P within the Brick 3D tensor and the corresponding representation of the broadcasted product.

To combine these two matrices into a single threedimensional tensor, we perform the following product: ×N ×D e N ×N ×↑ ⊙ P e N ×↑×D BfNirst = A | {z } | {z } (a)

(4)

(b)

where the ↑ symbol denotes the axis along which broadcasting occurs, expanding the tensor to match the dimensions of the corresponding operand. Concretely, as shown in Figure 2b-c, the operation can be understood in two steps: (a) Broadcast A: replicate the N ×N matrix D times along e N ×N ×D (Figure 2b-(a)). the z-axis, obtaining a tensor A (b) Broadcast P: replicate the N ×D matrix N times along e N ×N ×D (Figure 2b-(b)). the x-axis, obtaining a tensor P Then, the Brick is given by the element-wise product of the ×N ×D e N ×N ×D ⊙ P e N ×N ×D two broadcast matrices: BfNirst =A (Figure 2c). We refer to this as a broadcasted product: a fundamental operation of our tensor-based Viterbi algorithm.

C. Brick Update After constructing the static Brick Bf irst , we enter the main loop over time-steps. At each step t, two additional matrices must be incorporated: the past delta values ∆past and the emission probability product E, both of which depend strictly on t. This phase corresponds to lines 5–7 of Algorithm 2. 1) Extracting Past Delta Matrix: The past delta values are the simpler of the two time-dependent matrices. We collect from the stored ∆ matrix a window of size N × D obtaining ×D ∆N past .   δ1 (t − D) · · · δ1 (t − 1)  δ (t − D) · · · δ (t − 1)  2   2 N ×D . ∆past =  .. .. ..     . . . δN (t − D) · · · δN (t − 1) As shown in Figure 3, this window starts at t − D (or 1 if t < D) and ends at t − 1, reversed so that the first column corresponds to the most recent past step (duration d = 1). The ×D resulting matrix ∆N past is aligned on the x- and z-axes of the Brick (lower face).

leading us to propose an optimized caching strategy as detailed in Sec. IV.

Fig. 4: Computing Emission Product to obtain the Emission Probability Matrix. Once the past delta matrix ∆past and the emission probability matrix E have been computed for a given time step t, we combine them with the precomputed Brick B through two successive broadcasted products: N ×N ×D ×N ×D e ↑×N ×D ⊙ E e N ×↑×D (5) B(t) = BfNirst ⊙ ∆ past (t) (t) | {z } | {z } (c)

(d)

(c) ∆past lies on the x–z plane, and broadcast the N × D e N ×N ×D . matrix N times along the y-axis, obtaining a tensor ∆ past

Fig. 3: Extracting from ∆ past [t−D, t) values to obtain ∆past . 2) Computing Emission Probability Matrix: The emission probability matrix is the more complex of the two factors. Recall from the term (d) of Equation 3 that, for a given state sj and duration d, we need the product of the emission probabilities over the d most 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 entire N × D matrix of cumulative emission at each time step.  Qd−1 Q1 b1 (ot ) k=0 b1 (ot−k ) . . . k=0 b1 (ot−k ) Qd−1 Q1    b2 (ot ) k=0 b2 (ot−k ) . . . k=0 b2 (ot−k )  N ×D  . E = . .. .. ..  .  ..  . . Q1 Qd−1 bN (ot ) k=0 bN (ot−k ) . . . k=0 bN (ot−k ) 

As shown in Figure 4, we extract the D most recent observation indices, look up the corresponding emission prob for all states inside BN ×M , reverse the order (so that the first column corresponds to duration d = 1). Then, we apply a cumulative product along the duration axis obtaining EN ×D (Figure 4). Our tensor formulation highlights significant redundancies,

(d) E lies on the y–z axis, and broadcast the N ×D matrix e N ×N ×D . N times along the x-axis, obtaining a tensor E N ×N ×D The resulting tensor B is the fully populated Brick for time step t: each entry B[j, i, d] encodes the likelihood of transitioning from state si to state sj with duration d, given the observations up to time t. Having obtained all combinations of (si , d) for each sj , we must now identify the one that yields the maximum likelihood for each sj . D. Maximum Extraction We now need to extract, for each target state sj , the combination (s∗i , d∗ ) that maximizes the tensor entry. Recall that sj is indexed along the y-axis; the corresponding slice is therefore a two-dimensional N ×D matrix. For each sj -slice we seek the maximum value and its associated coordinates (si , d) as shown in Figure 5. This yields the resulting maximum values vector ∆(t)N . This phase corresponds to lines 8–9 of Algorithm 2.

Fig. 5: Argmax computation in each sj -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. E. Initialization Phase As described in Sec. II-A, the initialization phase covers the first D time steps (1 ≤ t ≤ D), accounting for the possibility that state sj has persisted since t = 0 with no prior transition. In the tensor formulation, the entire initialization is expressed as a single broadcasted product of three matrices (lines 1–2 of Algorithm 2): ∆[∀j, 1 : D]N ×D = π N ×↑ ⊙ PN ×D ⊙ EN ×D

(6)

Here, (a) π N is the initial state probability vector, broadcasted along the duration axis; (b) PN ×D is the duration probability matrix; and (d) EN ×D is the cumulative emission matrix, where each entry E[j, d] accumulates the emission probabilities of the first d observations under state sj . Since no entry depends on any other, all N ×D values are computed with no loop-carried dependencies. During these first D steps, both initialization and induction contribute; the final δt (j) is taken as the maximum of the two (lines 8–9 of Algorithm 2). IV. I MPLEMENTATIONS AND O PTIMIZATIONS This section details the optimization strategies for the Tensor-Based Viterbi algorithm from Sec. III. 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: T ENS -P Y, T ENS -1 C, T ENS - MC and T ENS - GPU. T ENS -P Y is a direct transcription of Algorithm 2 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 logspace. Since the quantities involved lie in the interval [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], [13]. The practical consequence is straightforward: products become sums, cumulative products become cumulative sums, and broadcasted products become broadcasted sums. TABLE II: Summary of Implementations Naming

Architecture

Language and Tools

T ENS -P Y

CPU Single-Core

Python, Numpy

T ENS -1 C

CPU Single-Core

C++

T ENS - MC

CPU Multi-Core

C++, OpenMP

T ENS - GPU

GPU

CUDA, HIP

A. NumPy Version T ENS -P Y follows 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 the Brick Computation and Brick Update phases. The Past Delta Extraction phase is performed using NumPy slicing operations, and the argmax is computed using NumPy’s argmax function independently for each destination state sj . 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 warmup phase (t > D), the full emission buffer is obtained by a single element-wise addition and a cache shift, reducing the per-time-step cost from O(DN ) to amortized O(N ) and eliminating the loop-carried dependency on the duration axis. This optimization has been critical in the C++ versions. B. CPU Implementation T ENS -1 C and T ENS - MC are developed in C++ without relying on any tensor library, since neither BLAS [22] nor frameworks such as xTensor [23] provide broadcasted sums natively. Implementing the operations explicitly also enabled us to fuse the Brick Update and Argmax phases, avoiding the need to store the fully populated Brick in 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. 1) Memory Access Patterns: We introduced two complementary flat layouts designed for spatial locality. The ∆past buffer uses a time-major layout (t · N + j) so that all states at a given time-step are contiguous. The ∆ and Ψ arrays use a state-major layout (j · T + t) as required by the Backtracking Phase. The Brick tensor uses a (j, d, i) layout where, for a fixed state j, the entire D × N block 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 states si within a fixed (sj , d) slice, achieving contiguous access. 2) Cached Emissions Computing: We build a 2D emission buffer E indexed by (d · N + j), storing the cumulative sum of emission log-probabilities over each observation window. For t ≤ D, the cumulative sum is computed from scratch. For t > D, the buffer is updated incrementally using an emission cache Ecache : the new observation log-probability at time t is added to a right-shifted copy of Ecache , and the result is saved back to Ecache for the next time-step. This reduces the perstate emission update from O(D) to amortized O(1) after the

warm-up phase, and makes the emission values for all (d, j) pairs available independently before the Argmax loop begins. 3) Fused Brick Update and Maximum Extraction: The two most intensive phases of the algorithm can be fused in this implementation, performing a fused Brick Update + Argmax with 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 ), iterations are fully independent across all three axes, making each axis independently parallelizable or vectorizable. 4) Multi-Core Implementation: The T ENS - MC implementation parallelizes the single-core version using OpenMP. A single #pragma omp parallel region spawns a persistent thread team, using implicit barriers between t steps to avoid repeated fork/join overhead. Pre-computation phases exploit full independence across their iteration spaces. The Brick Construction, for instance, distributes N · D · N independent element-wise sums across all three axes via collapse(3). Similarly, the Cached Emission Computing exposes D · N independent work units via collapse(2). The key parallelization challenge lies in the fused Brick Update and Maximum Extraction, which in the single-core version offers only N independent tasks. We decompose it into two phases: Phase A distributes work over (sj , d) pairs, where each thread sweeps over all source states si to find the local maximum, exposing N · D independent tasks. Phase B then reduces over d per destination state sj . This decomposition enables full thread utilization even when N alone is smaller than the available core count. C. GPU 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. The T ENS - GPU implementation is developed in CUDA and ported to HIP via hipify [24], maintaining a single codebase for both NVIDIA and AMD architectures. The primary computational bottleneck lies in constructing (Bf irst ) and updating (B) the Brick tensor. To maximize hardware utilization, we employ a fine-grained work distribution where each individual thread is responsible for computing a single Brick element. Then, after obtaining the final Brick version within each time-step iteration, the threads operating on elements of the same sj -slice must cooperate to perform the maximum extraction. To implement this mapping, we decompose the N × N × D tensor into N × N vectors of size D aligned along the z-axis, as illustrated in Figure 6. This data layout is mapped onto a 2D grid of N × N thread blocks, where each block manages a single vector. The D elements 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 duration D.

Fig. 6: Mapping the Brick to GPU Thread Blocks. 1) Memory Access Patterns and Coalescing: The emission cache (ordered j · D + d), the brick Bf irst (ordered j · N · D + i · D + d), and ∆ (indexed i · T + (t − 1 − d)) all maintain d-contiguity for consecutive and coalesced thread access. Meanwhile, the emission probability E for the current observation ot and state sj is shared by all threads within a block, allowing the GPU to serve it via a single broadcast read. 2) Cached Emissions Computing: Naively computing E at each time-step requires O(D) work per thread. To eliminate redundant operations, we employ the same cached strategy explored in Sec. IV-B2. This reduces complexity to O(1) via a double-buffering scheme. Thread (j, i, d) reads E[:, d − 1] from Ecache , adds bj (ot ), and writes the result to Ef inal . The value is kept in a register for the current iteration, and the two buffers alternate roles across time-steps. 3) Separated Brick Update and Maximum Extraction: As shown in Sec. III, the tensor formulation decomposes each time-step into two phases: Brick Update and Maximum 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 Bf irst is precomputed once before the induction loop. a) First Kernel: each thread (sj , si , d) computes one element of the final Brick 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) coordinates in O(log D) steps: cross-warp steps use __syncthreads(), while the final log2 (warpsize) steps switch to warp-shuffle instructions (__shfl_down_sync), keeping the running maximum entirely in registers. b) Second Kernel: aggregates the per-block local maximums to determine the global arg max over the full (si , d) plane for each sj . Splitting the computation into two kernels avoids grid-wide synchronization, which would otherwise require cooperativegroups designs with additional occupancy constraints. V. E XPERIMENTAL R ESULTS 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) and validate our implementations (Sec. V-B). Then, we evaluate the performance of both CPU (Sec. V-C) and GPU (Sec. V-D)

1) Baseline: As a baseline for validation and performance comparison, we selected hsmmlearn [16], a C++ library (with a Python API) for HSMMs with explicit duration distributions originating from the R hsmm package [15], 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], 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 as BASE -1 C. We also considered edhsmm [25], which is inspired by hsmmlearn, but its Viterbi algorithm is implemented in Cython, potentially limiting low-level compiler optimization compared to a native C++ implementation. As hsmmlearn is 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 fournested-loops formulation with fully independent iterations. BASE - MC represents the maximum parallelism extractable from the sequential formulation without significant algorithmic restructuring, serving as our primary baseline for both multicore and GPU comparisons. It exhibits near-linear scaling with thread count, provided the number of threads does not exceed N . 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. 2) Problem Size: To evaluate our formulation under realistic conditions, we selected problem sizes guided by computational genomics. We tested configurations with a number of states (N ) ranging from 10 to 75, spanning prokaryotic gene finders [26] at the lower end, chromatin state annotation [27] in the mid-range (15–25), and eukaryotic gene finders such as AUGUSTUS [28] at the upper end. For the sequence length (T ), we tested from 103 to 107 time steps, covering typical gene-finding invocations [5], [29] up to T = 106 [11]. For the maximum duration (D), we adopted values from 100 to 10,000, ranging from the explicit intron duration cutoff of SNAP [5] to stress-test scenarios capturing the longest gene structure features in the human genome [30]. 3) Architectures: We evaluated our implementations on a representative set of high-performance CPU and GPU architectures, summarized in Table III.

CPU

A. Evaluation Environment

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

AMD NVIDIA GPU GPU

implementations, including an investigation into the impact of varying hardware architectures (Sec. V-E). Last, we analyze the energy consumption (Sec. V-F) and conduct a stress test using extreme-scale inputs (Sec. V-G).

Processor

Units

Memory

Compiler

AMD EPYC 7A53 Grace Xeon 8480+

64 Cores 72 Cores 2×56 Cores

512 GiB DDR4 480 GB LPDDR5x 512 GiB DDR5

Cray clang v19.0.0 GCC 14.2 ICX v2024.1.0

A100 SXM H100 SXM H200 SXM

108 SMs 132 SMs 132 SMs

80 GB HBM2e 80 GB HBM3 141 GB HBM3e

CUDA v12.2 CUDA v11.8 CUDA v12.4

MI250X (1 GCD) MI300X

110 CUs 304 CUs

64 GB HBM2e 192 GB HBM3

ROCm v6.3 ROCm v5.7

B. Validation Results All implementations produce identical output to hsmmlearn across every tested configuration, achieving 100% decoding accuracy. To ensure exact equivalence, we include the tail adjustment phase without further optimization, matching hsmmlearn’s boundary handling. Notably, T ENS P Y, our direct NumPy transcription of the tensor formulation (Algorithm 2), already achieves a 4.5× speedup over BASE -1 C. 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. C. CPU Speedup Analysis We begin by comparing T ENS -1 C against BASE -1 C on Intel Xeon 8480+. Figure 7a reports the speedup for T = 105 across N ∈ {10, 15, 25, 50, 75} and D ∈ {100, 250, 500, 1000}. T ENS -1 C is consistently faster, with speedups ranging from 8.7× (N = 75, D = 1000) to 14.1× (N = 10, D = 100). The speedup decreases as either N or D grows, reflecting the point at which D × N exceeds L2 cache capacity.

(a) T ENS -1 C speedup over BASE -1 C and T ENS -1 C runtime (in parentheses; s: seconds; m: minutes). Xeon 8480+, T = 105 .

(b) Baselines and Tensors profiling metrics on Xeon 8480+ CPU with ICX compiler.

Fig. 7: Speedup and runtime of T ENS -1 C and profiling report. We profiled the execution using LIKWID [31], which reveals that the tensor reformulation reduces retired instructions by 36× (from 610B to 16.9B). Figure 7b highlights the resulting hardware efficiency: the vectorization ratio increases from < 0.001% in BASE -1 C to 30.9% in T ENS -1 C, while branch misprediction overhead drops from 6.9% to 0.46%. The

figure also shows a reduced unique DRAM footprint, falling from 2.63 GB to 1.42 GB. 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 larger N and D. Nevertheless, even at the largest configuration, the speedup remains significant: T ENS -1 C completes in 12.2 minutes, whereas BASE -1 C requires over 1.76 hours. We now compare T ENS - MC against BASE - MC on Xeon 8480+. Figure 8 reports the speedup for T = 104 and T = 105 across the same value of D (duration) and N (number of states) considered before. The trend is clear: T ENS - MC becomes increasingly effective as N and D grow. At T = 105 the speedup over BASE - MC reaches 11.4× (N = 50, D = 1000), while at T = 104 a peak of 11.5× is observed at the same configuration. This behavior is the direct consequence of the parallelization strategy. BASE - MC can only distribute work over the N destination states sj , since the sequential formulation carries loop dependencies across durations and source states.

we will evaluate both T ENS - MC and T ENS -1 C across additional CPU architectures, demonstrating that they consistently outperform their respective baselines. Moreover, it is worth remarking that when compiled with GCC, the T ENS - MC runtime for D = 100, N = 10, and T = 104 drops to 0.34 seconds. 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.

Fig. 9: Speedup of T ENS - GPU on H100 over BASE - MC on Xeon 8480+. Cells show T ENS - GPU speedup and runtime. D. GPU Speedup Analysis

Fig. 8: Speedup of T ENS - MC over BASE - MC on Xeon 8480+. The cell shows the speedup and the runtime of T ENS - MC. In contrast, T ENS - MC reformulates the computation as broadcasted sums with no loop-carried dependencies, exposing N · D independent tasks. This exposes sufficient parallelism to fully saturate the 112 available threads, even when N alone would leave most cores idle. LIKWID profiling confirms this efficiency: T ENS - MC retires only 9.9B instructions compared to 21.9B for BASE - MC, while driving L3 bandwidth to nearpeak utilization of the memory hierarchy. Furthermore, Figure 7b shows that T ENS - MC achieves an 80.4% vectorization ratio (versus 0.003% for BASE - MC), reduces bad speculation from 6.25% to 1.16%, and lowers DRAM volume from 4.59 GB to 1.43 GB. At small problem sizes, the trend reverses: for T = 104 with N = 10 and D = 100, the T ENS - MC speedup falls to 0.1×. In this regime, pre-computation and barrier overhead dominate because the per-step workload is insufficient to amortize them. This deficit shrinks to 0.5× at T = 105 and vanishes, becoming a 14× improvement, at T = 106 (as we will show in Fig. 10a), confirming that the overhead is successfully amortized over the longer sequences typical of genomic applications. In Sec. V-E,

Since no GPU implementation of the HSMM Viterbi exists in the literature we use BASE - MC, the multi-core variant we developed by integrating OpenMP parallelization into the sequential C++ backend of hsmmlearn, as our reference. We rely on BASE - MC as a baseline for two reasons: in addition to the total lack of reference GPU implementations, the loopcarried dependencies in the standard sequential formulation across durations and states (si ) 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 reports the resulting speedup of T ENS - GPU (H100) over BASE - MC (Xeon 8480+) for T ∈ {104 , 105 }. T ENS - GPU achieves speedups across the entire parameter space, peaking at 36.8× (N = 25, D = 1000, T = 105 ). Performance scales with both N and D: larger N increases the N × N thread-block grid, better saturating the H100’s 132 SMs, while larger D provides longer reduction vectors per block, improving warp-shuffle efficiency. Even at the smallest configuration (N = 10, D = 100), the speedup remains 3.6– 4.0×, despite only 100 thread blocks being insufficient to fully occupy all SMs. For large N (N = 75), the speedup plateaus (e.g., 17.7× at D = 1000, T = 105 ): the per-block shared memory footprint grows with D, and the second reduction kernel over si states becomes a bottleneck as N increases. At T = 106 the advantage grows further, reaching 54.3× at N = 15, D = 500: BASE - MC requires over 5.7 minutes on 112 CPU cores, while

(a) Speedup of T ENS -1 C over BASE -1 C.

(b) Speedup of T ENS - MC and T ENS - GPU over BASE - MC for N=50.

Fig. 10: Speedup of T ENS -1 C, T ENS - MC, and T ENS - GPU over different CPUs and GPUs, for T=1,000,000.

T ENS - GPU completes the same decoding in 6.4 seconds. Notably, the absolute GPU runtimes remain sub-second for most configurations and never exceed 10 seconds even at the largest tested (N = 75, D = 1000, T = 105 ), 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.

on a H200 completes the same decoding in 24 seconds, a reduction of 570×. 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.

E. Architecture Comparison

F. Energy Consumption

Figure 10a compares T ENS -1 C against BASE -1 C at T = 106 with N ∈ {10, 15, 25} on the three CPUs described in Table III. T ENS -1 C delivers consistent speedups of 10–11× on Grace, 11–14× on Xeon, and 11–12× on EPYC across all analyzed configurations, confirming that the gains are portable across architectures. In absolute terms, for N = 25 and D = 1000, BASE 1 C requires 2.3 hours on Xeon, which T ENS -1 C reduces to 12.4 minutes. For N = 50 and N = 75 (not shown), speedups remain consistent. At N = 75 and D = 500, T ENS -1 C reduces the runtime from 9.2 hours to 56 minutes on Xeon; for D = 1000, BASE -1 C timed out (exceeding 15 hours), whereas T ENS -1 C completed in 2 hours. Grace exhibits the lowest runtime, due to its (almost 2×) higher memory bandwidth. Figure 10b extends the comparison to T ENS - MC and T ENS GPU at T = 106 with N = 50 across all the CPUs and GPUs introduced in Table III. For T ENS - MC, speedups are reported over BASE - MC on the same CPU. For T ENS - GPU, speedups are reported over the CPU where BASE - MC is fastest (Grace). On the CPU side, at D = 1000 T ENS - MC shows a speedup over BASE - MC of 2× on EPYC and 8× on Xeon. On the GPU side, the speedup over BASE - MC grows steadily with D, as increasing the duration expands the per-block workload and improves SM occupancy. At D = 1000, the speedups over BASE - MC range from 5× on MI250X to 12× on 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 at D = 500 remain below 30 seconds for a million-step sequence. For comparison, the original unmodified hsmmlearn baseline BASE -1 C requires 3.8 hours on Grace for this configuration; T ENS - GPU

Figure 11a reports the energy consumption of all implementations at N = 50, T = 104 , normalized to BASE -1 C. For the sake of space, and because we observed a similar trend for the other values of D, we only report the data for D ∈ {100, 1000}. For the CPU versions, we measure the energy on the EPYC 7A53, whereas for T ENS - GPU we measure the energy on the MI250X. In both cases, energy is monitored through the Cray Power Management (PM) counters [32]. 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 over (b) T ENS - GPU performance on a BASE -1 C for N = 50, T = 104 stress test case with T = 107 , on EPYC 7A53 and MI250X. N = 100, D = 104 .

Fig. 11: Energy consumption and stress test. T ENS -1 C reduces energy consumption by ∼10× over BASE -1 C across all tested D values, consistent with its singlecore speedup. In the multi-core regime, T ENS - MC consumes 525.5 J at D = 1000, a 3× reduction compared to BASE - MC

(1.5 kJ), demonstrating that the tensor reformulation translates its runtime advantage into proportional energy savings. T ENS GPU achieves the lowest energy footprint, requiring only 327.7 J at D = 1000 and 85.2 J at D = 100 (just 2% of BASE 1 C). Although the GPU exhibits higher instantaneous power draw, its shorter execution time more than compensates, rendering it the most energy-efficient platform for D > 100. These results confirm that the tensor formulation not only accelerates HSMM Viterbi decoding but also enables a significantly more energy-efficient profile. G. Stress Test: Beyond Current Workloads To demonstrate the practical impact of our formulation, we evaluate T ENS - GPU on an extreme-scale configuration: N = 100 states, D = 10,000 maximum duration, and T = 107 time steps. This scale is entirely inaccessible to the baseline; BASE -1 C triggers a memory allocation failure (std::bad_alloc) for configurations exceeding T = 106 and D = 1,000, precluding direct measurement. By extrapolating from runtimes measured at N = 75, D = 100, T = 106 (where BASE -1 C requires 2 h and BASE - MC 6.7 min) using the O(T · N 2 · D) theoretical complexity, we estimate this extreme configuration would require approximately 148 days for BASE -1 C and 2.1 days for BASE - MC on 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 11b reports the per-iteration runtime of T ENS - GPU on five GPUs. The H200 leads at 53.7 minutes, followed by the MI300X at 57.3 minutes, the H100 at 1.7 hours, the A100 at 2.0 hours, and the MI250X at 3.4 hours. 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, the N × N = 10,000 thread-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 its 110 CUs per GCD, is bottlenecked by its HBM2e bandwidth, the lowest among the five. These results demonstrate that T ENS - GPU reduces 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. VI. R ELATED W ORK The Viterbi algorithm has been fundamental in highperformance 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 par-

allelization difficulty helps explain why performant, hardwareaware HSMM decoders remain absent from the literature. A. Hidden Markov Models (HMMs) HSMMs generalize standard HMMs by introducing explicit state-duration distributions, raising computational complexity from O(T N 2 ) to O(T N 2 D). The Viterbi algorithm has been extensively accelerated for the simpler HMM formulation, including SIMD-vectorized CPU frameworks [18], [20], CUDAbased GPU implementations [19], [33]–[35], hardware– software co-design and domain-specific approaches [17], [36], and distributed computing [37]. However, the additional duration dimension cannot simply be wrapped around existing HMM accelerators: it introduces a cumulative emission product over the d most 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. B. Hidden Semi-Markov Models (HSMMs) Several statistical frameworks implement major HSMM algorithms in R or Python, with performance-critical routines in C/C++ [13]. Domain-specific solutions also exist, such as biomvRhsmm [14] for genomic segmentation, and Pertsinidou and Limnios [38] that propose Viterbi algorithms based on the backward recurrence Markov chain formulation. All of these implementations, however, are sequential and singlethreaded, and none explicitly targets modern high-performance CPUs or GPUs. The most closely related work is Lu et al. [39], 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. C. Summary 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.

VII. D ISCUSSION We now discuss the main design choices behind our formulation, the trade-offs they involve, and the technical directions they leave open. A. Alternative 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 from O(N 2 D) to O(N D + 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) layout that keeps the working set resident across the duration sweep. On GPU, it splits a single joint maximization over the (si , 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. B. Mapping 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. C. Sequence-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.

D. Model Assumptions and Algorithmic Scope Our formulation uses one global maximum duration D for all states, which keeps the tensor dense and the work per thread uniform. Giving each state its own bound Dj 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 ForwardBackward and Baum-Welch procedures iterate over the same (sj , si , 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. VIII. C ONCLUSIONS 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 library tensor-hsmm.1 Across three CPU and five GPU architectures, our implementations achieve speedups of up to 14× on a single core, over 200× with multi-core, and over 570× on GPU with respect to the sequential BASE -1 C baseline, while producing output identical to hsmmlearn in every tested configuration. The gains are structural rather than platform-specific: a direct NumPy transcription of the formulation already outperforms the compiled sequential baseline by 4.5×, and profiling attributes the compiled speedups to a 36× reduction in retired instructions together with vectorization ratios rising from below 0.001% to 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 as 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. 1 https://github.com/HLC-Lab/tensor-hsmm/

R EFERENCES [1] L. Gabriel, T. Brůna, K. J. Hoff, M. Ebel, A. Lomsadze, M. Borodovsky, and M. Stanke, “Braker3: Fully automated genome annotation using rna-seq and protein evidence with genemark-etp, augustus, and tsebra,” Genome research, vol. 34, no. 5, pp. 769–777, 2024. [Online]. Available: https://doi.org/10.1101/gr.278090.123 [2] S. Qin, Z. Tan, and Y. Wu, “On robust estimation of hidden semiMarkov regime-switching models,” Annals of Operations Research, 2024. [Online]. Available: https://doi.org/10.1007/s10479-024-05989-4 [3] H. Zen, K. Tokuda, T. Masuuko, T. Kobayasih, and T. Kitamura, “A hidden semi-markov model-based speech synthesis system,” IEICE TRANSACTIONS on Information, vol. E90-D, no. 5, pp. 825–834, May 2007. [Online]. Available: https://doi.org/10.1093/ietisy/e90-d.5.825 [4] S.-Z. Yu, “Hidden semi-markov models,” Artificial Intelligence, vol. 174, no. 2, pp. 215–243, 2010, special Review Issue. [Online]. Available: https://doi.org/10.1016/j.artint.2009.11.011 [5] I. Korf, “Gene finding in novel genomes,” BMC bioinformatics, vol. 5, no. 1, p. 59, 2004. [Online]. Available: https://doi.org/10.1186/ 1471-2105-5-59 [6] J. Ernst and M. Kellis, “ChromHMM: Automating chromatin-state discovery and characterization,” Nature Methods, vol. 9, no. 3, pp. 215–216, 2012. [Online]. Available: https://doi.org/10.1038/nmeth.1906 [7] R. Durbin, S. R. Eddy, A. Krogh, and G. Mitchison, Biological Sequence Analysis: Probabilistic Models of Proteins and Nucleic Acids. Cambridge University Press, 1998. [Online]. Available: https://doi.org/10.1017/CBO9780511790492 [8] L. Rabiner, “A tutorial on hidden markov models and selected applications in speech recognition,” Proceedings of the IEEE, vol. 77, no. 2, pp. 257–286, 1989. [Online]. Available: https: //doi.org/10.1109/5.18626 [9] L. E. Baum, T. Petrie, G. Soules, and N. Weiss, “A maximization technique occurring in the statistical analysis of probabilistic functions of markov chains,” The Annals of Mathematical Statistics, vol. 41, no. 1, pp. 164–171, 1970. [Online]. Available: http: //www.jstor.org/stable/2239727 [10] A. P. Dempster, N. M. Laird, and D. B. Rubin, “Maximum likelihood from incomplete data via the em algorithm,” Journal of the Royal Statistical Society: Series B (Methodological), vol. 39, no. 1, pp. 1–22, 1977. [Online]. Available: https://doi.org/10.1111/j.2517-6161. 1977.tb01600.x [11] A. Lomsadze, V. Ter-Hovhannisyan, Y. O. Chernoff, and M. Borodovsky, “Gene identification in novel eukaryotic genomes by self-training algorithm,” Nucleic acids research, vol. 33, no. 20, pp. 6494–6506, 2005. [Online]. Available: https://doi.org/10.1093/nar/gki937 [12] Y. Guédon, “Estimating hidden semi-markov chains from discrete sequences,” Journal of Computational and Graphical Statistics, vol. 12, no. 3, pp. 604–639, 2003. [Online]. Available: https: //doi.org/10.1198/1061860032030 [13] C. Bérard, M.-J. CROS, J.-B. DURAND, C. Lothodé, S. Plancade, R. Trepos, and N. Vergne, “Review of hsmm r and python softwares,” A Comprehensive Guide to HSMM: Theory, Software, and Advanced Extensions, pp. 47–77, 2025. [Online]. Available: https://doi.org/10.1002/9781394427581.ch2 [14] Y. Du, E. Murani, S. Ponsuksili, and K. Wimmers, “biomvrhsmm: Genomic segmentation with hidden semi-markov model,” BioMed Research International, vol. 2014, no. 1, p. 910390, 2014. [Online]. Available: https://doi.org/10.1155/2014/910390 [15] J. Bulla and I. Bulla, hsmm: Hidden Semi Markov Models, 2013, r package, archived May 2022. [Online]. Available: https: //cran.r-project.org/package=hsmm [16] J. Vankerschaver, “hsmmlearn: A library for hidden semi-Markov models with explicit durations,” 2021, c++/Cython with Python interface. Wraps the C++ code from the R hsmm package. Archived January 2023. [Online]. Available: https://github.com/jvkersch/ hsmmlearn [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., “Aphmm: Accelerating profile hidden markov models for fast and energy-efficient genome analysis,” ACM Transactions on Architecture and Code Optimization, vol. 21, no. 1, pp. 1–29, 2024. [Online]. Available: https://doi.org/10.1145/3632950

[18] S. R. Eddy, “Accelerated profile hmm searches,” PLoS computational biology, vol. 7, no. 10, p. e1002195, 2011. [Online]. Available: https://doi.org/10.1371/journal.pcbi.1002195 [19] L. Yu, Y. Ukidave, and D. Kaeli, “Gpu-accelerated hmm for speech recognition,” in 2014 43rd International Conference on Parallel Processing Workshops. IEEE, 2014, pp. 395–402. [Online]. Available: https://doi.org/10.1109/ICPPW.2014.59 [20] H. Jiang, N. Ganesan, and Y.-D. Yao, “Cudampf++: A proactive resource exhaustion scheme for accelerating homologous sequence search on cuda-enabled gpu,” IEEE Transactions on Parallel and Distributed Systems, vol. 29, no. 10, p. 2206–2222, Oct. 2018. [Online]. Available: http://dx.doi.org/10.1109/TPDS.2018.2830393 [21] S. Hassan, S. Sarkka, and A. Garcia-Fernandez, “Temporal parallelization of inference in hidden markov models,” IEEE Transactions on Signal Processing, vol. 69, p. 4875–4887, 2021. [Online]. Available: http://dx.doi.org/10.1109/TSP.2021.3103338 [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., “An updated set of basic linear algebra subprograms (BLAS),” ACM Transactions on Mathematical Software, vol. 28, no. 2, pp. 135–151, 2002. [Online]. Available: https://doi.org/10.1145/567806.567807 [23] xtensor-stack, “xtensor: Multi-dimensional arrays with broadcasting and lazy computing,” https://github.com/xtensor-stack/xtensor, 2025, c++14 header-only library. [24] AMD, “HIPIFY: Convert CUDA to portable C++ code,” 2024, accessed: 2026-04-06. [Online]. Available: https://github.com/ROCm/HIPIFY [25] poypoyan, “edhsmm: An(other) implementation of explicit duration hidden semi-Markov models in Python 3,” 2024, python/Cython. Archived May 2024. [Online]. Available: https://github.com/poypoyan/ edhsmm [26] A. V. Lukashin and M. Borodovsky, “Genemark. hmm: new solutions for gene finding,” Nucleic Acids Research, vol. 26, no. 4, pp. 1107–1115, 1998. [Online]. Available: https://doi.org/10.1093/nar/26.4.1107 [27] A. Kundaje, W. Meuleman, J. Ernst, M. Bilenky, A. Yen, P. Kheradpour, Z. Zhang, A. Heravi-Moussavi, Y. Liu, V. Amin et al., “Integrative analysis of 111 reference human epigenomes,” Nature, vol. 518, no. 7539, p. 317, 2015. [Online]. Available: https://doi.org/10.1038/ nature14248 [28] M. Stanke, “Gene prediction with a hidden markov model and a new intron submodel,” Bioinformatics, 2003. [Online]. Available: https://doi.org/10.1093/bioinformatics/btg1080 [29] C. Burge and S. Karlin, “Prediction of complete gene structures in human genomic dna,” Journal of molecular biology, vol. 268, no. 1, pp. 78–94, 1997. [Online]. Available: https://doi.org/10.1006/jmbi.1997. 0951 [30] M. K. Sakharkar, V. T. Chow, and P. Kangueane, “Distributions of exons and introns in the human genome,” In silico biology, vol. 4, no. 4, pp. 387–393, 2004. [Online]. Available: https://doi.org/10.3233/ISB-00142 [31] J. Treibig, G. Hager, and G. Wellein, “LIKWID: A lightweight performance-oriented tool suite for x86 multicore environments,” in Proceedings of PSTI2010, the First International Workshop on Parallel Software Tools and Tool Infrastructures, 2010, pp. 207–216. [Online]. Available: https://doi.org/10.1109/ICPPW.2010.38 [32] HPE Cray, “Cray Performance and Analysis Tools cray pm,” 2026, accessed: 2026-04-06. [Online]. Available: https://cpe.ext.hpe.com/docs/ 24.03/performance-tools/index.html [33] M. HoseinyFarahabady and A. Y. Zomaya, “Gpu-accelerated out-ofcore hmm inference with concurrent cuda streams,” in International Conference on Computational Science. Springer, 2025, pp. 369–376. [Online]. Available: https://doi.org/10.1007/978-3-031-97635-3 44 [34] A. Mohammadidoost and M. Hashemi, “High-throughput and memoryefficient parallel viterbi decoder for convolutional codes on gpu,” arXiv preprint arXiv:2011.09337, 2020. [Online]. Available: https: //doi.org/10.48550/arXiv.2011.09337 [35] V. Roubtsova, “Parallel algorithm for a hidden markov model with an indefinite number of states and heterogeneous observation data.” in IWOCL, 2023, pp. 31–1. [Online]. Available: https: //doi.org/10.1145/3585341.3587954 [36] L. Hummelgren, V. Palmkvist, L. Stjerna, X. Xu, J. Jaldén, and D. Broman, “Trellis: A domain-specific language for hidden markov models with sparse transitions,” in Proceedings of the 17th ACM SIGPLAN International Conference on Software Language Engineering, 2024, pp. 196–209. [Online]. Available: https://doi.org/10. 1145/3687997.3695641

[37] I. Sassi, S. Anter, and A. Bekkhoucha, “Paradist-hmm: A parallel distributed implementation of hidden markov model for big data analytics using spark,” International Journal of Advanced Computer Science and Applications, vol. 12, no. 4, 2021. [Online]. Available: http://dx.doi.org/10.14569/IJACSA.2021.0120438 [38] C.-E. Pertsinidou and N. Limnios, “Viterbi algorithms for hidden semimarkov models with application to dna analysis,” RAIRO-Operations Research, vol. 49, no. 3, pp. 511–526, 2015. [Online]. Available: https://doi.org/10.1051/ro/2014053 [39] Z. Lu, L. T. Yang, A. Azman, F. Zhou, S. Zhang, and X. Fu, “Tensorbased hidden semi-markov model for cpss user activity analysis and services,” IEEE Transactions on Services Computing, 2025. [Online]. Available: https://doi.org/10.1109/TSC.2025.3618011