PaperlayerView exact revision on arXiv

Source-linked AI summary

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

Lorenzo Piarulli, Elia Belli, Daniele De Sensi

arXiv:2609.16500v1cs.LGcs.DCcs.DS

TL;DR

HSMM Viterbi implementations have remained sequential, despite the added duration dimension that makes decoding costly. This paper reformulates the algorithm as tensor operations and implements it on CPUs and GPUs, reporting speedups of up to 570× over a sequential baseline while matching tested outputs.

  • Problem

    HSMM Viterbi implementations remain sequential, and the added duration loop raises computational complexity to O(TN^2D), making large-scale decoding costly.

    Source: §I · p. 2

  • Method

    The paper restructures the Viterbi algorithm's inner loops into dense tensor operations and provides optimized single- and multi-core CPU and GPU implementations.

    Source: §I · p. 2, §VIII · p. 12

  • Results

    Up to 14× on a single core, over 200× with multi-core, and over 570× on GPU versus the sequential BASE-1C baseline across three CPU and five GPU architectures.

    Source: §VIII · p. 12

  • Takeaways & Limitations

    The implementations produced output identical to hsmmlearn in every tested configuration, with 100% decoding accuracy.

    Source: B. Validation Results · p. 8, §VIII · p. 12

    Maximum-based reductions do not directly map to the multiply-accumulate pipelines of specialized accelerators and require a dedicated mapping strategy.

    Source: B. Mapping onto Specialized Accelerators · p. 12

Abstract

from arXiv · show

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 14x on a single core, over 200x with multi-core, and over 570x on GPU over the state-of-the-art sequential baseline, establishing a new performance baseline for large-scale HSMM decoding.

I. INTRODUCTION

HSMMs explicitly model variable state durations, but their Viterbi algorithm adds a costly duration dimension and has remained confined to sequential CPU execution. The paper introduces a tensor formulation and optimized implementations to address this gap.

I. INTRODUCTION

  • HSMMs model variable state durations explicitly, unlike standard HMMs' geometric, memoryless state occupancy.The model also includes state transitions and state-specific emissions.Source: §I · p. 1
  • Viterbi decoding finds the most likely hidden state sequence and can be iterated with parameter re-estimation for ab initio gene prediction.The paper identifies Viterbi as the most widely adopted of the three fundamental HSMM algorithms because it serves both decoding and ab initio learning.Source: §I · p. 2
  • O(TN^2D), with D ranging from hundreds to thousands, makes HSMM Viterbi computationally costly, while existing implementations remain scalar, single-threaded CPU code.The duration loop distinguishes HSMM complexity from the O(TN^2) complexity of standard HMM Viterbi.Source: §I · p. 2
  • The tensor formulation restructures the algorithm's inner loops as dense operations suited to SIMD and SPMD, with optimized CPU and first-time GPU implementations.The authors report speedups of up to 570× over state-of-the-art HSMM frameworks.Source: §I · p. 2

II. SEQUENTIAL HSMM VITERBI ALGORITHM

  • The sequential HSMM Viterbi algorithm computes the likelihood of the best path ending in each state and proceeds through initialization, induction, and backtracking.Both sequential and tensor formulations maintain these likelihoods in a state-by-time matrix.Source: §II · pp. 2–3

A. Initialization Phase

Initialization scores paths that may have occupied a state since the sequence began; induction then evaluates transitions from earlier states across admissible durations.

A. Initialization Phase

  • For the first D time steps, initialization scores paths that began in state j at time 1 using its initial probability, duration probability, and observation emissions.For t > D, the path that has persisted since the beginning is no longer considered.Source: A. Initialization Phase · p. 3

B. Induction phase

  • Induction tests prior states and durations, combining the prior path likelihood, transition probability, duration probability, and emissions, then stores the best likelihood and (d*, i*) for each target state.For t ≤ D, the algorithm takes the greater of the initialization case and the induction case.Source: B. Induction phase · p. 3

C. Backtracking Phase

The tensor formulation restructures induction with tensor operations, while backtracking remains sequential because it is not a computational bottleneck. A 3D Brick organizes target states, source states, and durations for the induction computation.

C. Backtracking Phase

  • Backtracking starts from the final time step's highest-delta state and proceeds backward to recover the most likely state chain.The recovered path typically has states persisting across multiple consecutive time steps, reflecting explicit duration modeling.Source: C. Backtracking Phase · p. 3

III. TENSOR-BASED VITERBI ALGORITHM

  • The tensor formulation introduces new formulations for initialization and induction, but leaves backtracking sequential because it is not a computational bottleneck.Optimization focuses on induction, where tensor reformulation yields the greatest benefit.Source: §III · p. 3

From Loops to Tensors

  • Each time step considers every previous-state and duration pair for each current state, requiring four nested loops in the sequential algorithm.The three inner loops iterate over current states, durations, and previous states.Source: From Loops to Tensors · pp. 3–4

A. The Brick Representation

  • The Brick maps target states, source states, and durations onto a 3D tensor, with each target-state slice containing all source-state and duration combinations.Brick construction is precomputed once; update and maximum extraction repeat at each time step.Source: A. The Brick Representation · p. 4

B. Brick Construction

Brick construction combines the time-independent transition and duration probabilities in a precomputed tensor. Broadcasting aligns both matrices with the Brick dimensions before their element-wise product.

B. Brick Construction

  • The transition matrix A and duration matrix P are time-independent, so their contribution to the Brick is precomputed before the time-step loop.A aligns with the y- and x-axes; P aligns with the y- and z-axes.Source: B. Brick Construction · p. 4
  • The broadcasted product expands matrices along absent axes so their dimensions match for combination into one three-dimensional tensor.The operation is presented as two broadcasts followed by an element-wise product.Source: B. Brick Construction · p. 4
  • Broadcast A across the duration axis and P across the source-state axis, then multiply the expanded tensors element-wise to construct the Brick.Each expanded matrix has dimensions N×N×D.Source: B. Brick Construction · p. 4

C. Brick Update

At each time step, the Brick is updated with past delta values and cumulative emission probabilities. These matrices are aligned and broadcast into the Brick before extracting the likelihoods for each target state.

C. Brick Update

  • The past-delta matrix is a reversed window ending at t−1, with its first column representing the most recent past step, corresponding to duration 1.The window starts at t−D, or at 1 when t < D, and has size N×D.Source: 1) Extracting Past Delta Matrix · p. 5
  • Because tensor formulation eliminates the duration loop, it computes an N×D matrix of cumulative emission probabilities at every time step.It looks up emissions for the D most recent observations across states, reverses their order, and takes a cumulative product along the duration axis.Source: 2) Computing Emission Probability Matrix · p. 5
  • The emission computation exposes redundancies, motivating an optimized caching strategy described in Sec. IV.Source: 2) Computing Emission Probability Matrix · p. 5
  • The time-dependent past-delta and emission matrices are broadcast into the Brick, whose entries encode transition likelihoods for each source state, target state, and duration.The past-delta values align on the Brick's x–z plane and the emission values on its y–z axis.Source: 2) Computing Emission Probability Matrix · p. 5, C. Brick Update · p. 5

D. Maximum Extraction

Maximum extraction selects the best source-state and duration combination for each target state. Initialization handles the first D steps with a broadcasted product, while initialization and induction both contribute to the final delta values during that period.

D. Maximum Extraction

  • For each target-state slice, maximum extraction returns the largest likelihood and its source-state and duration coordinates.The resulting maximum values form the delta vector for that time step.Source: D. Maximum Extraction · p. 5, Figure 5 · p. 5

E. Initialization Phase

  • Initialization covers the first D time steps, accounting for a state that may have persisted since time zero without a prior transition.Source: E. Initialization Phase · p. 6
  • Initialization combines initial-state probabilities, duration probabilities, and cumulative emissions in a broadcasted product; during the first D steps, delta takes the maximum of initialization and induction.The N×D entries are computed without loop-carried dependencies.Source: E. Initialization Phase · p. 6

IV. IMPLEMENTATIONS AND OPTIMIZATIONS

The implementations exploit the tensor formulation through cache-aware computation and parallelism, from a NumPy reference to optimized C++ CPU versions. The CPU implementations improve memory access, reuse emissions, and parallelize independent work.

IV. IMPLEMENTATIONS AND OPTIMIZATIONS

  • The tensor formulation exposes data reuse and independent computation axes that support cache-friendly access, SIMD, and multithreaded execution.Source: §IV · p. 6

A. NumPy Version

  • TENS-PY provides an unoptimized NumPy reference implementation for expressing and validating the tensor algorithm.It uses broadcasted sums, slicing for past-delta extraction, and NumPy argmax by destination state.Source: §IV · p. 6, A. NumPy Version · p. 6
  • After warm-up, a rolling emission cache reduces per-time-step emission accumulation from O(DN) to amortized O(N).The cache update decouples emission accumulation from the main recurrence and removes the duration-axis loop-carried dependency.Source: A. NumPy Version · p. 6

B. CPU Implementation

  • The C++ CPU versions implement tensor operations directly and fuse Brick Update with Argmax to avoid storing the full Brick.Source: B. CPU Implementation · p. 6
  • The CPU layouts place time-step states contiguously in ∆past, arrange ∆ and Ψ for backtracking, and keep each state’s Brick block contiguous.The Brick layout supports contiguous access across source states within a fixed destination-state and duration slice.Source: B. CPU Implementation · p. 6
  • TENS-MC uses a persistent OpenMP team and splits fused maximum extraction into N·D local tasks before reducing over durations.This exposes enough work to use available cores even when the number of states is below the core count.Source: 4) Multi-Core Implementation · p. 7

C. GPU Implementation

The paper presents its first GPU-accelerated HSMM Viterbi implementation, mapping Brick computation to thread blocks and splitting update and maximum extraction across kernels. The experimental-results section then outlines validation, performance, architecture, energy, and stress-test analyses.

C. GPU Implementation

  • The authors present TENS-GPU as, to their knowledge, the first GPU-accelerated Viterbi implementation for HSMMs, using CUDA and a HIP-portable codebase.Source: C. GPU Implementation · p. 7
  • Each GPU thread computes one Brick element, while threads within a destination-state slice cooperate to extract maxima.Source: C. GPU Implementation · p. 7
  • The N × N × D Brick is arranged as N × N vectors of length D, with each thread block handling one vector and distributing its D elements among threads.The mapping uses an N × N grid of blocks, each handling a temporal slice.Source: C. GPU Implementation · p. 7, Figure 6 · p. 7
  • The GPU reuses cached emissions and separates Brick Update from maximum extraction into two kernels, with a block-level reduction followed by a global argmax.The first kernel finds local maxima and their coordinates; the second aggregates those maxima for each destination state.Source: 1) Memory Access Patterns · p. 7, 3) Separated Brick Update and Maximum… · p. 7

V. EXPERIMENTAL RESULTS

  • The evaluation covers implementation validation, CPU and GPU performance, architecture variation, energy consumption, and extreme-scale inputs.Source: §V · pp. 7–8

1) Baseline:

The evaluation uses hsmmlearn as its compiled single-core baseline and tests state, sequence-length, and duration ranges motivated by computational genomics. Implementations match hsmmlearn’s outputs, while the CPU tensor versions show substantial speedups and profiling differences.

1) Baseline:

  • The single-core baseline is hsmmlearn, a C++ library for explicit-duration HSMMs; BASE-MC adds OpenMP parallelization for multicore and GPU comparisons.BASE-MC’s parallelism is limited to destination states and saturates or degrades slightly when thread count exceeds N.Source: 1) Baseline · p. 8

2) Problem Size:

  • The tested ranges are N = 10–75, T = 10^3–10^7, and D = 100–10,000, selected to reflect computational-genomics applications and stress tests.Source: 2) Problem Size · p. 8

3) Architectures:

B. Validation Results

  • All implementations match hsmmlearn’s output across tested configurations, achieving 100% decoding accuracy.The comparison includes the tail adjustment phase to match hsmmlearn’s boundary handling.Source: B. Validation Results · p. 8

C. CPU Speedup Analysis

  • 14.1× speedup is reported for TENS-1C over BASE-1C at N = 10, D = 100 on Xeon 8480+; its speedups range from 8.7× to 14.1× across the tested configurations.At the largest configuration, TENS-1C takes 12.2 minutes versus over 1.76 hours for BASE-1C.Source: Figure 7 · p. 8, C. CPU Speedup Analysis · pp. 8–9
  • TENS-1C reduces retired instructions by 36× and raises vectorization from < 0.001% to 30.9% versus BASE-1C.The profiling passage also reports branch-misprediction overhead dropping from 6.9% to 0.46%.Source: C. CPU Speedup Analysis · p. 8
  • 11.5× is the reported peak speedup for TENS-MC over BASE-MC at N = 50, D = 1000, T = 10^4 on Xeon 8480+.At T = 10^5, the speedup at that configuration is 11.4×.Source: Figure 8 · p. 9, C. CPU Speedup Analysis · p. 9
  • TENS-MC exposes N·D independent tasks and reaches 80.4% vectorization versus 0.003% for BASE-MC, but small workloads can be slower when setup and barrier overhead dominate.At T = 10^4, N = 10, D = 100, the speedup is 0.1×; at T = 10^6 it becomes a 14× improvement.Source: C. CPU Speedup Analysis · p. 9

D. GPU Speedup Analysis

TENS-GPU outperforms the multicore CPU baseline across the tested parameter space, and the architecture comparison reports large gains across GPUs and CPUs. The authors note that the broad comparison uses shared codebases rather than architecture-specific specialization.

D. GPU Speedup Analysis

  • At T = 10^5, TENS-GPU speedups range from 3.6–4.0× at N = 10, D = 100 to 36.8× at N = 25, D = 1000 over BASE-MC.At N = 75, speedup plateaus as shared-memory footprint and the second reduction kernel become bottlenecks.Source: D. GPU Speedup Analysis · p. 9
  • 54.3× is TENS-GPU’s reported speedup over BASE-MC at N = 15, D = 500, T = 10^6; the GPU takes 6.4 seconds while BASE-MC takes over 5.7 minutes.Source: D. GPU Speedup Analysis · pp. 9–10

E. Architecture Comparison

  • At D = 1000, TENS-MC is 2× faster than BASE-MC on EPYC and 8× on Xeon, while GPU speedups range from 5× on MI250X to 12× on H200.For TENS-1C at T = 10^6, speedups are 10–11× on Grace, 11–14× on Xeon, and 11–12× on EPYC.Source: E. Architecture Comparison · p. 10
  • 570× is the reported reduction over BASE-1C for a million-step decoding on H200: 24 seconds versus 3.8 hours on Grace.Across GPUs, the reported runtime at D = 500 remains below 30 seconds for a million-step sequence.Source: E. Architecture Comparison · p. 10
  • The architecture comparison uses one CUDA/HIP GPU codebase and OpenMP-only multicore parallelization, leaving architecture-specific memory and instruction tuning as possible further specialization.Source: E. Architecture Comparison · p. 10

F. Energy Consumption

The tensor implementation reduces energy consumption, and the stress test shows that a workload estimated to take over a month on one core can finish in under an hour on one GPU. The related-work discussion explains why HSMM decoding has remained difficult to accelerate.

F. Energy Consumption

  • TENS-GPU uses 327.7 J at D=1000 and 85.2 J at D=100, just 2% of BASE-1C energy at D=100.TENS-1C reduces energy by ∼10× across tested D values; TENS-MC uses 525.5 J at D =1000, a 3× reduction versus BASE-MC.Source: Figure 11 · p. 10, F. Energy Consumption · pp. 10–11

G. Stress Test: Beyond Current Workloads

  • TENS-GPU reduces a 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 practical.The extreme-scale configuration uses N = 100 states, D = 10,000 maximum duration, and T = 10^7 time steps; BASE-1C triggers a memory allocation failure for configurations exceeding T = 10^6.Source: G. Stress Test: Beyond Current Workloads · p. 11

VI. RELATED WORK

  • HSMM acceleration is difficult because duration handling adds a computational burden and makes the inner loops difficult to parallelize, helping explain the lack of performant hardware-aware decoders.The related-work discussion notes that HSMM decoders remained sequential and single-threaded, unlike the acceleration efforts for standard HMMs.Source: §VI · p. 11

A. Hidden Markov Models (HMMs)

  • Standard HMM accelerators do not extend directly to HSMMs, whose duration dimension requires cumulative emission products, variable-depth delta access, and joint state-duration maximization.The added duration dimension raises complexity from O(TN 2) to O(TN 2D).Source: A. Hidden Markov Models (HMMs) · p. 11

B. Hidden Semi-Markov Models (HSMMs)

The paper distinguishes its exact, acceleration-focused HSMM Viterbi formulation from prior tensor-based modeling work, and discusses algorithmic and hardware trade-offs and scope boundaries. Across CPU and GPU architectures, it reports large speedups over the sequential baseline while preserving outputs identical to hsmmlearn.

B. Hidden Semi-Markov Models (HSMMs)

  • Earlier HSMM frameworks remained sequential; Lu et al.'s T-HSMM models correlated entities and collapses duration to expected values, unlike this paper's exact decoding.The passage says a direct head-to-head performance comparison is not meaningful because the methods solve different problems.Source: B. Hidden Semi-Markov Models (HSMMs) · p. 11

C. Summary

  • The paper identifies a gap in high-performance HSMM Viterbi implementations and proposes a tensor reformulation as a basis for accelerator implementations and application-specific mappings.Source: C. Summary · p. 11

A. Alternative Algorithmic Formulations

  • A factorized formulation could reduce per-step work from O(N 2D) to O(ND + N 2), but the authors reject it because of hardware-efficiency costs.On GPU it requires two sequential reductions and an extra grid-wide synchronization per time step; on CPU it removes the Brick precomputation and contiguous layout.Source: A. Alternative Algorithmic Formulations · p. 12

B. Mapping onto Specialized Accelerators

C. Sequence-Level and Distributed Parallelism

  • The implementations decode one sequence at a time; chunked sequence decoding and multi-node execution are proposed extensions, while batching independent sequences could keep devices busier on small inputs.Source: C. Sequence-Level and Distributed… · p. 12

D. Model Assumptions and Algorithmic Scope

  • The formulation uses one global maximum duration and assumes discrete emissions; state-specific duration bounds and continuous emissions would require different handling.Forward-Backward and Baum-Welch procedures are described as reusing the same tensor operations with a sum replacing the maximum.Source: D. Model Assumptions and Algorithmic… · p. 12

VIII. CONCLUSIONS

  • The formulation replaces three sequential inner loops with dense tensor operations to expose optimization opportunities, and supports optimized single- and multi-core CPUs and the first HSMM GPU implementation.The implementations are released as the open-source library tensor-hsmm.Source: §VIII · p. 12
  • Up to 14× single-core, over 200× multi-core, and over 570× GPU speedups over sequential BASE-1C, with outputs identical to hsmmlearn in every tested configuration.The conclusion reports evaluation across three CPU and five GPU architectures.Source: §VIII · p. 12
Loading 2609.16500v1…