Discovering group dynamics in coordinated time series via hierarchical recurrent switching-state models
Authors: Michael Wojnowicz
Kaitlin Gili
Preetish Rath
Eric Miller
Jeffrey W. Miller
Clifford Lee Hancock
Meghan O'Donovan
Seth Elkin-Frankston
Tad Brunye
Michael C Hughes
Venue: Accepted by TMLR ()

Abstract: We seek a computationally efficient model for a collection of time series arising from multiple interacting entities (a.k.a. "agents"). Recent models of temporal patterns across individuals fail to incorporate explicit system-level collective behavior that can influence the trajectories of individual entities. To address this gap in the literature, we present a new hierarchical switching-state model that can be trained in an unsupervised fashion to simultaneously learn both system-level and individual-level dynamics. We employ a latent system-level discrete state Markov chain that provides top-down influence on latent entity-level chains which in turn govern the emission of each observed time series. Recurrent feedback from the observations to the latent chains at both entity and system levels allows recent situational context to inform how dynamics unfold at all levels in bottom-up fashion. We hypothesize that including both top-down and bottom-up influences on group dynamics will improve interpretability of the learned dynamics and reduce error when forecasting. Our hierarchical switching recurrent dynamical model can be learned via closed-form variational coordinate ascent updates to all latent chains that scale linearly in the number of entities. This is asymptotically no more costly than fitting a separate model for each entity. Analysis of both synthetic data and real basketball team movements suggests our lean parametric model can achieve competitive forecasts compared to larger neural network models that require far more computational resources. Further experiments on soldier data as well as a synthetic task with 64 cooperating entities show how our approach can yield interpretable insights about team dynamics over time.

Published in Transactions on Machine Learning Research (08/2025)
Discovering group dynamics in coordinated time series via hierarchical recurrent switching-state models
Michael T. Wojnowiez Kaitlin Gili', Preetish Rath', Eric Miller!, Jeffrey Miller?, Clifford Hancock!, Meghan O'Donovan', Seth Blkin-Frankston! 1.4, Tad T. Brunye' and Michael C. Hughes
Tufts University, Medford, MA, USA
Haroard Universily, Boston, MA, USA
Montana State University, Bozeman, MT, USA
U.S. Army Combat Capabilities Development Command Soldier Center (CCDC SC), Natick, MA, USA
Reviewed on OpenReview: https: // openreview. net/ permn? id= Lcchhthth
Abstract
We seek computationally efficient model for collection of time series arising from multiple interacting entities (a.k.a. "agents"). Recent models of temporal patterns across individuals fail to incorporate explicit system-level collective behavior that can influence the trajectories of individual entities. To address this gap in the literature, we present a new hierarchical switching-state model that can be trained in an unsupervised fashion to simultancously learn both system-level and individual-level dynamics. We employ latent system-level discrete state Markov chain that provides top-down influence on latent entity-level chains which in turn govern the emission of each observed time series. Recurrent feedback from the observations to the latent chains at both entity and system levels allows recent situational context to inform how dynamics unfold at all levels in bot tom-up fashion. We hypothesize that including both top-down and bottom-up influences on group dynamics will improve interpretability of the learned dynamics and reduce error when forecasting, Our hierarchical switching recurent dynamical model can be learned via closed-form variational coordinate ascent updates all latent chains that scale lincarly the number of entities. This asymptotically no more costly than fitting separate model for each entity. Analysis of both synthetic data and real basketball team movements suggests our lean parametric model can achieve competitive forecasts compared to larger neural network models that require far more computational resources. Further experiments on soldier data as well as synthetic lask with 64 cooperating entities show how our approach can yield interpretable insights about team dynamics over time.
1 Introduction
We consider the problem of jointly modeling collection of multivariate time series arising from individual entities that can influence each other's behavior over time and might share goals. Each series in the collection describes the evolution of one entity, sometimes also called an "agent" (Yuan et al., 2021). All entities are observed over the same time period within a shared environment or system. Our work motivated by the need to capture an essential property of such data in many applications: the temporal behaviors of the individual entities are coordinated in a systematic but fundamentally latent (i.e., unobserved) manner.
As motivating example, consider the dynamics of team sport like basketball (Terner & Franks, 2021). To accomplish team goal, one player might set "screen," physically blocking a defender to allow teammate an open drive to the basket. This screen could be preplanned or arise players act on an in-moment opportunity. As another example, consider marching band players that practice moving logether in a coordinated fashion
Open-source code: https://github.. com/ tff tonfft-mll testttssmmmmaameeeeeeeee www.m//ddddaanics---e------ time- sime-seriis/

---

Published in Transactions on Machine Learning Research (08/2025)
across a field. Most movements are planned out by coach. However, in some cases, situational adjustments are needed: when enough individuals make mistakes in their own trajectories, the coach may decide to rein in all of the players and reset. Our experiments directly cover basketball and marching band case studies later in Sec. 5.2 and Sec. 5.3, respectively. Similar group dynamics arise in many other domains, such as businesses in a common economic system (van Dijk et al., 2002), animals in a shared habitat. (Sun el al., 2021), biological cells sharing copy-number mutations (Babadi et al., 2023), or students a collaborative problem-solving session (Odden & Russ, 2019; Earle-Randell al., 2023).
In all these examples, the dynamics of the individuals are far from independent. Instead, the observed trajectories exhibit "top-down" patterns of coordination that are planned in advance or learned from extensive training together. They also exhibit "bottom-up" adaptations of the individuals and the group to evolving situational demands. We seek to build a model that can infer how group dynamics evolve over time given only entity-level sensory measurements, while taking into account top-down and bottom-up influences.
While modeling individual time series has seen many recent advances (Linderman et al., 2017; Gu el al., 2022; Farnoosh et al., 2021), there remains need for improved models for coordinated collections of time series. To try to model a collection of time series, a simple approach could repurpose models for individual time scries. As a step beyond this, some efforts pursue personalized models that allow custom parameters that govern each entity's 's dynamics while sharing information ween entities via common priors on these parameters (Severson et al., 2020; Linderman et al., 2019; Alaa & van der Schaar, 2019), often using mixed effects (Altman, 2007; Liu et al., 2011). But personalized models allow each sequence to unfold asynchronously without interaction. In contrast, our goal is to specifically model coordinated behavior within the same time period. Others have pursued this goal with complex neural architectures that can jointly model "multi-agent" " trajectories (Zhan et al., 2019; Alcorn & Nguyen, 2021; Xu et al., 2022). Instead, we focus on parametric methods that are easier to interpret and more likely to provide sample-eficient quality fits in applications with only a few minutes of available data (such as the data described in Sec. 5.4).
One potential barrier to modeling coordination across entities is computational complexity. For instance, a model with discrete hidden states which allows interactions among entities has a factorial structure with inference that scales exponentially in the number of entities (see Sec. 3). In this paper we present a tractable framework for modeling collections of time scries that overcomes this barrier. All estimation can be done with cost. linear in the number of entities, making our model's 's asymptotic runtime complexity no more costly than fitting separate models to each entity.
The first key modeling contribution is an explicit representation of the hierarchical structure of group dynamics, modeled via a latent system-level state that exerts "top-down" " influence on each individual time series. We use well-known switching-state models (Rabiner, 1989) as a building block for both system-level and individual-level dynamics. As shown in Fig. 1, our model posits two levels of latent discrete state chains: system-level chain shared by all entities and an entity-level chain unique to each entity. We assume that the system-level chain is the sole mediator of cross-entity coordination; each entity-level chain is conditionally independent of other entities given the system-level chain. Our model achieves "top-down" coordination via the system-level chain's influence on entily-level state transition dynamics. In turn, an emission model produces each entity's observed time series given the entity-level chain.
The second key modeling contribution to allow "bottom-up" adjustments to recent situational demands at all levels of the hierarchy. In our assumed generative model in Fig. 1, the next system-level state and entity-level state both depend on feedback from per-entity observed data at the previous timestep. In the basketball context, this feedback captures how a basketball player driving to the basket switches to another behavior after reaching their goal. Critically, with our model that situational change in one player can quickly influence the trajectories of other players. We refer to observation-to--aten feedback as recurrent, following Linderman et al. (2017). Note that this recurrent feedback also allows the model to learn state duration distributions that are more expressive than the geometric distributions that would be implied otherwise (A (Ansari et al., 2021). Previously, Linderman et al. (2017) incorporated such recurrent feedback into model with flat single-level of switching states. We show how recurrent feedback can inform a two-level hicrarchy of system-level and entity-level states, so that entity-level observations can drive system-level transitions in a computationally efficient manner.
2

---

Published in Transactions on Machine Learning Research (08/2025)
Discrete states for system
Disurete states for entity
Observations
D
212
Figure 1: Graphical model representation of our hierarchical switching recurvent dynamiiaa model (HSRDM) for system of J interacting entilies. Colored edges highlight key insights behind our flexible transition models of system-level hidden states and entity-level hidden states 2. Transitions to the next system state depend via blue arrows on the current system state and recurrent feedback from observations of all entities (up-diagonal). Transitions to the next entity state depend via red arrows the next system-level state (down arrow), current entity state (horizontal), and recurrent feedback from that entity's observations (up-diagonal).
Our overall contribution is thus a proposed framework hierarchical switcching recurrent dynamical models by which our two key modeling ideas provide a natural and cost-effective solution the problem of unsupervised modeling of coordinated time series. Unlike other models, framework allows each entity's next-step dynamics to be driven by both system-level discrete state ("top-down" influence) and recurrent feedback from previous observations ("bottom-up" influence). Optional exogenous features (e.g. the ball position in basketball) can also be incorporated. We further provide variational inference algorithm for simultancously estimating model parameters and approximate posteriors over system-level and entity-level chains. Each chain's posterior maintains the model's temporal dependency structure while remaining affordable via efficient dynamic programming that incorporates recurrent feedback. We conduct experiments on two synthetic datasets well as two real-world tasks to demonstrate the model's superior capability in discovering hidden dynamics as well as its competitive forecasts obtained at low computational cost.
2 Model Family
Here we present a family of hierarchical switching recurrent dynamical models (HSRDMs) to describe a collection of time series gathered from entitics that interact over a common time period (discretized into timesteps E (0, 1,2,.... T)) and in a common enviroument or system. For each entily, indexed by (1,..., J}, we observe a time series of feature vectors {x( RR, 0,1,
Our HSRDM represents the j-th entity via two random variables: the observed or "emitted" fcatures oor and hidden entity-level discrete state sequence < {2{ {1. K2), J k 0, 0..., T). We further assume system-level latent time scries of discrete states (s E f1...., L), t= , The complete joint density of all random variables across all J entities factorizes according to the graphical model in Figure as
T
(1:J) (1:J) p( "0:T" 8O.T ()
p(so
II --11 (1:3)
II j=l
p 20:7' So:r; 0)
(2.1)
system state init.
[=1 system state transitions
per-entity states and emissions
Here, 0 denotes all model parameters, and superscript. (1:J) denotes the concatenation of variables over all entities. For a specific entity index j, its individual states and emitted data factorize
T P( 2 8t,
T PT m0 22
(1) (0) p( Porr sO.T, 0) = p(=0
80:
p(a =0
1
eutity state init.
entity state transitions
emission init.
emission dynamics
The design principle of HSRDMs is coordinate the switching-stard state dynamics of multiple entities so they receive top-down influence from system-level state well bottom-up influence via recurrent feedback from
3

---

Published in Transactions on Machine Learning Research (08/2025)
entity observations (diagonal up arrows in Fig. 1). Under the generative model, the next entity-level state depends on the interaction of three sources of information: the next state of the system, the current state of the entity, and the curent entity observation. Likewise, the next system state depends on the current system state and observations from all entities.
Transition models. To instantiate our two-level discrete state transition distributions, we use categorical generalized linear models incorporate each source of information via additive utilities
$1-1, (1:J) ~ Cat-GLM&
A gv (2 (1:J) ,
(2.2a)
transition preferences
recurrent ferdback
(
(1) 2-11 8118 Cat-GLMK
(P) (
(2.2b)
recurrent foedback
transition preferences
Across levels, common sources of information drive these utilities. First, the state-to-state transition term selects an appropriate log transition probability vector from matrices II, P via one-hot vector eş indicating the previous stale k. Second, recurrent feedback governs the next term, via leaturization functions for the system gv : DJ + R and for entities fs RDD + RRD with parameters v & (known or learned) and weights A, y.
Some applications may benefit from using available ezogenous covariates to inform recurrent fcodback, such as using the location of the ball itself the remaining time on the game clock better model future basketball player movement. Our general framework denotes such covariates at the system-level # or entity-level as - When available, these covariates can drive transition probabilities; no such variables exist they can be left out. Note that inference (Sec. 3) applies not merely to Equation 2.2, but to arbitrary instantiations.
Emission model. We generate the next observation for entity j via a state-conditioned autoregression:
(1) Hc. where ( (=<(0(21)=)))
(2.3)
Users can select the emission distribution иl to match the domain of observed features a(") our later experiments Gaussians for real-valued vectors and Von-Mises distributions for angles. The parameter ( < 2())) of the chosen H depends on the previous observation a.(2) and current entity-level state ( We focus on lag 1 autoreerrssioon here, though extensions that condition on more than just one previous timestep are possible.
Priors. The Appendix describes prior distributions p(0) on parameters assumed for the purpose of regular- ization. We use "sticky" Dirichlet prior (Fox et al., 2011) to obtain smoother system-level segmentations.
Specification. To apply HSRDM to a concrete problem, a user must select the number of system states L and entity states K as well as functional forms of g. We assume that g can be evaluated in O(J).
Special cases. If we remove the top-level system states $0.7 (or equivalently set L 1), our HSRDM reduces to separate models for each of J entities, where each per-entity model is recurrent autoregressive HMM (rAR-HMM) as in Linderman et al. (2017). Tf we remove the recurrence from our HSRDM, we obtain a multi-level autoregressive HMM that we refer to as a hierarchical switching dynamical model (HSDM). Later on, we compare our model to both ablations in several experiments.
3 Inference
Given observed time scries ForT we now explain how to simultancously estimate parameters 0 and infer approximate posteriors over hidden states sp.r for the system and hidden states *o. Because all system-level and entity-level states are unobserved, the marginal likelihood p(Zur ()) a natural objective for parameter estimation. However, exact computation of this quantity, by marginalizing over all
(1:J)
(1:J) for all J entities.
4

---

Published in Transactions on Machine Learning Research (08/2025)
hidden states, is intractable. Given L system-level states and K entity-level states, computing p( |0) naively via the sum rule requires a sum over (I.KJ)² values. While the forward algorithm (Rabiner, 1989) resolves the exponential dependence in time, the exponential dependence in the number of entities persists: TLK21 operations are required to do forward-backward on HSRDMs. This exponential dependence remains prohibitively costly even in moderate settings; for instance, when (T,J, K) = (100, 10, 4), a direct. application of the forward algorithm requires around 220 trillion operations.
Instead, we will pursue structured approximation g to the true (intractable) posterior over hidden states. Following previous work Slameda-Pinedae et al., 2021; Linderman et al., 2017), we define
( #0:T g(5o:T) Q(=0.r)).
(3.1)
intending ((Surr 20.T ) p(surr, PoTT o:T ,0). Each factor of q is paramcterized separately, without dependence on 0 or other random variables in the model. Each factor retains temporal dependency structure,
(1:3)
(1:))
avoiding the simplistic independence assumptions of complete mean-field inference (Barber et al., 2011).
Using this approximate posterior q, we can form a variational lower bound on the marginal log likelihood (1:J) ,(1:J)
VLBO log p(x(2r) (1:J) 0), defined VLBO [6, g = Eg log p( 001 -0:T Su-T | [6]] [q(=0))
(1:J), -oor))] As shown in the Appendix, computation of this bound scales as O(TJLL crucially linear rather than exponential in the number of entities J. This reduces the approximate number of operations required for inference on the earlier moderate example setting (T'=100, J=10, 2 K=4) from 220 trillion to 64 thousand.
To estimate A and q given data Ty:T; we pursue coordinate ascent variational inference (CAVI; (Blei et al., 2017)) on the VLBO, known as vaniational expectation maximization (Beal, 2003) when is approximated with a point mass. Given a suitable initialization, we alternate between specialized update steps to each variational posterior parameter:
{
q(SO.T) x exp
y ( log (( (17)) F0r sO. 01]
(1:J) (1:))
(3.2)
q (1:1)) x exp
E Eotsorr pl (o ooT SO:T {1]]
(1:J) (1:))
log logp(2007 (1)) (1:J) 2T. s0:T 0) log p(0)
0 = arg max
E (1:
Dat )a(sor)) 90-
The updates above define the variational E-S step (VES step), variational E-Z step (VEZ step), and M-step, respectively. The first two formulas in equation 3.2 are derived by following the well-known generic variational recipe for optimal updates (Blei al., 2017). The M-step allows the inclusion of an optional prior on some all parameters. We've worked out efficient ways to achieve the optimal update for each step, described below. Full details about each step, well as recommendations for initialization, are in the Appendix. We also share code via the link at bottom of page 1, built upon JAX for efficient automatic differentiation (Bradbury et al., 2018).
VES step for system-level state postcriors. We can show the VES step reduces to updating the posterior of a surrogate Hidden Markov Model with J independent autoregressive categorical emissions. Optimal variational parameters for this posterior can be computed via a that extends classic forward-backwart for an AR-HMM to handle recurrence. The runtime required is O[T + KD KL + KM) + +T?)).
algorithm
VEZ step for entity-level state posteriors. We can show that the VEZ update reduces to updating the posterior of separate surrogate Hidden Markov Models for each entity j with autoregressive categorical emissions which recurrently feedback into the transitions. Given a fixed system-level factor q(su:r), we can update the state posterior for entity j independently of all other entities. This means inference linear in the number of entities J, despite the fact that the HSRDM couples entities via the system-level sequence. The linearity arises even though our assumed mean-field variational family of Equation 3.1 did not make an outright assumption that q(<(1)2)) = II5- ((2(7)) Optimal variational parameters for this postcrior can again be computed by dynamic programming that extends the forward- backwardd algorithm. The runtime required to update each entity's 's factor is [ 2 KD2 KT+ +KM).
5

---

Published in Transactions on Machine Learning Research (08/2025)
M step for transition/emission parameters. Updates to some parameters, particularly for emission model parameters when H has exponential family structure (such as the Gaussian or Von-Mises AR likelihoods we use throughout experiments), can be done in closed-form. Otherwise, in general, we optimize 0 by gradient ascent on the VLBO objective. This gradient ascent has the same per-itcration cost as the computation of the the recurrence function parameters v or Φ are learnable, these can also
VLBO, with runtime O(T). If be updated in the M step.
Like many variational methods, the alternating update algorithm in equation 3.2 provides useful guarantee. Each successive step will improve the VLBO objective (or if incorporating priors, the modified objective VLBO + p(f) until convergence to local maximum. This improvement assumes any M step that uses gradient ascent relies on a suitable implementation that guaranlees a non-decrease in utility.
4 Related Work
Below we review several threads of the scientific literature in order to situate our work.
Continuous reprosentations of individual sequences. Other efforts focus on latent continuous repre- sentations of individual time series. These can produce competitive predictions, but do not share our goal of providing a segmentation at the system and entity level into distinct and interpretable discrete regimes. Probabilistic models with continuous latent state representations are often based on classic linear dynamical system (LDS) models (Shumway & Stoffer, 1982). Deep generative models like the Deep Markov Model (Krishnan et al., 2017) and DeepState (Rangapuram et al., 2018) exlend the LDS approach with more flexible transitions or emissions via neural networks.
Discrete state representations of individual sequences. Our focus on discrete state representations which provide interpretable segmentations of available data, line of work that started with classic approaches to entity-level-only sequence models like hidden Markov models (HMM) (Rabiner, 1989) and autoregressive hidden Markov models (AR-HMM) (Ghahramani & Hinton, 2000), later extended to include both discrete and continuous latent representation ationss as in switching-state linear dynamical systems (SLDS) ( Slameda-Pinedde et al., 2021). Recent efforts such as DSARF (Farnoosh et al., 2021) and DS3M (Xu et al., 2025) have extended such base models to non-linear transitions and emissions via neural networks. All of these efforts still represent each time series with only entity-level (not system-level) dynamics, and do not incorporate recurrent feedback to guide the evolution of latent variables.
Discrete states via recurrence on continuous observations. Linderman et al. (2017) add a notion of recurrence to classic AR-HMM and SLDS models, increasing the flexibility in each timestep's transition distribution by allowing dependence on the previous continuous features, not just the previous discrete states. Later work has extended recurrence idcas in several dircctions that improve entity-level sequence modeling, such as multi-scale transition dependencies via the tree-structured construction of the T-SLDS (Nassar el al., 2019), recurrent transition models via SNLDS (Dong al., 2020), or recurrent transition models that can explicitly model state durations via RED-SDS (Ansari et al., 2021). To model multiple recordings of worm neural activity, Linderman et al. (2019) pursue recurrent state space models that are described as hierarchical because they encourage similarity between each worm entity's s custom dynamics model via common parameter priors in hierarchical Bayesian fashion. Their model assumes only entity-level discrete states.
Multi-level discrete representations. Stanculescu et al. (2014) developed a hierarchical switching linear dynamical system (HSLDS) for modeling the vital sign trajectories of individual infants in an intensive care unit. The root level of their directed graphical model assumes discrete state sequence (analogous to our s) indicating whether disease was present or absent in the individual over time, while lower level discrete states (analogous to our 2) indicate the occurrence of specific "factors" representing clinical events such as brachycardia or desaturation. While their graphical model also contains a multi-level discrete structure, we emphasize three key differences. First, they require fully-supervised data for training, where each timestep is labeled with top-level and factor-level states. In contrast, our structured VI routines to simultaneoull estimate parameters and hidden states in the unsupervised sctting are new. Second, their model does not incorporate recurrent feedback from continuous observations. Finally, they model individual time series not multiple interacting entities.
6

---

Published in Transactions on Machine Learning Research (08/2025)
More recently, hierarchical time series models composed of Recurrent Neural Networks (RNNs) have been proposed for dynamical systems reconstitation (Brenner et al., 2025). Different from our work, this framework is not designed for modeling entity interactions within a single system, but modeling shared properties among multi-domain dynamical systems for time series transfer learning.
Tastly, Hierarchical Hidden Markov Models (HHMMs) (Fine et al., 1998) and their extensions (Bui et al., 2004; Heller et al., 2009) describe single entity's observed sequence with multiple levels of hidden states. The chief motivation of the HHMM is to model different temporal length scales of dependency within an individual sequence. While HHMMs have been applied widely to applications like text analysis (Skounakis et al., 2003) or human bchavior understanding (Nguyen et al., 2005), to our knowledge HHMMs have not been used to coordinate multiple entities overlapping in time.
Models of teams in sports analytics. Terner & Franks (2021) survey approaches to player-level and team-level models in basketball. Miller & Bornn (2017) apply topic models to tracking data to discover how low-level actions (e.g. run-lo-basket) might co-occur among leammates during the same play. Metulini et al. (2018) model the convex hull formed by the court positions of the 5-player team throughout a possession via one system-level hidden Markov model. In contrast, our work provides a coordinated two-level segmentation reprosenting the system as well as individuals.
Personalized models. Several switching state models assume each sequence in a collection have unique or personalized parameters, such custom transition probabilities or emission distributions (Severson et al., 2020; Alaa & van der Schaar, 2019; Fox et al., 2014). In this style of work, entity time series may be collectod asynchronously, and entities related by shared priors on their parameters. In contrast, we focus entities that are synchronous in the same environment, and relate entities directly via system-level discrete chain that modifies entity-level state transitions.
Models of coordinated entities. Several recent methods do jointly model multiple interacting entities or "agents", often using sophisticated neural architectures. Zhan et al. (2019) develop a variational RNN where trajectories are coordinated in short time intervals via entity-specific latent variables called "macro- intents". Yuan et al. (2021) develop the AgentFormer, a transformer-inspirdd stochastic multi-agent trajectory prediction model. Alcorn & Nguyen (2021) develop baller2vec++ a transformer specifically designed to capture correlations among basketball player trajectories. Xu et al. (2022) introduce GroupNet to capture pairwise and group-wise interactions. Unlike these approaches, ours builds upon switching-state models with closed-form posterior inference and produces discrete segmentations. Our approach may also be more sample efficient for applications with only a few minutes of data, as in Sec. 5.4.
Models of interaction graphs. Some works (Wu et al., 2020; Löwe et al., 2022) seek to learn an interaction graph from many entity-level time series. In such a graph, nodes correspond to entities and edge existence represents a direct, pairwise interaction between entities. Other works (Kipf et al., 2018; Webb et al., 2019) assume collecee graph, but can learn to annotate each edge with different discrete type representing different kinds of interaction. One possible type may be hard-coded to mean "non-edge" for no interaction. Intentional priors can control the graph sparsity (frequency of non-edge labels). Recently, the GRAph Switch- ing dynamical Systems (GRASS) approach (Liu et al., 2023) models interactions via a latent graph whose edges change dynamically over time.
For some applications, discovering possible pairwise interactions or interaction types that influence data is explicit goal. In our chosen applications (e.g. basketball player movements). domain expertise indicates the graph is fully-connected. Moreover, when the graph is fully connected, interaction graph approaches burdensomely require runtimes that are quadratic in the number of entities J. In contrast, our approach models system-level dynamics explicitly with a more affordable runtime cost that is linear in J. Very recent work by Wang & Pang (2024) offers a different route to scaling to J in the hundreds. Their approach can estimate directed graph structure by combining a variational dynamics encoder with partial correlation ideas.
5 Experiments
We now demonstrate our model's utility across several experiments. As we aim to highlight our model across two tasks- multi-step-ahead forecasting and interpretability of system dynamics- we show one of each on
7

---

Published in Transactions on Machine Learning Research (08/2025)
synthetic and real datasets. In all but the first task, the data generation model is either unknown not the same as the HSRDM. Within each task, we compare our model's performance to task-specific competitive baselines. We also compare to ablations that remove either the top-down influence of system-level hidden states or bottom-up recurrent feedback from observations. Overall, we see that our compact and efficient HSRDM is able to oulperform alternatives with respect to discovering hidden system dynamics and maintain similar or better prediction performance than computationally intensive neural network methods.
5.1 FigureEight: Synthetic task of forecasting coordinated 2D trajectories
To illustrate the potential of our HSRDM for the purpose of high-quality forecasting of several coordinated entitics, we study a synthotic dataset we call FigureEight. In the truc generative process (dotailed in App. D.3), each entity switches welween clockwise motion around a top loop and counter-clockwin motion around a bottom loop. The overall observed entity-level 2D spatial trajectory over time approximates the shape of an "8" as in Fig. 2. Each entity has two true states, one for cach loop, with specific Gaussian vector autoregression process for each. Transition between these loop states depends on both top-down and bottom-up up signals in the data-generating process. For top-down, a binary system-level state sets which loop is favored for all entities at the moment. For bottom-up, switches between loop states are only probable when the entity's 's current position is near the origin, where the loops intersect. Though coordinated, entity trajectories are not perfectly synchronized, varying due to individual rotation speeds and initial positions.
We generate data for 3 entities over 400 total timesteps. For training, we provide complete data for the first two entities (times 1-400) and partial data for the last entity (times 1-280). The prediction task of interest is aatity-specific partial forecasting: estimate the remaining trajectory of entity j = 3 for times 281-400 (120 time steps), given partial information (timesteps 1-280) for that entity and full information (1-400) for other entities. The true trajectory for this heldout window is illustrated in Fig. 2: we see a smooth transition over time from the top loop to the bottom loop of the '8', We wish to compare our method to competitors at estimating this heldout trajectory.
Baseline selection. To show the benefits of modeling system-level dynamics, we compare to baseline methods that can produce high-quality forecasts with discrete latent variables, yet only model entity-level (not system-level) dynamics and can do partial forecasts. First, we select the deep switching autoregressioe factorization model (DSARF; (Farnoosh et al., 2021)). This model represents recent state-of- state-of-the-art forecasting performance and uses deep neural networks to flexibly define transition and emission structures yet still can produce discrete segmentations. Second, we compare to a recurrent autoregressive hidden Markov model (HAR-HMM; (Linderman et al., 2017)). This is an ablation of our method that removes our system-entity hierarchy. This experiment requires partial forecasting of one entity at times 281-400 given context from other entities. Some neural net baselines like Agentformer and GroupNet considered in later experiments do not easily handle such partial forecasting in released code, so exclude them here on this task.
Entity-to-system strategy. Since DSARF and 1AR-HMM baselines each only model entity-level dynamics, for cach one we try three different entity-to-syytemm strategies to convert any entity-only model to handle a system of entities. First, complete Independence fits a separate model to each entily's data only, with no information flow between entities. Next, complete Pooling fits a single model on N' N total sequences, treating each sequence n from each entity j as an i.i.d. observation. Finally, Concatentation models a multivariate time series of expanded dimension D' = J. D constructed by stacking up all entity-specific feature vectors .r(")) ... ......")) at each time t.
Method configuration. For the HSRDM, we set L= 2 system states and K 2 model, we pick Gaussian autoregressive to match the true process. We do not use any system-level recurrence 9 as the true data generating process does not have this feedback. We set entity-level recurrence f to radial basis function indicating distance from the origin. While we do not expect this f to improve training fit, we do expect. it to improve forecasting, as it capt.res a key aspect of the true process: that switches between loops are only probable near the origin.
2 entity states. For emission
For all baselines (DSARF and HARHHMM), we set the number of entity states at K to match the intended ground truth. The FAR-HMM uses the same emission model and entity-level recurrence f our HSRDM.
8

---

Published in Transactions on Machine Learning Research (08/2025)
Ground Truth
HSRDM
DSARF
15
15
42
Independont
Pooling
Concatantioon
-302
Model
MSE
Avg. MSE
seeple,
Mda. MSE
wnaphe
TAR-HMM
HSRDM
0,0025
0.105 (0.002)
0.01
-307
No Recurnntt Feraltack
0.5
2.41 (0.78)
2.41
Tireeton
340
EAR-HIM Independent
0.05
2,08 (0.55)
2,08
1AM-IMM Pooling
L,7 (0.19)
2,56
-303
AAR-HMM Concalatation
(0.3)
1.19
DSARF Independens
031
L.73 (0:27)
2.12
35
110
60
DSARI Pooling
0.91
189 (0.19)
1.91
lidependeit
Pooling
Concatonation
DSARF Conestenation
115
4.4((1.50)
3.01
Figure 2: Predictions for the heldout time segment of one entity in Fight task. Models were trained all data (1-400) from entities 1-2 and times 1-280 for entity 3, then asked to forecast the heldout period (times 281-400) for entity 3. Colors correspond to individual time steps, shown in the bar the bottom left. X-marks in each panel show the ground truth trajectory; solid dols show the model forecast. Top left: Our best HSRDM prediction closely matches the ground truth at all time points. Top right: Best overall sample for the DSARF baseline under each strategy. Bottom right: Best overall sample for the T ARR-MMM baseline, under each strategy (Indep., Pool, and Concat., defined in Sec. 5.1) for adapting an entity-only model to our hierarchical setting. rAR-HMM Independent is "no system state" ablation of our model. Bottom left: Table reports for each model: best mean- setet--uuuuddrrrreee (MSE) across independent trials and forecast samples, average MSE (with standard error in parens) across the samples from the trial with the best MSE, and median MSE across sample averages from all trials. This table includes ablation results for version of our HSRDM with no recurrent feedback. forecasting plot can be found App. D.3.
Hyperparameter tuning. Optimal hyperparameters for each model and entity-to-systed strategy are determined independently. For DSARF, we tune its number of spatial factors and the lags indicating how the next timestep depends on the past. We don't tune any specific hyperparameters of our HSRDM or the AAR-HMM.
Training. For each tested method (where method mcans model and (if needed) entity-to-systen strategy), at each hyperparameter train via 5 separate trials with different random seeds implying distinct inilializations of parameters. This helps avoid local optima common to models with latent discrete state. DSARF is intentionally allowod more trials (10) to stress fair ovaluation to external basclincs. All DSARF models were trained with 500 epochs and a learning rate of 0.01. The HSRDM and its ablations are trained with 10 CAVI iterations. All models required similar training time on this small dataset: < minute per run. Reproducible details for all methods (including specific hyperparameters) are in the Appendix and released code.
Forecasting. To forecast, from each trial's fit model draw 5 samples of the heldout trajectory 2281-110 for target entity j = 3 over the time period of interest. We keep the "best" sample, meaning the sample with lowest mean squared error (MSE) compared to the true (withheld) trajectory. Each method can be represented by the trial and hyperparameter setting with lowest best-sample MSE, or via summary statistics across samples, trials, and hyperparameters.
Results: Quantitative error. For each method, we report the best-sample MSE in Fig. 2 (bottom left). Our HSRDM outperforms others by wide margin, scoring best-sample MSE less than 0.003 compared to the next-best value of 0.05 and values of 0.31 1.75 among others. To indicate reliability across forecast samples, we further report the average sample MSE from the trial that produces the best sample (with standard errors to indicate variation). This average is useful to distinguish between a trial that consistently produces quality forecasts from a trial that "gets lucky." While the TAR-HMM Independent model achieves a relatively low best MSE, the average MSE across all 5 samples from that trial is much worse than our method (2.08 vs. 0.105). Note that the FAR-HM Independent model is a "no system state" ablation of the regular HSRDM. Our approach also outperforms a "no recurrent feedback" ablation (qualitative results in App. D.3).
9

---

Published in Transactions on Machine Learning Research (08/2025)
To understand whether MSE is consistent across trials, we also report the median across trials of the average- sample MSF. For our HSRDM, the median-over-trials MSE is quite low (0.01), indicating that the HSRDM is typically reliable. The HSRDM's best-sample trial happened to also yield one of the worst samples, which is why that trial's average-sample MSE is higher than the median. Comparing the best column to the median column, all methods show wide gaps indicating variation in quality across trials. There a clear need to avoid local optima by considering many random initializations.
Visual quality. Fig. 2 also shows the best sampled trajectory from each method. The forecast from our proposed HSRDM looks quite similar to the true trajectory. In contrast, every competitor struggles to reproduce the truth. The closest competitor method is the rAR-HMM with the Independent entity-to-ssstem strategy. The main difference lies in how that model produced a larger outer circle for the bottom loop.
Results:
5.2 NBA Basketball: Real task of forecasting 2D player trajectories
We next aim to evaluate the HSRDM's lorecasting capabilities on real movements of professional basket team. Specifically, we model the 5 players of the NBA's Cleveland Cavaliers (CLE), together with their 5 opponents, across multiple games in an open-access dataset (Linou, 2016) of player positions over time recorded from CLE's 2015-2016 championship season. To better evaluate the ability to model specific entities, we focused exclusively on 29 games involving one of CLE's four most common starting lineups. We randomly assigned these games to training (20 games), validation (4 games), and test (5 games) sets.
We split each game into non-overlapping basketball event segments, typically lasting 20 seconds to 3 minutes. Event segments contain periods of uninterrupted play (e.g. shot block - rebound offense + shot made) from the raw data, ending when there is an abrupt break in player motion or a sampling interval longer than the nominal sampling rate. Each event segment gives 2D court positions over time for all 10 players that form multi-entity emission sequence (1:J), Each such sequence is modeled as an independent draw from our proposed HSRDM or competitor models. We standardized the court so that CLE's offense always faces the same direction (left), and downsampled the data to 5 Hz.
Baseline selection. We compare the forecasting performance of the HSRDM to multiple competitors. Because of DSARF's poor prediction performance on even our synthetic task, here we instead explore SNLDS (Dong et al., 2020) as a neural network baseline which can provide flexible model of complex time series. Because the SNLDS baseline is restricted to entity-level (not system-level) dynamics, we fit an independent SNLDS model to cach player. To exemplify neural network methods that can prodict trajectorios of groups dircctly, we select GroupNet (Xu et al., 2022) and Agentformer (Yuan el al., 2021), given GroupNet and Agentformer's leading reputations for high-quality forecasting of multiple entities (Xu et al., 2022; Yuan et al., 2021).
We further consider two ablations of our method. First, removing system- level states yiclds an independent rAR-HMM HMM (Linderman et al., 2017) for each player. Second, we remove recurrent feedback. As in Yeh et al. (2019), we also try simple but often competitive fized velocity baseline.
Model configuration. The ground truth number of states is unknown; we pick K = 10 entity states and L = 5 system states. For our HSRDM and its ablations, our emissions distribution is a Gaussian vector autoregression with entity-stato-dependent parameters (sec Soc. E). System-level recurrence g reports all player locations at (I:J) to the sestem-level transition function, allowing future latent states to depend on player locations. Following Linderman et al. (2017), our entity-level recurrence function f reports an individual player's location a!") and out-of-bounds indicators to that player's 's entity-level transition function, allowing each player's next state probability to vary over the 2D court. We use a sticky prior for system-level state transitions with a 1 and 50 (see Sec. E.2 for details).
Hyperparameters. Optimal hyperparameters for each model and strategy are determined independently. For GroupNet, we used hyperparameters recommended by Xu et al. (2022) from their own basketball models. For Agentformer, we use the recommended architecture and training settings from Yuan et al. (2021), adjusting learning rates and number of epochs for our data (see App. E). For SNLDS, we use the architocture and training sottings from Ansari et al. (2021) model an electricity datasot, adjusting certain hyperparameters to match the HSRDM well (see App. E). We don't tune any hyperparameters of our HSRDM or its ablations. Ablations inherit hyperparameters like number of states where applicable.
10

---

Published in Transactions on Machine Learning Research (08/2025)
Table 1: Quantitative evaluation on Basketball.
20
Grouplict
11.6 (0.5)
12.2 (0.4) /
12.4 (0.3) /
Directional Varéation
X Tn Rowads
HSRIX (ours)
15.5 (-1
13.8 [ [-]
13:4 (-)
HSRDN (ours)
.449
.54
No rystem state (raumees)
15.9 (0.3)
15.4 (0:4) X
15.9 [0.3)
X
No system state (rAREWNs)
.631 (.013)
X
906 (.008)
X
No fredback (HSDN)
X
X
16.0 (0.3)
No recurrent feedback (ISDN)
.499 (.017)
881 (.008)
X
SuLDEa
17.2 (0.7) X
15.0 (0.5) X
16.1 [0.5]
Pised velocity
17.0 X
17.0 (0.5)
X
17.0 (0.5)
Agentforzer
35.1 (0.4)
X
21.3 (0.7)
X
25.9 (0:4)
(b) Statistical comparisons to ablations.
(a) Forecast error (in feet)) vs. train set size (num. games 12).
Lebron James
Kevin Love
T. Mozgov
J. Smith
M. Williams
GroupNet
HSRDM
No system state
4
4
(*ARHMMs)
No recurrent feedback (HSDM)
SNLDSs
Agentformer
Figure Sample forecasts of NBA player location trajectories. Grey tracks show true player trajectories from the 6-second forecasting window of the test event on which our model had median forecasting error. In color are three sampled forecasts from each model; purple/red/gree indicate 1st/5th/Iue 75th/10th best forecasting error (from 20 samples). Time runs from light to dark.
Training time. Our HSRDM is more computationally efficient for this task. Training HSRDM 2023 Macbook with Apple M2 Pro chip Training Agentformer on an Intel Xeon Gold 6226R CPU took 1.5, 6, and 13 hours, respectively. That is, our HSRDM was 17-45 times faster to train. We saw similar gains with the other neural network baselines; on the same 2023 Macbook, HSRDM trained 120 times faster than GroupNet and 10-50 times faster than SNLDS.
1, 5, 20 training games took 2, 15, and 45 minutes, respectively.
Model size. Our compact HSRDM has 9,930 parameters. In contrast, GroupNet has over 3.2 million parameters (320x larger) and Agentf ormer has over 6.5 million (650x larger). The collection of SNLDS models has about 0.5 million parameters (50x larger), even though this approach does not model cross-entity interactions.
Evaluation procedure. We randomly select 6 second forecasting window within each of the 75 test set events. Preceding observations in the event are taken as context, and postceding observations are discarded. We sample 20 forecasts from each method for the 5 starting players on the Cavaliers. We report the mean distance in feet. from forecasts to ground truth, with the mean taken over all events, samples, players, timesteps, and dimensions. We perform paired t-tests on the per-event differences mcan distances between our model vs competitors, using Benjamini & Hochberg (1995)'s correction for multiple comparisons over positively correlated tests. All tests were performed at the .05 significance level for two-sided tests.
Results: Quantitative crror. Tab. la reports the mean distance in feet from forecasts to ground truth for each method. Methods whose forecasting error are siguificantly better (worse) than HSRDM according to the hypothesis tests are marked with a green check (red x). The standard error of the difference in means is given in parentheses. Our first key finding that HSRDM provides better forecasts than ablations, supporting the utility of incorporating multi-level recurrent feedback and top-level system states into switching autoregressive models of collective behavior. Second, HSRDM provides better forecasts than some neural network baselines
11

---

Published in Transactions on Machine Learning Research (08/2025)
(SNLDSs and Agentformer), as well as the fixed velocity baseline. Third, a neural network method designed to model coordinated entities (GroupNet) does provide the best overall forecasts in terms of error alone. Yet. we consider our results "competitive" in that on basketball court measuring 94x50 feet, our n = 20 forecasts are on average only one foot further off in terms of error while using a simple interpretable model that is roughly 320x smaller and 120x faster to train.
Results: Qualitative forecasts. Fig. shows sampled forecasts from our model and baselines. As a first key finding, HSRDM forecasts are qualitatively more similar to GroupNet than the forecasts of other baselines. Second, system-level switches appear to help coordinate entities; players move in more coherent directions under HSRDM than the no system state ablation. Third, our multi-level recurrent feedback supports location-dependenl state transitions; players are more likely to move to feasible in-bounds locations under HSRDM than without recurrence.
Results: Statistical comparisons. To corroborate the above conclusions from visual inspection, we examine Iwo statisties on the entire test. set: Directional Varialion, which measures the coherence of movements by basketball team via the variance across players of the movement direction on the unit circle between the first and last timesteps in the forecasting window, and % In Bounds, the mean percentage of each forecast that is in bounds. We report these statistics for the HSRDM and its ablations in Tab. 1b. Methods significantly different from the HSRDM baseline according to hypothesis testing are marked with a red "x." These results corroborate the intended purpose of each of our modeling contributions. Removing the system-level states significantly increases the directional variation across players, indicating lack of top-down coordination. Removing recurrence significantly reduces the percentage of forecasts that remain in bounds. Removing system states also significantly reduces the in bounds percentage, although to lesser extent. We suggest that greater coordination in the player movement also helps keep the players in bounds.
5.3 MarchingBand: Synthetic task of interpreting marching band routines
Beyond the forecasting evaluations of carlier subsections, we now evaluate the HSRDM's capabilities to discover useful and interpretable discrete system dynamics in systems with many (J 64) entities. To do this, we introduce synthetic dataset, MarchingBand, consisting of individual marching band players ("entities") moving across a 2D field in a coordinated routine to visually form sequence of letters. Each observation is a position 2.(5) IR22 of player j a time within the unit square centered at (0.5, 0.5) representing the field.
By design, player movement unfolds over time governed by both top-down and bottom-up signals. The team's goal is to spell out the word 'LAUGH'. An overall discrete system state sends the lop-down signal of which particular letter, one of 'L', 'A', "U', "G", or 'H", the players form on the field via their coordinated movements. Each entity has an assigned vertical position on the field (when viewed from above), and moves horizontally back and forth to "fill in" the shape of the current letter. Example frames are shown in Fig. 4. Each state is stable for 200 timesteps before transitioning in order to the next state.
Each entity's position over time on the unit square follows the current letter's top-down prescribed movement pattern perturbed by small-scale i.i.d. zero-mean Gaussian noise. When reaching field boundary, typically the player is reflected back in bounds instantaneously. However, there is small chance an entily will continue out-of-bounds (OOB, In1 & (0, 1)). When enough players go OOB, this bottom-up signal triggers the system state immediately to a special "come together and reset" state, denoted 'C'. This reset state does not form the letter "C". Instead, in state "C" all players move to the center for the next 50 timesteps, then return to repeat the most recent letter before continuing on to remaining letters.
Altogether, we build a dataset of N 10 independently sampled sequences of 'LAUGH' each with J 64 marching band players. Each sequence contains a different mumber of time-steps, ranging from 1000 to 1100, depending on how many reset segments occur. To trigger reset state "C", we use a threshold of 11 players OOB as this generates a moderate number of 6 "C" segments distributed across the 10 sequences.
Research goal. We seek to understand whether our proposed model or competitors can discover the overall discrete system-level dynamics, as measured by discrete segmentation quality. The true data-generating process, while similar to the HSRDM in using top-down and bottom-up signals, is not strictly an HSRDM generative model because each state has fixed duration. Thus, all methods are somewhat. misspecified here.
12

---

Published in Transactions on Machine Learning Research (08/2025)
L
G
DSARF Concatenation
2000
ACOO
0000
2000
22000
Ground Truth
HSRDM No System State (FAR-HMM)
2000
5002
5000
HSRDM
HSRDM No Recurrent Feedback
100.00
Figure 4: System-level segments for all models the MardhingBand task. Top left: An example snapshot in time of ground truth entity observations for each system-state letter (spells "LAUGH"). Bottom left: The ground truth segmentation across all 10 sequences well as the HSRDM segmenlation. Right side: Results for ablations (no system state (rAR-HMM), no recurrent leedback) as well as DSARF.
Baseline selection. We focus on methods that can produce an estimated discrete segmentation at the system-level. Using this criteria, our primary external competitor is DSARF (Farnoosh et al., 2021). Neither GroupNet nor Agentformer can produce segmentation or offer interpretable discrete system states.
Additionally, compare to two ablations of our method: removing the system-level state (leaving an FAR-HMM) and removing the recurrent feedback (leaving a multi-level HMM). For both DSARF and FAR-HMM, we obtain system-level segmentations in two steps. First, we use the Independent strategy ( SSec. 5.1) to obtain entity-level segmentations. Next, for each timestep we concatenate the one-hot indicator vectors of each entity to form a longer vector of size J x K. These per-time features are clustered via k-means to obtain a system-level segmentation. This should favor clustering timesteps that share entity state assignments.
Model configuration. For all methods, we fairly provide knowledge of the true number of system states, 6= KL, A, U, G, H, C)I. Models with system states (HSRDM and its recurrent ablation) have L = 6 system states and K = 4 entity states. Models with only entity-level states (DSARF and rAR-HMM) are fit with K = 6 entity states, imagining one state per letter, followed by k-means to discover 6 system-level states. For HSRDM, we set system-level recurrence g to count the number of entities out of bounds. For HSRDM and FAR-HMM, we set entity-level recurrence f to the identity function. The emission model set to a Gaussian Autoregressive.
Hyperparameter tuning. For DSARF, we again select the number of spatial factors and the lags, via grid search (see Appendix F. F.4). Ultimately, select 25 spatial factors and the set of lags {1, (1.... 200), suggesting long-range dependency is useful to compensate for lack of recurrence. For our HSRDM and ablations, no specific hyperparameter search was done. Methods were trained with similar epochs/iteratie /iiteaateiinnnss as the Eight data, with reasonable convergence verified by trace plots. To avoid local optima, for each method we take the best trial (in terms of segmentation accuracy) of 5 possible seeds controlling initialization.
Results: segmentation quality. Fig. 4 al- lows visual comparison of the ground truth system-level segmentation and the estimated segmentations from various methods. Esti- mated states are aligned to truth by minimizing classification accuracy (Hamming distance of one-hot indicators) via the Munkres algorithm (Munkres, 1957). The visualized segmentations depict each model's best classification accuracy
Table 2: Accuracy of system-level segmentation on Morching- Band, compared to ground-truth after alignment.
Method
Best Acc.
Median Acc.
HSRDM
88%
83%
No Recur. Feedback (HSDM)
80%
69%
No System State (НAR-HMM)
48%
42%
DSARF
38%
32%
13

---

Published in Transactions on Machine Learning Research (08/2025)
over multiple training trials with different random initializations. DSARF was given 10 chances and all others 5 in an attempt to provide DSARF with a more-than-fair chance. As reported in Tab. 2, our HSRDM obtains 88% accuracy overall, with clear recovery of each of the 6 states in visuals. In contrast, other methods struggle, delivering 38- - 80% accuracy and notably worse segmentations in Fig. 4. Even with respect to median accuracy across trials (also in Tab. 2), the HSRDM outperforms its competitors. Inspecting Fig. 4, the HSRDM recovers each true state most of the time. In contrast, the next-best method, the no recurrence ablation, completely misses capturing the "U" system-level state across examples even in its best trial. This experiment highlights the HSRDM as a natural model to recover hidden system dynamics when many individuals influence collective behavior. We provide supplemental results for accurate estimation of system states for a larger number of entities (J= = 200) and a larger number of dimensions (D 303) in the Appendix F.3.
5.4 Soldier Training: Real task of interpreting risk mitigation strategies
For our final experiment, we model a squad of active-duty U.S. Army soldiers during training excrciso designed to improve leam coordination. Our goal is to demonstrate the HSRDM's capability to reveal useful and interpretable group dynamics in a real application. For simplicity and to highlight our model's capability as an explanatory tool, we focus on interpreting one model's fit rather than comparisons to alternatives.
In this training exercise, the squad's task to maintain visual security of their entire perimeter while simulated enemy fire comes primarily from the south. Focus on the south creates potential blindside to the north. If this blindside is left unchecked for sufficiently long time, this leaves the squad vulnerable to a blindside attack. Among several goals, the squad was instructed that a primary goal was mitigating overall risk with visual security as a key sub-task. As a way to mitigate risk, the squad should periodically plan to have at least one soldier briefly turn their head to the north to regain visibility and reduce vulnerability to a blindside atlack. Strong performance at this sub-task requires coordination ion across all soldiers the squad. Our modeling aim is to utilize the HSRDM to interpret the soldiers' risk mitigation strategy with respect to checking their blindside they complete the overall task.
The data consists of univariate time series of heading direction angles xx{") recorded at 130 Hz from each soldier's helmet inertial measurement unit (IMU), downsampled to 6.5 Hz. We have one 12 minute recording of one squad of 8 soldiers. Raw data from the first minute of contact is illustrated in Fig. 5. Due to privacy concerns, the dataset is not shareable. This study was approved by the U.S. Army Combat Capabilities Development Command Armaments Center Institutional Review Board and the Army Human Rescarch Protections Office (Protocol # 18-0033.
Model configuration. Our HSRDM captures the risk mitigation strategies of the group by setting the system-level recurrence function g to the normalized clapsed time since any one of the J soldiers looked within the north quadrant of the circle. No entily-level recurrence I was included. Soldier headingg must remain on the unit circle throughout time, so we use a Von Mises antoregression as the emission model Hç for the k-th state of the j-th soldier:
a000) 1
22 1 t2t (2())) =k) ~ VM (3 kg,,. where ((,((((2))) = = 2221 + dj
(5.1)
The Von Mises distribution (Banerjee et al., 2005; Fisher & Lee, 1994), denoted VM k), is a exponential family distribution over angles on the unit circle, governed by mean μ and concentration λ > 0. Here, a,,x is an autoregressive coefficient, 65,* is a drift term, and Ki..k a concentration for entity j in state k.
We set the number of entity- and system-level states to K = 4 and L = 3, based on a quick exploratory analysis. We use sticky Dirichlet prior for system-level transitions, as in Equation B.2, with k == 50.0, so that the prior puts most probability mass on self-transition probabilities between .90 and .99.
1.0 and
Training. We applied the CAVI inference from Sec. 3, observing suitable convergence after 10 iterations.
Results. Fig. 5 visualizes key results. Inspection of the inferred system-level states, entity-level states, and Iearned transition probabilitics suggests that the model learns a special risk mitigation strategy for Soldier 6. Specifically, the model learns a "turn north" state (blue) for Soldier 6 that is particularly probable when the entire squad reaches system-level "high elapsed time since any blindside check" state (red).
14

---

Published in Transactions on Machine Learning Research (08/2025)
(N)
0.30
135
45
0.35
180
(9)
Orean
vodio
Cen Arber hol "ur hoh
0.54
0.26
0.96
mot
WAy
0:32
0.1
M.M4
0.21
0:33
0.59
09:185 100000 09:1915 09:1922 06:18:33
06:1651
Figure 5: Modeling the heading directions of squad of soldiers engaged én simalated battle. Top left: Heading directions (in degrees) of squad of soldiers over time. Each color represents different soldier. Time moves from center of circle to boundary. Top right: Inferred squad-level states So.7 (colors) superimposed over black curve representing the squad's cumulative security risk the north direction (elapsed time since any soldier checked their blindside) as function of time. The learned red squad state seems to indicale high security risk. These squad states can modulate voldier-leve heading dynamics. Bottom right: Inferred entity-level states 20:7 (colors) for Soldier 6, superimposed observed time series of heading direction from that soldier's helmet IMU. The light blue state's autoregressive emission dynamics produce a rapid turn to the north. Twice this state persisted long enough for the soldier to reduce security risk in the north (around 19:03 and 19:20). Bottom left: The learned probability that Soldier 6 turns to the north from various sppeier- states (z, rows) depends upon the squad-level states (s, columns). The soldier is most likely to persist in turning north when the squad has security vulnerability (s is red).
6 Discussion & Conclusion
We have introduced cost-cfficient family of models for capturing the dynamics of individual entities evolving in coordinated fashion within shared environment over the same lime period. These models admit elficient structured variational inference in which coordinate ascent can alternate between E-step dynamic programming routinos similar to classic forward-backwant recursions to infer hidden state posteriors at both system- and entity-levels and M-step updates to transition and emission parameters that also use closed-form updates when possible. Across several datasets, we have shown our approach represents a natural way to capture both top-down system-to-entity and bottom-up entity-to-system coordination while keeping costs linear in the number of entities.
Limitations. Our method intentionally prioritizes discrete latent variables for interpretability and efficiency. For some tasks, this choice may be less flexible than more rich continuous representations. We further assume interactions between entities are moderated by a top-level discrete variable. This choice allows runtime to be linear in the number of entities J, but may be less expressive than direct. pairwise interaction. Several coordinate ascent steps in any per-entity FAR-HMM with Gaussian emissions scale quadratically in D due to our choice to parameterize via full covariance matrix. Scaling beyond a few dozen features may require diagonal low-rank parameterizations. Furthermore, the parametric forms of both transitions and emissions in our model allow tractability but clearly limit expressivity compared to deep probabilistic models that integrate non-linear neural nets (Krishnan et al., 2017). Scaling to many more entities would require extensions of our structured VT to process minibatches of entities (Hoffman et al., 2013; Hughes et al., 2015). Scaling to much longer sequences might require processing randomly sampled windows (Foti et al., 2014).
When to favor HSRDM over alternatives? Readers may wonder when to expect our HSRDM has advantages more flexible non-linear models. We suggest that our approach may be favorable when there is clear
15

---

Published in Transactions on Machine Learning Research (08/2025)
top-down coordination among 5-200 entities that could be expressed via simple discrete segmentations by domain expert, and there is interest in both the forecasting and explanatory purposes of the model. When bottom-up coordination alone dominates, interaction graph approaches or neural net methods like GroupNet or AgentFormer may be better. Our basketball case study highlights this: GroupNet slightly outperforms our model in forecasting accuracy. We think this is in part because basketball play is driven by dynamically responding to in-moment conditions (pairwise interactions) rather than rigidly following the prescribed plays called by a coach. Our parameter-effcient approach also has advantages when number of available time series is relatively limited (say few hundred recordings less). When abundant data exists, flexible function approximation methods for forecasting may be preferred.
Future directions. For some applied tasks, it may be promising to extend our two-level system-enttty hierarchy to even more levels (e.g. to represent nested structures platoons, squads, and individual soldiers all pursing the same mission). Additionally, we could extend from rARHMMs to switching linear dynamical systems by adding an additional latent continuous variable sequence between discretes 2 and observations æ in the graphical model, or avoid first-order Markov dependency in generaling state sequences via more flexible distributions parameterized by recurrent neural nets. Lastly, we can add selective and adaptive recurrence functions for modeling entity interactions where current observations influence the type of signal that would be useful feedback for the hidden states.
Broader Impacts
Our work attempts to provide a fundamental framework for modeling interacting entities via explicit mechanisms for top-down and bottom-up influences on group dynamics. As with many technological innovations, for specific applications there may be positive and negative downstream consequences. It is important for researchers to take an active role to avoid misuse and harm.
Many possible applications involve human subjects. We strongly recommend all researchers follow appropriate local regulations and best practices for human subjects research. Each application must carefully consider ethical principles such as respect for persons and beneficence (maximizing benefits while minimizing risks). When using this model for human subjects research, please use care when making specific decisions about what covariates to include. Consider de-identification whenever possible to preserve the privacy of individuals.
Some potential applications of our work involve modeling human teams working in military applications. Our present work focuses exclusively on deidentified data from simulated training exercises via research plan that was approved by a U.S. Army institutional review board. We strongly recommend similar ethics oversight and care in future use cases.
Acknowledgments
This research was sponsored in part by the U.S. Army DEVCOM Soldier Center, and was accomplished under Cooperative Agreement Number W911QY-19-2-0003. The views and conclusions contained in this document are those of the authors and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Army DEVCOM Soldier Center, or the U.S. Government. The U. S. Government is authorized reproduce and distribute reprints for Government purposes notwithstanding any copyright. notation hereon.
Authors KG, EM, and MCH also acknowledge funding from the U. S. National Science Foundation under Growing Convergence Research grant #2428640. MCH is further supported by NSF CAREER award #2338962 and a gift from Apple, Inc.
References
Ahmed M Alaa and Mihacla van der Schaar. Attentive state-space modeling of discase progression. Advances in neural information processing systems, 32, 2019.
16

---

Published in Transactions on Machine Learning Research (08/2025)
Xavier Alameda-Pineda, Vincent Drouard, and Radu Patrice Horaud. Variational inference and learning of piecewise linear dynamical systems. TEEF Transactions on Newral Networks and Learning Systems, 33(8): 3753 3764, 2021.
Michael A Alcorn and Anh Nguyen. baller2vec+++: A look-ahead multi-entity transformer for modeling coordinated agents. arXiv preprint arXiv:2104.11990 2021. URL http:///arxiv.orggbbbb// 111p//0000...110...1111 11980.
Rachel MacKay Altman. Mixed hidden markov models: an extension of the hidden markov model to the longitudinal data setting. Journal of the American Statistical Association, 102(477):201 210, 2007.
Abdul Fatir Ansari, Konstantinos Benidis, Richard Kurle, Ali Canor Türkmen, Harold Soh, Alexander J Smola, Yuyang Wang, and Tim Januschowski. Deep Explicit Duration Switching Models for Time Series. In Advances in Neural Information Processing Systems (NeurIPS), 2021. URL https: http://aaxvvvvoo org/obgg//o abs/ 2110. 13878.
Mehrtash Babadi, Jack M Fu, Samuel K Lee, Andrey N Smirnov, Laura D Gauthier, Mark Walker, David Benjamin, Xuefang Zhao, Konrad J Karczewski, Isaac Wong, al. Gatk-genv enables the discovery of rare copy number variants from exome sequencing data. Nature genetics, 55(9):1589 1597, 2023.
Arindam Banerjoe, Inderjit S Dhillon, Joydcep Ghosh, Suvrit Sra, and Greg Ridgoway. Clustcring on the unit hypersphere using von mises-fisher distributions. Journal of Machine Learning Research, 6(9), 2005.
David Barber, A Taylan Cemgil, and Silvia Chiappa. Bayesian time series models. Cambridge University Press, 2011.
Matthew James Beal. Variational algorithms for approzimate Bayesion inference. University of London, University College London (United Kingdom), 2003.
Yoav Benjamini and Yosef Hochberg. Controlling the false discovery rate: practical and powerful approach to multiple testing. Journal of the Royal statistical socicty: series B (Methodological), 57(1):280-300 -000, 1995.
David M Blei, Alp Kucukelbir, and Jon D McAuliffe, Variational inference: review for statisticians. Journal of the American statistical Association, 112((18):8599- 8777 2017.
James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. JAX: Composable transformations of Python+NumPy programs, 2018. URL http://github com/ gomm//oool/jaaa
Manuel Brenner, Elias Weber, Gcorgia Koppe, and Daniel Durstewitz. Learning interpretable hicrarchical dynamical systems models from time series data. In International Conference on Learning Representations (ICLR), 2025. URL https://axxvvvoor//aaa
04814.
Hung H Bui, Dinh Q Phung, and Svetha Venkatesh. Hierarchical Hidden Markov Models with General State Hierarchy. In Proceedings of the National Conference on Artificial Intelligence (AAAI), 2004.
Zhe Dong, Bryan A Seybold, Kevin P Murphy, and Hung H Bui. Collapsed amortized variational inference for switching nonlinear dynamical systems. In International Conference on Machine Learning, 2020.
Toni V Farle-Randell, Josoph B Wiggins, Julianna Martinez Ruiz, Mehmet Celepkolu, Kristy Elizaboth Boyer, Collin F Lynch, Maya Israel, and Eric Wiebe. Confusion, conflict, consensus: Modeling dialogue processes during collaborative learning with hidden markov models. In International Conference Artificial Intelligence in Education, 2023.
Amirreza Farnoosh, Bahar Azari, and Sarah Ostadabbas. Deep Switching Anto-Regressive Factoriza- tion:Application to Time Series Forecasting. In AAAI Conference on Artificial Intelligence, 2021. URL http://arxiv. org/ abs////s./00.....////// 05135.
Panna Felson, Patrick Luccy, and Sujoy Ganguly. Where will they go? predicting fino-grained adversarial multi-agent motion using conditional variational autoencoders. In Proceedings of the European Conference on Computer Vision (ECCV), 2018.
17

---

Published in Transactions on Machine Learning Research (08/2025)
Shai Fine, Yoram Singer, and Naftali Tishby. The Hicrarchical Hidden Markov Model: Analysis and Applications. Machine Learning, 32, 1998. URL https://link. hatpss////........iegerrrr cww/coooc/aoo/oo http/ootent//d///// 1023/A: 1007469218079, .df..
I. Fisher and A. J. Lec. Time series analysis of circular data. Journal of the Royal Statistical Society: Series B (Melhodological), 56(2):327 339, 1994.
Nick Foti, Jason Xu, Dillon Laird, and Emily Fox. Stochastic variational inference for hidden Markov models. Tn Advances in Neural Information Processing Systems, 2014. URL https://papers: .ips.. cc/paper_ mills/papers /tt///////66566e4444444444511110
Emily B. Fox, Erik B. Sudderth, Michael I. Jordan, and Alan S. Willsky. A sticky HDP-HMM with application to speaker diarization. Annals of Applied Statistics, 5(2A):1020 1056, 2011.
Emily B. Fox, Michael C. Hughes, Erik B. Sudderth, and Michael I. Jordan. Joint modeling of multiple time series via the beta process with application to motion capture segmentation. Annals of Applied Statistics, 1308 4747 4747.pdf.
(3):1281-1313, 2014. URL
Zoubin Ghahramani and GeofTrey E Hinton. Variational learning for switching state-space models. Newral computation, 12(4):831-865) 864, 2000.
Albert Gu, Karan Goel, and Christopher Ré. Efficiently Modeling Long Sequences with Structured State Spaces. In International Conference on Learning Representations (ICLR). arXiv, 2022. URL http: wwt.:/rrrbbb/.......//
James D Hamilton. Regime switching models. Mucrocconomelricc and lime series analysis, pp. 202 209,
2010.
James Douglas Hamilton. Time series analysis. Princeton university press, 1994.
Katherine A. Heller, Yee W. Teh, and Dilan Görür. Infinite hierarchical hidden Markov models. In Artificial Intelligence and Statistics, 2009. URL https://proceeddngs .mlr
http://ccatgga.....rere//////aaa prooss//Ehhhlleer00//aal1rrr0000 pdf.
Matt Hoffman, David Blei, Chong Wang, and John Paisley. Stochastic Variational Inference. Journal of Machine Learning Research, 14(1), 2013.
Michael C Hughes, William T Stephenson, and Frik Sudderth. Scalable adaptation of state complexity for nonparametric hidden markov models. In Advances in Neural Information Pro- cessing Systems (NeurIPS), 2015. URL https://papers. nips. haaaper_iileesppperr//aaa Abstract html.
Thomas Kipf, Ethan Fetaya, Kuan-Chich Wang, Max Welling, and Richard Zemel. Neural Relational Inference for Interacting Systems. In Proceedings of the 35th International Conference on Machine Learning, pp. PMLR, 2018. URL https://proceeednngg mlr. http:///vv///iif/88////
Rahul G. Krishnan, Uri Shalit, and David Sontag. Structured Inference Networks for Nonlinear State Space Models. In AAAI Conference on Artificial Intelligence, 2017. URL httt:://arxiv.orr/abbbba
09869.
Scott Linderman, Matthew Johnson, Andrew Miller, Ryan Adams, David Blei, and Liam Paninski. Bayesian learning and inference in recurrent switching linear dynamical systems. In Artificial Intelligence and Slatislics, pp. 914-99999 922. PMLR, 2017.
Scott Linderman, Annika Nichols, David Blei, Manuel Zimmer, and Liam Paninski. Hierarchical recurrent. state space models reveal discrete and continuous dynamics of neural activity in c. clegans. BioRriv, pp- 621540, 2019. URL https://wwrt http://wwwbbbbb...//////
biorxiv.
Kostya Linou.
NBA
Player Movements,
2016.
URL https://github. http:::/ggttuub.ommmiooooo cotp:/noomm//n///////////
Movements.
18

---

Published in Transactions on Machine Learning Research (08/2025)
Dacheng Liu, Tao Lu, Xu-Feng Niu, and Hulin Wu. Mixed-effects state-space models for analysis of longitudinal dynamic systems. Biometrics, (7(2):476-485, 2011.
Yongtuo Liu, Sara Magliacane, Miltiadis Kofinas, and Efstratios Gavves. Graph Switching Dynamical Systems. International Conference on Machine Learning (ICML), 2023. URL http://arxiv, arrg/bs/222 00370.
Sindy Löwe, David Madras, Richard Zemel, and Max Welling. Amortized Causal Discovery: Learning Tnfer Causal Graphs from Time-Series Data. In Conference on Causal Learning and Rensoning (CLeaR), 2022. URL http://arxiv.. http://axxvv..../aaeaee htrg/abs/2006 10833.
Patrick Lucey, Alina Bialkowski, Peter Carr, Stuart Morgan, Iain Matthews, and Yaser Sheikh. Representing and discovering adversarial leam behaviors using player roles. In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR), 2013.
Rodolfo Metulini, Marica Manisera, and Paola Zuccolotto. Modelling the dynamic pattern surface area in basketball and its effects on team performance. Journal of Quantitative Analysis in Sports, 14(3):117-130, 2018. URL https:// http://www..gruutte..eeee degruyter.... comm/aocument//ii/11o0oooo (11551//jjaa2211-00000
Andrew C Miller and Luke Bornn. Possession Sketches: Mapping NBA Strategies. In MIT Sloan Sports Analytics Conference, 2017.
James Munkres. Algorithms for the assignment and transportation problems. Journal of the Society for Industrial and Applied Mathematics, 5(1), 1957. URL https://doi.onn org/ 10. 1137/0105003.
Josue Nassar, Scott W. Linderman, Monica Bugallo, and Memming Park. Tree-Structured Recurrent Switching Linear Dynamical Systems for Multi-Scale Modeling. In International Conference on Learning Representations (ICLR), 2019. URL http://arxiv http:///aaxxi.../aabsaassee arggabb//811....
N.T. Nguyen, D.Q. Phung, S. Venkatesh, and H. Bui. Learning and detecting activities from movement trajectories using the hierarchical hidden Markov model. In IEEE Compuler Sociely Conference Computer Vision and Pattern Recognition (CVPR), 2005.
Tor Ole Odden and Rosemary Russ. Vexing questions that sustain sensemaking. International Journal of Science Education, 41, 2019. doi: 10.10800/09500093.20113 201.00999...11.1119553.
Dominik Raabe, Reinhard Nabben, and Daniel Mommert. Graph representations for the analysis of multi-agent spatiotemporal sports data. Applied Intelligence, 53(4))3783 3803, 2023.
L. R. Rabiner. A Tutorial Hidden Markov Models and Selected Applications in Speech Recognition. Proc. of the IEEE, (7(2):257-286, 1989.
Syama Sundar Rangapuram, Matthias Seeger, Jan Gasthaus, Lorenzo Stella, Yuyang Wang, and Tim Januschowski. Deep State Space Models for Time Series Forecasting. In Advances in Neural Information Processing Systems (NeurIPS), 2018.
Kristen A Severson, Lana Chahine, Luba Smolensky, Kenney Ng, Jianying Hu, and Sounya Ghosh. Personalized input-output hidden markov models for discase progrossion modeling. In Machine Learning for Healthcare Conference. PMLR, 2020.
R. H. Shumway and D. S. Stoffer. An Approach to Time Series Smoothing and Forecasting Using the EM Algorithm. Journal of Time Series Analysis, 3(4):253 264, 1982. URL /ttps://doi.org// 1111/j. 1467-9892. 1982. tb00349..
Robert H Shumway and David Stoffer. Tine series analysis and its applications, volume 3. Springer, 2000.
Marios Skounkiis, Mark Craven, and Soumya Ray. Hierarchical Hidden Markov Models for Information Extraction. In International Joint Conferences on Artificial Intelligence (IJCAT), 2003.
19

---

Published in Transactions on Machine Learning Research (08/2025)
Ioan Stanculescu, Christopher K Williams, and Yvonne Freer. A Hierarchical Switching Linear Dynamical System Applied to the Detection of Sepsis in Neonatal Condition Monitoring. In Uncerlainly in Artificial Intelligence, 2014. URL https://homepagee inf ed. ac b.c.u//cki
ed.ac.
Jennifer J. Sun, Tomomi Karigo, Dipam Chakraborty, Sharada P. Mohanty, Benjamin Wild, Quan Sun, Chen Chen, David J. Anderson, Pictro Perona, Yisong Yue, and Ann Kennedy. The Multi-Agent Bchavior Dataset: Mouse Dyadic Social Interactions, 2021. URL http://arxiv http:://axxvv.orr/bb///// arg/abs/2104. 02710.
Zachary Terner and Alexander Franks. Modeling Player and Team Performance in Basketball. An- nual Review of Slalistics and Its Application, 8(1):1 2021. URL https:///bii.org/10.
Dick van Dijk, Timo Teräsvirta, and Philip Hans Franses. Smooth Transition Autoregressive Models a Survey of Recent Developments. Econometric Reviews, 21(1):1-47, 2002. URL https://sww. tandfonline. com/doi/ale http//abb////.001 10. 1081 /EEEEEITT-112022
Aoran Wang and Jun Pang. Structural inference with dynamics encoding and partial correlation coefficients.
In International Conference on Learning Representations (ICLR), 2024.
Ezra Webb, Ben Day, Helena Andres-Terre, and Pietro Lió. Factorised Neural Relational Inference for Systems. arXiv, 2019. URL http://arxiv.. htt/a/b//1199 08721.
Tailin Wu, Thomas Breuel, Michael Skuhersky, and Jan Kautz. Discovering Nonlinear Relations with
Minimum Predictive Information Regularization. arXiv, 2020. URL http://arxiv /tt////aaaaarrrrr//////1//1//111 0r98085520011.0885 01885.
Chenxin Xu, Maosen Li, Zhenyang Ni, Ya Zhang, and Siheng Chen. GroupNet: Multiscale hypergraph neural networks for trajectory prediction with relational reasoning. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pottern Recognition, 2022.
Xiuqin Xu, Hanqiu Peng, and Ying Chen. Deep Switching State Space Model (DS*M) for Nonlinear Time Series Forecasting with Regime Switching. International Journal of Forecasting, 2025. URL
https://wwer sciencedirect. http::///ww..siineddreee... http////ccccccccc...caaa/aa////a
Raymond A Yeh, Alexander G Schwing, Jonathan Huang, and Kevin Murphy. Diverse generation for multi-agent sports games. In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pp. 4610-4619, 2019.
Ye Yuan, Xinshuo Weng, Yanglan Ou, and Kris M. Kitani. AgentFormer: Agent-Aware Transformers for Socio-Temporaa Multi-Agent Forecasting. In Proceedings of the IEEE/CVE Internalional Confer- ence on Computer Vision, pp. 9813-908 118-9823, 2021. URL https: /ttee://opeacceees thecvf wwt.cooccoooooo content/ ICCV2021/ttm//uaa AgentFormer_Ageeeee -er_AAgnettaaa Aware transformers forr Socio- Terccor-eemprall... Mauoorr-l Multi-Agent_ Forecasting_Iccy _021_002e1.______1_1111111111111 paper. paper.html.
Eric Zhan, Stephan Zheng, Yisong Yue, Long Sha, and Patrick Lucey. Generating Multi-Agent Trajectories using Programmatic Weak Supervision. In International Conference on Learning Representations, 2019.
URL
wtwp:eeeeenetteeeeeeeeeeeeee bet//forumiddrroo-bbae
20

---

Published in Transactions on Machine Learning Research (08/2025)
Contents
A ARHMM Method Details
21
B HRSDM Method Details
25
C Forccasting Methodology Details
30
D FigureEight: Experiment Details, Settings, and Results
31
E Basketball: Experiment Details, Settings, and Results
33
F Marching Band: Experiment Details, Settings, and Results
38
G Soldiers:
Experiment Details, Settings, and Results
41
Code
Source code for running our proposed HSRDM and reproducing experiments in this paper can be found https://github. com/tuffs-mmm//ean www.://eaaaaaannss-ttttaae- series/
A ARHMM Method Details
As detailed below in Sec. A.1, recurrent autoregressive Hidden Markov Model (HAR-HMM) generalizes a standard Hidden Markov Model by adding autoregressive and recurrent edges to the probabilistic graphical model. Although AAR-HMM models have been previously proposed in the literature (Linderman et al., 2017), we do know of any explicit proposition (or justification) describing how to perform posterior state inference for these models. The literature provides such proposition for (non-recurrent) autoregressive Hidden Markov Model (ARHMMs; e.g., (Hamilton, 1994)), but not for rARHMMs. Hence, provide the missing propositions with proofs here; see Props. A.2.1 and A.2.2. We believe that these explicit propositions can be useful when composing recurrence into more complicated constructions. Indeed, we use them throughout the supplement in order to derive inference for our HSRDMs; for example, see the VES step in Sec. B.2 or the VEZ step in Sec. B.3. In fact, we also utilize the proofs of these propositions when describing how to perform inference with HSRDMs when the dataset partitioned into multiple examples; see Soc. B.7.
A.1 Model
The complete data likelihood for (K, m, no-order recurrent AR-HMM (HAR-HMM) is given by Radon-Nikodým density
p(aLT,211 ()) p(a1 0)p(z1 21,
II p(ee 24- ((-))(()))((111))) 0) 1=2
(1-
(A.1)
Initialisation
transitions
emissions
where æ1.7 are the observations, #1.T E (1,...,KK K} are the discrete latent states, and 0 are the parameters. The FAR-HMM generalizes the standard HMM (Rabiner, 1989), which contains neither autoregressive emissions (blue) nor recurrent feedback (red) from emissions to states. The (K, m, n)-order rAR-HMM gives (K, n ooddr autoregressive HMM (ARHMM) in the special case where
p(a R1-1, ma 0) = p(z. 1._1,
(A.2)
See Fig. A.1 for a probabilistic graphical model representation the special case of first-order recurrence (m 1) and autoregression (n = 1).
21

---

Published in Transactions on Machine Learning Research (08/2025)
Discrete states
Observations
F2
IF
aF's
(a) Recurrent Autoregressive Hidden
Markov Model
(b) Autoregressive Hidden Markov Model (ARHMM).
(TAR-HMM).
Figure A.1: Probabilistic graphical model representation of a Recurrent Autoregressive HMM (HAR-HMM), and its special case, Autoregressive HMM (ARHMM) For simplicity the illustration recurrence, but higher-order dependencies can also be accomodated (see Equation A.1 and Props. A.2.1 and A.2.2). Autoregressive edges shown in blue and recurrent edges are shown red.
first-order autoregression and
Remark A.1.I. (On generalizing HMM with autoregressive emissions and recurvent state transitions.) Let us highlight how a rAR-HMM model generalizes a conventional HMM:
Recurrence: A lookback window of previous observations Tt-mut1 can influence the transitions structure for the current state zt-
pa(::) (=(-1 (=((=))(((1)))
recurrent
Autoregression: A lookback window m previous observations zt-nd-1 can influence the emissions structure for the current observation Zt.
pa(æs) = [ze] U (d(t-n)((t-1))
Autoregressive
Note in particular that each node (observation or state =)) can have many parents among previous observation variables Eitt, but only one parent among state variables 2ıa (namely, the closest in time from the present or past). This assumption will be important when deriving the smoother in Sec. A.2.
A.2 State Estimation
Here we discuss state estimation for the FAR-HMM. We begin with some notation.
Notation A.2.1. Given sequence of observations up to time t, define the conditional probability of the state 2g at a target time & E (1,2,...T) T} via the probability vector Es E AK-1 C R The k-th element this vector is given by po(2, = T1.)). That is,
f - pa(=s ma
pa(=s = maa) pa(=s = K m1:))
Using this notation, we can define three common inferential lasks:
Filtering. Infer the current state given observations & = Po x..))
2. Smoothing. Infer past state given observations = wia), where < t.
Prediction. Predict future state given obscrvations, &. I po Tia), where u > t.
Now we can give Props. A.2.1 and A.2.2, which parallel the presentation of the Kalman filter and smoother in the context of state space models (Shumway & Stoffer, 2000; Hamilton, 1994). In particular, we will present the forward algorithm terms of measurement update (which uses the observation Φ, to transform
What if wanted to relax the specification that the emissions could depend on finite number M of previous states p(æ: 2t 1, 4t- 24-,, 0)) This situation be handled by simply redefining the states terms of tuples #1 (11...11.... such that 27 takes KM possible values, each sequence the back window (Hamilton, 2010, pp.8).
22

---

Published in Transactions on Machine Learning Research (08/2025)
( into ( e) and time update (which transforms into Ee+1 t without requiring observation). These propositions show that fillering and smoothing can be done using the same recursions used in a classical HMM (Hamilton, 1994), except that the variable interpretations differ for both the emissions step and transition step. In the statements and proofs below, we continue to use the same color scheme as was used in Equation A.1 and Fig. A.1, whereby blue designates autoregressive edges and red designates recurrence odges the graphical model. These colors highlight differences from classic HMMs, which lack both types of edges.
Proposition A.2.1. (Filtering a Recurrent Autoregressive HMM.) Filtered probabilities Er (2 x1.t) for Recurrent Autoregressive Hidden Markov Model can be obtained by recursively updating some initialization E1j0 by
Measurement updale.
(E: (-1 ① E:)
&
(T (511 00)
Time update.
Ex+1 | A, Er|
where here E, = (cn1....ex)) (po 24 = k x1=-1)) the (K x 1) vector whose k-th element the emissions density, A, represents the (K K) transition matric whose (k.k')tth element is po(2e41 k' k k,0, represents (K 1) vector of 1s, and the symbol denotes element- beynent-y--elemed multiplication (Hadamard product).
Proof.
Measurement update.
1 po( (ra))
Notation
po(2e, are (01:t-1)
Conditional density Chain rule of probability
po(z (iit-11
pa(ar Re, #1:-))
E 14-1
[ O ct)
E. (T @ Et)
Normalize
Time update.
E++1 + Po(z441 3ia)
Def.
20 k
Law of Total Prob.
po(z 2z a.a)) po(za wxx)
Chain rule of probability
A
E 8th element 1
Notation
Rss
Ath row Ac
As Ee
Def. matrix multiplication
Remark A.2.1. (Inätializing the filtering algorithm in Prop. A.2.1.) Inspired by Hamilton (1994, pp.693), we provide some suggestions for initializing the filtering algorithm of Prop A.2.1. In particular, we can set fo to
23

---

Published in Transactions on Machine Learning Research (08/2025)
Any reasonable probability vector, such the uniform distribution K'1
The maximum likelihood estimate.
The steady state transition probabilities, if they exist.
Proposition A.2.2. (Smoothing a Recurrent Autoregressive Hidden Markov Model.) Smoothed probabilities & po(zt (1:)) for a Hidden Markov Model can be obtained by the recursion
T = { { [5++11 (+) 6+1111
where the formula is initialized by Er (obtained from the filtering algorithm of Prop. A.2.1) and is then iterated backwards for T 1 2 ..... 1, in step analogous the backward pass of the classic forward- backward recursions for plain HMMs (Rabiner, 1989). Here, At represents the (K K) tronsition matrix whose (k,, k')-th element is po(=-11 k/ 2t ki a, the symbol denotes element-une wise mulliplicntion, and the symbol (+) denotes element-wise division.
Proof.
We proceed in steps:
Pa(=: 2e+1, 21:)) (2(2) 2411 2 112) That is, the current state depends future observations 2++1:7 only through the next state 21+1-
Step We show
Step la
We show
po(e #4+1, (1.4+1 po 2441, x1a))
(a #++1, po(2i) #e+1:
split off term from sequence
po(as, #e+1
Ret1,
conditional density
po(e+11 F1+1)
l se, #1:)) (441,
chain rule
Pole
p(ze (141
FPOBN
In the last line, the two canceled terms are equal by FPOBN (the Fundamental Property of Bayes Networks).
Step 1b We show pe(zz 2++1, 21.-2) = pa(z 12+1, 11.+11 By the same argument as in step la (splitting up the sequence, conditional density, chain rule), but replacing
#14+1 + #1442
#++1 2042
the proposition holds if
p(#++2 2t, 2411, x14+1) = p(#242 2441,01111)
that is if we get the same cancelation. And
p(a122 2i, 2411 11441)
Σ
2(2 11
LTP
k=
p( R642 k, 36, #1th, p(=e42 = k x 2441,11441
chain rule, FPOBN
p(#
#142 l 21-1, 2
undo chain rule
km
p(242 #1:11)
undo LTP
"Our proof inspired by the proof given by Hamilton (1994, (007000702) for the ARHMM (i.e, the special of ZAR-HNM in which there recurrent edges).
3The Fundamental Property Bayes Networks node independent of its non- den-descendonts given parents. particular, since 2 is non-descendent #c+1s it independent. of set1 given parents 2++1 and Ritt-
24

---

Published in Transactions on Machine Learning Research (08/2025)
Conclusion
The claim follows from Steps 1a and 1b by induction argument.
p 2t V1:t)
Step 2.
We show that p(z, 24+1
(1:T)
p
F1:T)
p(2++1
flltered
We have
PZ
11:t
smoothed pairwise
smoothed
transition
predicted
2441 #1.T)
p(z44
(1:T) p(a.
1ep #1:T)
chain rule Step
= *(#41 (1er) p(ae
Pe41, m1:e))
p(ze mrra)) p(ze+1 (e.
#1:T)
Bayes rule (on 2nd term)*
p(=++1 m1:a))
Step 3.
We prove the proposition.
K
p(z
FITT)
p(t, 26+11 = lo x1.))
Law Total Prob.
K
Ee
& I
ce+1
Step 2, Notation
1+11
Eet1
|rl Σ
A1
Pull out constant
S44114
E { Af [544 i (+) (+) E44114
Def. matrix multiplication
Remark A.2.2. As in Step the derivation of the smoother in Prop. A.2.2 relies the fact that while each node (observation state) can have many observation parents, it can have only one state parent (namely, the closest in time from the present or past).
A
Remark A.2.3. The filtering (Prop A.2.1) and smoothing (Prop A.2.2) formulae reveal that state estimation for TAR-HMM can be handled for
any order of recurrence and/or autoregression'
any functional form of emissions and transitions
Furthermore, although it was not explicitly represented here, the same formulae hold when there are
Modulation of transitions and emissions by erogenous covariates."
A
B HRSDM Method Details
We now review modeling and inference details for our proposed HSRDM, in the following sections
"To justify the application of Bayes rule, imagine that plays the role of the parameter and #t+1 plays the role of the observed data. The term æi is just conditioning set throughout.
II fact, the proof reveals that the order increase with timestop opening the door to constructions involving exponential weighted moving averages.
6 sequence of vectors {uz} considered be sequence of exogenous covariates if each u, contains information about ze that is not contained in #i:-1 (Hamilton, 1994, pp.692).
25

---

Published in Transactions on Machine Learning Research (08/2025)
B.1 Priors model parameters
The symbol 0 denotes all model parameters for our HSRDM. Using the structure of our model in Equation 2.1, we can expand 0 into constituent components: 0 = (0ss. Pes, Oe.. Qmic)) where tee are the parameters that govern the system-level discrcte state transitions, fes govern the entity-level discrete state transitions, 8e govern the entity-level emissions, and fimit govern the initial distribution for states and regimes. We define a prior over 0 whose iactorization structure reflects this decomposition:
p(0) = p(0za) p(Oea) p(Oce) p(0mms))
(B.1)
As we see in Sec. B.4, this choice of prior simplifies the M-step.
Prior on system-level state transition parameters aag. For the system-level transition probability matrix L x matrix whose entries are all non-negative and rows sum to one, we assume sticky Dirichlet prior (Fox et al., 2011) encourage self-tran transitooons so that in typical samples, one system state would persist for long segments. Concretely, for each row we set
II..... 1 ~ Dir(a..... .......,a K, a......a)
(B.2)
where all I entries have symmetric base value 1.0, and the added value that impacts the self-transition entry (the (j.j)-th entry of the matrix) is set to 10.0. We then set the log transition probability II to the wlement-wise log of II.
Prior on entity-level state transition parameters tes. In our experiments, we used a non-informative prior, p(0ea) x 1. The use of a sticky Dirichlet prior, as was used with the system-level transition parameters, could be expected to produce smoother entity-level state segmentations. Currently, the entity-level segmenta- tions are choppier than those at the system-leve (e.g., compare the bottom-left and top-right subplots of Fig. 5).
Prior on emissions P.c. In our experiments, we used a non-informative prior, p(0.a) 1.
Prior on initial states and observations Oinit. For initial states at both system and entity level, we use a symmetric Dirichlet with large concentration so that all states have reasonable probability a-priori. This avoids the pathology MI. estimation that locks into only one state as a possible initial state early inference due to poor initialization.
B.2 Updating the posterior over system-level states
In this section, we discuss the update the posterior over system-level states; that is, the variational E-S step of Equation 3.2. We find
(xxp
log (50) inis
log Mog
1.-1. (
(1.))
0)
awm log logr(=()) G 6 1 1 81,0)
9 "(
transitions
amissinns
T II
(sa)
plsc st 11 (L.D.
II exp
log Tay (2 k))g0 k)}
Init state
km
transitions
inicial emissions
T
K
exp
log logp(z)) = 1 2(12) = k 21) a(st (3) (1) -K
k)
(B.3)
romaining emilsions
This can be considered as the posterior of an input-output Hidden Markov Model with J independent autore- gressive categorical emissions. The evaluation the transition function is O LDJM.)), where M, is the dimension of the system-level covariates, and where we have assumed that the evaluation of the system-level recurrence function takes DJM, operations, as it would if ge (x( (1:J) v)-1) (z(")... -111, vz-1)", The evaluation of the emissions function is J + KDM))), where M, is the dimension of the entity-level covariates, and where we have assumed that the evaluation of the entity-level recurrence function f takes
DM, operations, it would if fa (1x22) (14 [ (ar u JT Thus, by Props. A.2.1 and A.2.2, filtering
26

---

Published in Transactions on Machine Learning Research (08/2025)
and smoothing can be computed with + (DM))) runtime complexity, under mild assumptions on the recurrence functions. As result, so can the computation of the unary and adjacent pairwise marginals necessary for the VES and M steps.
B.3 Updating the posterior over entity-level states
this section, we discuss the update to the posterior over entity-level states; that is, the variational E-Z step of Equation 3.2. We obtain
(1:1) g( a(>0r))
III
II
exp
log
( 12-21 t,0) = =)
initial state g
2m
transitions (E 2 1)xXX)
pe 2 (2
9) (E 2
(B.4)
initial emission
remaining emissinna (E (T x)
where we have defined 7 A + to denote all timesteps after accounting for the zero-indexing. This variational factor can be considered as posterior of J conditionally independent Hidden Markov Models with autoregressive categorical emissions which recurrently feedback into the transitions. As per Sec. B.2, the transition function can be evaluated with O(TJL(K² + KDM)) runtime complexity, where Me is the dimension of the entity-level covariates, and where we have assumed that the evaluation of the entity-level u(,,) recurrence function f takes DM, operations, it would if fo ( of the emissions function JTT assuming that the emissions distribution has a density that can be evaluated with O(D") operations at each timestep. Thus, by Props. A.2.1 and A.2.2, filtering and smoothing can be computed with O + KD + K KDM))) runtime complexity, under mild assumptions on the entity-level recurrence function and the emissions distribution. As a result, so can the computation of the unary and adjacent pairwise marginals necessary for the VEZ and M steps.
(z (2,))" The evaluation
B.4 Updating the parameters
The M-step updates the transition parameters and emission parameters our HSRDM given recent estimates of state-level posterior 4(5or)) and entity-level posteriors ((=0.r)
This update requires solving the following optimization problem
max (())
where (()) a
(100 0:)) ol Dale
(1:)), (1.) (1:)) log logp(= oor "0:T
[0)
+ logp(0)
(B.5)
log prico
expected complete data likelibood
Based on the structure of the model in Equation 2.1, we can decompose this into separate optimization problems over the different model pieces (fse: lee. Oee, Oinit) by assuming an appropriately factorized prior, as was done in Equation B.1.
To be concrete, for a HSRDM with transitions given by Equation 2.2 and Gaussian vector autoregressive (Gaussian VAR) emissions
S n ( 2 + 1,7 )
(19))
(B.6)
used in Secs. 5.1 and 5.2 we have
A. = (A, II,, tos ( P3 l.o (1A [[Ajk,bjk byk. 1
Using this grouping of the parameters along with the complete data likelihood specification of Equation 2.1 and the prior assumption in Equation B.1, we can decompose the objective as
((0) (0mi)) + Coo()
Σ clu (1)() ) c6: (00)
1=1
27

---

Published in Transactions on Machine Learning Research (08/2025)
We can then complete the optimization by separately performing M-steps for each of the subcomponents of 8. For example, to optimize the parameters governing the entity-level discrete state transitions 00d for cach entity j = .....,,, J,, we only need to optimize
Los (002)
E ofr (a(Korr)
logp(s? 5
( a 2 t-1 81, 0ma)
+ log ((ap(()) log peior
expocted log cutity discrete state transitious
Σ
Σ
log logn(z 10 = + k. k 1: e, (.000) =K,= =K,,) =k) -k)((= =
+logp(fee)
(B.7)
In particular, we do not require the variational posterior over the full entity-level discrete state sequence g(=0,r)) but merely the pairwise marginals q(2(") = N,, =اا = k), obtainable from the VEZ step in Equation 3.2. Similarly, we do not require the variational posterior over the full system-level discrete state sequence ((3erT), but merely the unary marginals q(s; = 6), obtainable from the VES step in Equation 3.2.
The other components of 0 are optimized similarly. In general, the optimization can be performed by gradient descent (e.g. using JAX for automatic differentiation (Bradbury et al., 2018)); although it can be useful to bring in closed-form solutions for the M substeps in certain special cases. For instance, when Gaussian VAR emissions are used as in Equation B.6, the entity emission parameters Oae {[Ajk. byr. Qj k 1 can be estimated with closed-form updates using the sample weights q(2{)) = k)) available from the VEZ-step.
B.5 Variational lower bound
A lower bound on the marginal log likelihood VLBO I log p( (1:3) (1:))
(1:J)
0), is given by (W-)
VLBO [0, = E, log oor so.r 0]] + H (a( (a(=0rr) sorr)
(B.8)
The energy term Eg logg((#))) o.T *0TT sorr [0] is identical to the relevant term in the objective function for the M-step given in Equation B.5. Based on the structure of the model assumed Equation 2.1, the energy term decomposes into separate pieces for initialization, system transitions, entity transitions, and emissions. For example, see Equation B.7 for the piece relevant entity transitions.
(1:J)
(1:J)
Now we consider computation of the entropy H H[((<00 (1:2) :sorr)] in Equation B.8. Since the variational factors ((Sur) and (a(=()))) given respectively by the VES step in Sec. B.2 and the VEZ step Sec. B.3 both have the form of rARHMMs, we can compute the entropy TH [[(=0011) , so:r) E ((( (i) } [((o))):) using for each individual term the entropy for HMMs provided in Eq. 12 of Hughes et al. (2015).
B.6 Smart initialization
We can construct "smart" (or data-informed) initialization of a HSRDM via the following two-stage procedure:
We fit J bottom-level rARHMMs, one for each of the J entities. In particular, the emissions for each bottom-level rAR-HMN are the emissions of the full HSRDM given in Equation 2.3, and the transitions are the entity-level transitions given in Equation 2.2b.
We fit one top-level ARHMM. Here, the J emissions are the entity-level transitions given in Equation 2.2b. The transitions are the system-level transitions given Equation 2.2a. The observations are taken to be the most-likely entity-level states inferred by the bottom-level rARHMMs.
Below we give details on these initializations. In particular, both the bottom-level and top-level models themselves need initializations. We use the term pre-initialization refer the initializations of those models.
28

---

Published in Transactions on Machine Learning Research (08/2025)
B.6.1 Initialization of bottom-level TARHMMS
Here we J bottom-level rARHMMs, one for each the J entities, independently. In particular, the emissions for each bottom-level FAR-HMM are the emissions of the full HSRDM given in Equation 2.3, and the transitions are the entity-level transitions given in Equation 2.2b.
The J bottom-level AARHMSs themselves need good (dala-informed) initializations. As example, we describe the pre-initialization procedure in the particular case of Gaussian VAR emissions, as given in Equation B.6. In particular, we focus on a strategy for pre-initializing these emission parame- ters {A, (k) b (k) ki since the higher-level parameters in the model can be learned via the two-stage initialization procedure.
Pre-initialization.
In particular, for each = 1;.....- J,
a) We assign the observations (3) F1. to one of K states by applying the K-means algorithm to either the observations themselves or to their velocities (discrete derivatives) #2.7 x I:T-1: depending upon user specification. We use the former choice in the FigureEight data, and the latter choice for basketball data.
b) We then initialize the parameters by running separate vector autoregressions within each of the K clusters. In particular, for each state k I,..,,,, K,
a) We find state-specific observation matrix A A and biases b/tt) by applying a (multi-outcome)) linear regression to predict w from the a whenever ag!") belongs to the k-th cluster.
(
b) We estimate the regime-speciice covariance matrices o from the residuals of the above vector antoregresssion.
We initialize the entity-level transition parameters (V, P,J/ to represent sticky transition probability matrix. This implies that we initialize V = 0 for all j.
After pre-initialization, we estimate the independent FARHMMs by using the expectation maximization algorithm. Posterior state inference (i.e. the E-step) for this procedure is justified in Sec. A.2. Note that the posterior state inference for these bottom-level rARHMMs can be obtained by reusing the VEZ step of Equation B.4 by sctting the mumber of system states to L=1.
Expeetation-Maxiiizttion
B.6.2 Initialization of top-level ARHMM
Here we fit top-level ARHMM. In particular, the emissions for the ARHMM are the entity-level transitions of the HSRDM given in Equation 2.2b, and the transitions of the ARHMM are the system-level transitions given in Equation 2.2a. We can perform posterior state inference for the top-level ARHMM by reusing the VES step of Equation B.3 with inputs being the posterior state beliefs on =
(1:J) from the bottom-level rARHMMs.
B.7 Multiple Examples
In some datasets, we may observe the same J entities over several distinct intervals of synchronous interaction. We call cach soparate interval of contiguons interaction an "example". For example, the raw basketball datasot from Sec. 5.2 is organized the end of one play and the beginning of the next, the players might have changed positions entirely, perhaps even having gone to the locker room and back for halftimc.
collection of separate plays, where each play one separate example. Between
Let E be the number of examples. Each example, indexed by e E (1,2,... F), starts at some reference time To and has T. total timesteps, covering the time sequence t E (Te, Te + I,..., ......,T + We'll model each per-example observation sequence x Te an iid observation from our HSRDM model.
(1:))
To efficiently represent such data, we can stack the observed sequences for cach example on top of one
(1:J) that covers all timesteps across all examples, defining another. This yields total observation sequence T
29

---

Published in Transactions on Machine Learning Research (08/2025)
T=I1+72t...77 This representation doesn't it waste any storage on unobserved intervals between examples, easily accommodates examples of arbitrarily different. lengths, and integrates well with modern vectorized array libraries in our Python implementation. As before in the single example case, our computational representation of (1:J) is as 3-d array with dimensionality JJ)
O:T
For properly handling this compact representation, bookkeeping is needed to track where one example sequence ends and another begins. We thus track the ending indices of each example in this stacked representations E (to, ti,.,,,, tg-1.gggggl........ hg) where -1= = to < t < t2 <... < t < tg = T, and where te = Te + Te is the last valid timestep observed in the e-th example for e = I,..., E.
By inspecting the inference updates gives above, including the filtering and smoothing updates for ARR-HMM (see Props. 1.2.1 and A.2.2), we find that
handle this situation follows:
E-steps (VEZ or VES): Whenever we get to a cross-cxample boundary, we replace the usual transition function with an initial state distribution. More concretely, the transition function for the VES step in Equation B.3 is modified so that any timestep t that represents the start of a new example sequence (that is, satisfies t-1 E E) is replaced with *,, and the transition function for the VEZ step in Equation B.4 at such timesteps is replaced with ズ2,, Similarly, the emissions functions such timesteps are replaced with the initial emissions. This maneuver can be justified by noting that for any timestep designating the onset of new example, the initial state distributions play the role of A, and the initial emissions play the rolc of c in Props. A.2.1 and A.2.2.
M-steps: Due to the model structure, the objective function L for the M-step can be expressed as a sum over timestep-specific quanities; for example, see Equation B.7. Thus, in the case of multiple examples, we simply adjust the set of timesteps over which we sum in the objective functions relative to each M substep. We update the entity emissions parameters Oee by altering the objective to sum over timesteps that aren't at the beginning of an example (so we sum over timesteps t where t t- 11 & E). We update the system state parameters tge and entity state parameters tee by altering the objectives to sum only over timesteps that haven't straddled an example transition boundary. That is, we want to ignore any pair of timesteps ( +11++ where i E, so we again sum only over timesteps where t- 1 o E. Finally, we update the initialization parameters Oinit by altering the objective sum over all timesteps that are at the beginning of example.
C Forecasting Methodology Details
Here we detail how we assess model fit (Sec. C.1) and compute forecasts (Sec. C.2). The primary difference between fitting and forecasting is that only the former has access to observations from evaluated entities over time interval of interest. Hence, a good fit is more easily attained. A good forecast requires predictions of the discrete latent state dynamics without access to future observations, whereas fitting can use the future observations to infer the discrete latent state dynamics. However, model fit. is still useful to investigate; for instance, it can be useful to determine if piecewise linear dynamics (including the choice of K, the number of per-entity states) provide a good model for a given dataset.
C.1 Model fit
To compute the fit of the model to (x .... a the j-th entity's observed time series over some slice of integer-valued timepoints [t,..., + us,, we initialize
p = =
And then forward simulate. In particular, for time T in [t...., + u], we do
K
g( ( = k) (
(C.1)
A=1
30

---

Published in Transactions on Machine Learning Research (08/2025)
where 74,k the conditional expectation of the emissions distribution from Equation 2.3 with Radon-Nikodle density p(x) |x = (1-1122))) For example, with Gaussian vector autoregressive (VAR) emissions, we have
(5)
9 - A,, 1 + bj.k
The resulting sequence {ug ..... Pit gives the variational posterior mean for the j-th entity's observed time scries over timepoints [t.....t + u).
C.2 Partial forecasting
By partial forecasting, we mean predicting {ag ................................ N...I.uuuee tubee the observed time series from some to- be-forecasted entities (with indices J C [I......]] over some forecasting horizon of integer-valued (i) hyee from the contextual entities J & timepoints tt. u], given observations {a( .... Fituljeye (j E [1..... J]: $ J] over that same forecasting horizon, as well as observations from all entities over carlier time slices fao (1) ..... ......(()))) -fbbee(
To instantiate partial forecasting, we must first adjust inference, and then perform forward simulation.
Inference adjustment. The VEZ step (Sec. B.3) is adjusted so that the variational factors on the entity-level states the forccasting horizon {q(2))) 2t=a)) )) computed only for the contextual entities (j E Jd}. Likewise, the VES step (Sec. B.2) is adjusted that the variational factor on the system-level states over the forecasting horizon (s.,...., Sttn) is computed from the (i)
observations {a( "e(")....(() 3tt and estimated entity-level states fg(2(") ..... 2l2,)) only from the contextual entities (j Jo}. As result, the M-step on the system-level parameters Oss automatically exclude information from the ne-be-forecaster entities over the forecasting horizon [t..... f u].
2. Forward simulation. Using the adjusted inference procedure from Step 1, we can use the Viterbi algorithm (or some other procedure) to obtain estimated system-states (s,,..., that do not depend on information from the to-be-forecasted entities J over the forecasting horizon [t,..., ,t u]. We then make forecasts by forward simulating. In particular, for time in [t...., ud), sample
1 ~ p(z()) (1) 12
(C.2)
x ~ p(a (9) (2 # 21,0)))000
(C.3)
for all to-be-forecasted entities j E
Note in particular that the dependence of Equation C.2 upon & allows our predictions about to-be-forecasted entities {j J to depend upon observations from the contextual entities {j E J over the forecasting horizon.
D FigureEight: Experiment Details, Settings, and Results
D.1 Data generating process
Example D.1.1. (PigureEight.) Consider model where we directly observe continuous observations (1:J) and where each a. E IR? lives in the plane (i.e. D = 2). We form "Figure Eights" by having the observed dynamics rotate around an "upper circle" C1 with unit radius and center (() & (0, 1)T and a "lower circle" Cz with unit radius and center (2) & (0, -1)T. Entities tend to persistently rotate around one of these circles; however, when the observation approaches the intersection of the two circles C₁ nC2 = {0, 0)}, recurrent feedback can shift the entity's 's dynamics into new state (the other circle). These shifts occur only when the system-level state has changed; these shifts not predictable from the entity-level time series alone. In
31

---

Published in Transactions on Machine Learning Research (08/2025)
particular, have
St 8t-1
h(st)
(D.la)
system transitions
1-1 Cat-GLMX (1 = =v(n))) ((x
())
(D.1b)
entlty transitions
recurrence
transitions
(N(A
2
(D.1c)
observatton dynamics
Here, the notation used follows that of Equation 2.2. Each line of this true data-generating process is explained the corresponding paragraph below.
System-level state transitions. We take the number of system states to be == 2. We set the system state chain {st]/= through deterministic process h which alternates states every 100 timesteps. We emphasize that in the true data-generating process, there is recurrent feedback from obscrvations to system states 8.
Entity-level state transitions. We set entity-specific baseline transition preferences to be highly sticky, p (1 (1-p) where p is close to 1.0 (conerctely, p = .999). By design, these preferences can be p
P
(1-p) (1((1-))
overridden when an entity travels near the origin. We choose the recurrence transformation f: TR. -> R
be the radial basis function f(x)= K exx(( 20 ( which returns large value when the observation ())
x;-1 is
close to the origin. Similarly, we set the weight vector for these recurrent features to nudge observations near
the origin to the system-preferre state. We set vr(t) E IR so entry (fr(#)) = thigh if the entity-level state k is preferred by the system-level state E, and ajow otherwise, with ahigh > alow- Concretely, We set ahigh = 2 and Alow = 2.
Emissions. To construct the entity-level emission distributions for each state (indexed by k), we choose A = Aj to be a rotation matrix with angle 0 = (-1)" 2x for all entity-level states k, where Tj is the entity-specific periodicity and r (0, 1} determines the rotation direction. We may use rotation matrix A to rotate the observation around center b, by constructing dynamics of the form A(x - b) + b; therefore, to construct circle centers that are specific to entity-level states using Equation D.Ic, we set by = (I - (I-A)) 6 (R) for all entities j and all entity-level states k. We set each of the observation noise covariance matrices q₂at) be diagonal, with diagonal entries equal to 0.0001.
7
(A))) to
We simulate data from the FigureFight model (Example D.1.1), where there J = =3 entitics, each with T 400 observations, where the periodicities for each entity are given by (Ti, 72, Ta) (5, 20,4444)
A
D.2 Models and Hyperparameter Tuning
We fit HSRDM with transitions given in Equation 2.2 and Gaussian vector autoregressive emissions in Equation B.6. We set I, = K = 2. We set the entity-level recurrence to a Gaussian radial basis function and no system-level recurrence g. We use a sticky Dirichlet prior on system-level transition parameters (a = 1 and ς = 10). For 'smart' initialization, the bottom-half of the model is trained for 5 iterations while the tolf is trained for 20 iterations. For the K-Means algorithm, the random state is set seed 120. The model then trained for 10 CAVT iterations with 50 iterations per M-Step. We don't. tune any specific hyperparameters of our HSRDM. The model is run across 5 independent initializations (initialization seeds 124). Five forecasting samples are generated per trial (sample seeds 120-124). The optimal MSF across all initialization trials and forecasting samples is recorded. The average MSE for the trial containing the best sample is taken across all samples to demonstrate diversity in sample generation.
HSRDM.
TAR-HMM. collection of rAR-HMM HMM models can be fit as a special case a HSRDM model where the number of system states is taken to be 1. We train three rAR-HMM models with separate strategies: (" ("Idee.", " a "Pool", and " "Concat."). For 'smart' initialization, the bottom-half of the models are trained for 5 iterations while the top-halves are trained for 20 iterations. For the K-Means algorithm, the random state is
32

---

Published in Transactions on Machine Learning Research (08/2025)
set to seed 120. The models are then trained for 10 CAVI iterations with 50 iterations per M-Step. We don't tune any specific hyperparameters of the rAR-HMM. The models are run across 5 independent initializations (seeds 120-124). Five forecasting samples are generated per trial (sample seeds 120-124). The optimal MSE across all initialization trials and forecasting samples is recorded for each model. The average MSE for the trial containing the best sample is laken across all samples from that trial demonstrate diversity in sample generation.
DSARF. We train the Deep Switching Autoregressive Factorization model (Farnoosh et. al., 2021) with several different parameters. We train with several different choices of lags (l) and spatial factors (K), where the number of discrete states (S = 2) is fixed to match our data-generation. For each strategy ("Indep.", ",, "Pool", and " *Concat."), we conduct a grid search across combinations of and K, specifically for 1 = {1,.... i), Vi E {1,-..,10 10} and K E (1,..., value are selected for further experimentation. DSARF Indep. is optimal with spacial factors and lags L: K = 3, L {1.2}: DSARF Pool is optimal with K L f1.2,3, 3, 4. 5, 6, 7}; and DSARF Concat. is optimal with K 5, (1.2).
(il...,00) 10). The hyperparametens that produce the optimal MSE
Once hyperparameters are selected, each model (and strategy) is trained with 500 epochs and learning rate of 0.01, across 10 independent random seeds (seeds 120-129) representing different initializations and forecasts. Two seeds produce clear outlier results (one in Indep. and one in Concat.) from very poor initializations and were excluded from the average results in Fig. 2, as they would have drastically shifted the mcan (in a direction that demonstrates worse perofrmance for the model). Post training for each model, we use the learned parameters to draw long-tery forecasts.
D.3 Recurrence Ablation
Fig. 2, the rAR-HMM baselines are system state ablations for the HSRDM. Here, we include a recurrence ablation for the Fight task. The best MSE across independent random initializations 0.05 and the Ave. MSE with standard error for the random initialization across forecasting samples 2047(..78). The Mdn. MSE across the independent trial sample averages is 2.41. HSRDM without recurrence still performs worse than the HSRDM with recurrence. See the trajectory below in Fig. D.1.
-000
2
390
Timeskage
0
-900
-300
-2
-1.0
-05
0.0
0.5
Figure D.1: FigureFight No Recurrent Feedback
E Basketball: Experiment Details, Settings, and Results
E.1
Dataset
Raw dataset. We obtain NBA basketball player location data for 636 games within the 2015-2016 NBA season from a publicly available repo (Linou, 2016). Each sample provides the quarter of the game, number of seconds left in quarter, time on shot clock, (x.yzz) location of ball, and the (x,y) locations and IDs for the players on the court. The court is represented as the rectangle (0, 96) [o. 50) in the space of squared feet.
Selection of games. We focus on modeling the dynamics in games involving the Cleveland Cavaliers (CLE), the 2015-2016 NBA champions. In particular, out of 40 available games containing CLE, investigate the 31
33

---

Published in Transactions on Machine Learning Research (08/2025)
games containing one of the four most common starting lineups: 1. K. Irving - L. James K. Love J.
Smith - T. Thompson; 2. K. Irving - L. James K. Love - T. Mozgov J. Smith; 3. L. James
- к. Love - T. Mozgov - J. Smith - M. Williams; 4. M. Dellavedova - L. James K. Love - T.
Mozgov J. Smith. Two games had data errors (lack of tracking or event data), which left a total of G = 29 games for analysis.
Downsampling. The raw data is sampled at 25 Hz. Following Alcorn & Nguyen (2021), we downsample to 5 Hz.
From plays to examples. The raw basketball dataset is represented in terms of separate plays (e.g. shot block, rebound offense, shot made). Following Alcorn & Nguyen (2021), we preprocess the datasct so that these plays are non-overlapping in duration. We also remove plays that do not contain one of CLE's four most common starting lineups. For the purpose of unsupervised time series modeling, we then convert the plays into coarser-grained observational units. Although plays are useful for the classification task pursued by Alcorn & Nguyen (2021), play boundaries needn't correspond to abrupt transitions in player locations. For example, the player coordinates are essentially contimuous throughout shot block rebound offense -> shot made sequence mentioned above. Hence, we concatenate consecutive plays from the raw dataset until there is an abrupt break in player motion and/or a sampling interval longer than the nominal sampling rate. These observational units are called events in the main body of the paper (Sec. 5.2). Functionally, these observational units serve as examples (Sec. B.7). That is, when training models, each example is treated as an i.i.d. sample from the assumed model. For the remainder of the Appendix, we refer to these observational units as examples.
By construction, examples have longer timescale than the plays in the original dataset. Examples typically last between 20 seconds and 3 minutes. For comparison, rebound offense play takes fraction of a second.
At the implementational level, we infer an example boundary whenever at least one condition below is met in a sequence of observations:
The wall clock difference between timesteps is larger than 1.2 times the nominal sampling rate.
2. The player's step size on the court (given by the discrete derivative between two timesteps) is abnormally large with respect to either the court's length or width, where abnormally large is defined as having an absolute z-score larger than 4.0.
Court rotation. The location of team's own basket changes at half time. This can switch can alter the dynamics on the court. We would like to control for the direction of movement towards the offensive and defensive baskets, as well as for player handedness. To control for this, we assume that the focal team (CLE)'s scoring basket is always on the left side of the court. When it is not, we rotate the court. 180 degrees around the center of the basketball court. (Equivalently, we negate both the x and y coordinates with respect to the center of the court.) Since the basketball court has a width of 94 feet and a length of 50 feet, its center is located at (5) when orienting the width horizontally. We prefer this normalization strategy to the random rotations strategy of Alcorn & Nguyen (2021), because the normalization strategy allows to learn different dynamics for offense (movement to the left) and defense (movement to the right).
Index assignments. Each sample from our dataset gives the coordinates on the court of 10 players. Here we describe how we map the players to entity indices. Recall that we only model the plays that consist of starters from focal team, CLE. We assign indices 0-4 to represent CLE starters, and indices 5-9 to represent opponents.
Index assignment for CLE is relatively straightforward. Although we model plays from the G games involving four different starting lineups, we can consistently interpret the indices as 0: Lebron James, 1: Kevin Love, J.R. Smith, 3: Starting Center, 4: Starting Guard. Depending on the game, the starting center was either T. Mazgov or T. Thompson. Similarly, the starting guard was either K. Irving, M. Williams, or M. Dellavedova.
34

---

Published in Transactions on Machine Learning Research (08/2025)
Index assignment for the opponents is more involved. The opponent teams vary from game to game, and even a fixed team substitutes players throughout game. There are numerous mechanisms for assigning indices in the face of such player substitutions (Raabe et al., 2023). Although role-based representations are popular (c.g. see (Felsen et al., 2018) or (Zhan et al.; 2019)) because they capture invariants lost within identity-based representations (Lucey al., 2013), we used a simple heuristic whereby we assign indices 5-9 based on the the player's typical positions. The typical positions can be scraped from Wikipedia. We let the model discover dynamically shifting roles for the players via its hierarchical discrete state representation.
One complication in assigning indices from these position labels is that the provided labels commonly blend together multiple positions (e.g. 'Shooting guard small forward' or 'Center power forward'). Should the second player be labeled as center or a forward? What if there are multiple centers? How do we discriminate between two forwards? To solve such problems, we proceed as follows, operating on a play-by-plar basis
1. Assign players to coarse position groups. We first assign players to coarse position groups (forward, guard, center). We assume that each play has 2 forwards, 1 center, and 2 guards. We use indices 5-6 to represent the forwards, index 7 to represent the center, and indices 8-9 to represent the guards. As noted above, a given player can be multiply classified into a coarse position group; however, a reasonable assignment for a player can be made by considering the position labels for the other players who are the court at the same time. To do this, form B, 5 binary matrix whose rows playcrs the team and whose columns represent the coarse position groups. An entry in the matrix is set to True if the player is classified into that position group. We start with the rarest position group (i.e. the column in B with the smallest column sum) and assign players to that position group, starting with players who have the least classifications (i.e. the players whose rows in B have the smallest row sum). Ties are broken randomly. We continue until we have satisfied the specified assignments (2 forwards, 1 center, and guards). Tf is not possible to make such coarse assignments, we discard the play from the dataset.
2. Order players within coarse position groups. This step is only needed for forwards and guards; there is only ordering of the single center. We use an arbitrary ordering of forward positions:
DORWARD_POSITION DRRDERE
ORDERED
"Small forvard / shooting guard",
"Small forward / point guard",
"Snall forward",
"Small forvard / power forward",
"Power forvard / small forward",
"Power forvard",
"Power forward / center",
"Center power forward",
"Shooting guard small forward",
]
and guard positions by
CUARD_ POSITIONS COREEEEESSDRDDEEDDDDEE =
"Small forward / shooting guard",
"Shoating guard small forward",
"Shooting guard",
"Shooting guard point guard", "Point guard shooting guard", "Point guard" "Combo guard"
]
For each players assigned to position group in {forward, guard), order the players in terms of their location of their position the above lists. Ties are broken randomly.
Normalization
To assist with initialization and learning of parameters, we normalize the player locations the court. from the rectangle [0, 96) x [0,50] in units of feet to the unit square [0, 1] x [0, 1].
35

---

Published in Transactions on Machine Learning Research (08/2025)
E.2 Models and Hyperparameter Tuning
Here we model 10 basketball player trajectories on the court with an HSRDM with Gaussian vector autoregressive emissions; that is, we use
HSRDM.
St $1-1. (1:J) Cat-GLML
Tt' ex
(1:J) A gy ( (0 vz-1)
(E.1)
system transitions
endogenous transition preferences
bias from recurrence and covariates
( (1) -t- 1-1 Cat-GLMM
(P) (P(ot))e (
V (())) (=
(E.2)
entity transitions
endogenous transition proferences
biss from recurrence and covaristes
N (14
(2(
( +
Q (1
x 1-
(F.3)
observation dynamics
where x, E (10,1) x [0, 1I) gives player j's location on the normalized basketball court at timestep t.
Our system-level recurrence gu (2/-1) (1:J) v--1) x reports all player locations at to the system-- level transition function, allowing the probability of latent game states to depend on player loca- tions. Inspired by Linderman et al. (2017), our entity-level recurrence function fo u = Ne < 0.0), > 100, <0.0), [e > 1.0)),, where x(1) the d-th coordinate of of and [[+] is the indicator function, reports an individual player's location ,,,2), (and out-of-bounds indi- cators) to that player's entity-level transition function, allowing each player's probability of remaining in autoregressive regimes to vary in likelihood over the court.
(1:))
We set the number of system and entity states to be L = 5 and K = 10 based on informal experimentation with the training set; we leave formal setting of these values based on the validation set to future work. For the sticky Dirichlet prior on system-level transitions, as given in Equation B.2, we set a = 1.0 and к = 50.0 so that the prior would put most of its probability mass on celf-transition probabilities between .90 and .99.
We initialize the model using the smart initialization strategy of Sec. B.6. We pre-initialize the entity emissions parameters Osn by applying the k-means algorithm to each player's 's discrete derivatives (so long as consecutive timesteps do not span an example boundary). We pre-initialize the entity state parameters Bes by setting P to be the log of sticky symmetric transition probability matrix with a self-transition probability of 0.90, and by drawing the entries of V i.i.d from a standard normal. We pre-initialize the system state parameters Oan by setting II to be the log of sticky symmetric transition probability matrix with self- tef-siransition probability of 0.95, and by drawing the entries of A i.i.d from a standard normal. We pre-initialize the initialization parameters Oinis by taking the initial distribution to be uniform over system states, uniform over entity states for each entity, and standard normal over initial observations for each entity and each entity state. We execute the two-stage initialization process via iterations of expectation-maccinzatioo maximization for the J bottom-half rARHMMs, followed by 20 iterations for the top-half ARHMM.
We our CAVT algorithm for 2 iterations, informal experimentation with the training set suggested this was sufficient for approximate ELBO stabilization.
FRRHMss. By ablating the top-level discrete "game" states (i.e., the system-level switches) in the HSRDM, we obtain independent TARHMMs (Linderman et al., 2017), one for each of the J = 10 players. More specifically, by removing the system transitions in Equation E.1 from the model, the entity transitions simplify 1 2 x st) = p() ()) mT(2), bocause the entity transition parameters simplify as pf++) P,
p(2(1) (1) (1)
and L (+)) J. As a result, the J bottom-level rARHMSS are decoupled. Implementationally, this procodure is equívalent to an HSRDM with L system states. Initialization and training is otherwise performed identically as with HSRDM.
36

---

Published in Transactions on Machine Learning Research (08/2025)
HSDM. By ablating the multi-level recurrence from the HSRDM, we obtain a hierarchical switching dynamical model (HSDM). This can be accomplished by setting Ie = 0 in Equation E.1 and So 0 in Equation F.2.
Initialization and training is otherwise performed identically as with HSRDM.
Agentformer. AgentFormer (Yuan et al., 2021) multi-agent (i.c. multi-entity) variant of a transformer model whose forecasts depend upon both temporal and social (i.e. across-entity) relationships. Unless otherwise noted, we follow Yuan et al. (2021) in determining the training hyperparameters. In particular, our prediction model consists of 2 stacks of identical layers for the encoder and decoder with a dropout rate of 0.1. The dimensions of keys, queries and timestamps for the agentformer are set to 16, while the hidden dimension of the feedforward layer is set to 32. The number of heads for the multi-head agent aware attention is 8 and all MLPs in the model have a hidden dimension of (512,256). The latent code dimension of the CVAE is set to 32, and the agent connectivity threshold is set to 100. Because the basketball training datasets have many more examples than the redestrian trajectory prediction experiments in Yuan et al. (2021) (which only have 8 examples), we train the agentformer model and the DLow trajectory sampler for 20 epochs each (rather than 100) to keep the computational load manageable. We therefore apply the Adam optimizer with learning rate of 10 rather than 10- to accommodate the reduced number of epochs. Also, to match the specifications of the evaluation strategy from Sec. F.3, we set the number of future prediction frames during training to 30, and the number of diverse trajectories sampled by the trajectory sampler to 20. We ensure convergence by tracking the mean-squared error.
SNLDS. At suggestion of an anonymous reviewer, we also compared to a recent non-linear model which does not coordinate interactions across entities, but provides discrete and continuous latent variables as well recurrent feedback: SNLDS (Dong et al., 2020). " SNLDS inherently models only entity-level time series by design; it is not originally designed to model many interacting entities. We tried the Independence entity-to-systen strategy (see Sec. 5.1 of our paper) to adapt it to model many entities. The precise hyperparameter settings can be found by inspecting the basketball. yamm configurations file used to run SNLDS (see our codebase for link). Unless otherwise noted, we follow Ansari et al. (2021)'s modeling of the electricity duration dataset to determine the training hyperparameters (e.g. trainable covariance matrices, transformer within the emissions network with the same network size, etc.). Paralleling the clectricity duration datasct, we set the dimensionality of the latent continuous state (4) to be double that of the observed variables (2). For consistency with the HSRDM, we set the number of entity states to be K = 10. We also allowed melti-level recurrence to both the continuous and discrete latent chains to parallel the multi-level recurrence of the HSRDM. We trained each basketball player's model for 20,000 iterations (using the first 1,000 as warmup iterations for setting the learning rate). We ensure convergence by examining trace plots to verify that the ELBO-based training objective converged after many epochs, and by verifying that reconstructions of past data looked reasonable.
E.3 Evaluation strategy
We divide the G' = 29 total games into 20 games to form candidate training set, 4 games to form validation set (for setting hyperparameters), and 5 games to form a test set. Of the first 20 games within our candidate training set, we construct small (1 game), medium (5 games), and large (20 games) training sets. The small, medium, and large training sets contained 20, 215, and 676 examples, respectively.
The test set contained 158 examples overall. However, we required that each example be at least 10 seconds long (i.e. 50 timesteps) be included in the evaluation run. This exclusion criterion left E = 75 examples. For each such example, we uniformly select a timepoint T E Tnnnn T Ieugth to demarcate where the context window ends. We set T -nin- 4 seconds (i.e. 20 timesteps) and
T lerhhhst-length = 6 seconds (i.e. 30 timesteps). The first [0, seconds are shown to the trained model as context, and forecasts are made within the forecasting window of : [T F=
+ 1,7" + +T Tforecsst-eenghh seconds.
We used the open source implementation of the SNLDS model provided by a successor package for RED-SDS (Ausari et al., 2021): attps://githut com/abdulfatir//EDDDs We chose this because (1) this the newest suggested codcbase, should be the most oo-ate-ff-the-rtt and (2) as the newest codebase, should be easler install and maintain. (Older code from 2017 is often tough to get working on newer hardware, such the non-Intel chips of newer Macbooks.)
37

---

Published in Transactions on Machine Learning Research (08/2025)
For fixed example forecasting sample s, player j, and forecasting method m, we summarize the error in forecasted trajectory by mcan forecasting crror (MFE)
Σ F
Σ d=D
(E.4)
where Re,.j. is the true observation on example e at time for player j on court dimension d, and Le,.jjddm.s is the forecasted observation by forecasting sample s using forecasting method m. So MFBme,,d gives the average distance over the forecasting window between the forecasted trajectory and the true trajectory.
To quantify performance of forecasting methods, we define a model's ezample- wise mean forecasting error
MPEm;e
MFE
(E.5)
SJ
Taking the mean of MFEme and its standard error lets us quantify a model's typical squared forecasting error on an example, well as our uncertainty, with
R
(MFBme
MFEm.)
MFE
MFEm;e:
(MEE..)
(E.6)
Although in Sec. E.1, we described normalization of basketball coordinates to the unit square for the purpose of model initialization and training, when evaluating models, we convert the forecasts and ground truth back to unnormalized coordinates, so that MFE has units of feet. That is, we represent observations Teaddd and forccasts Fe..jjd.mm. on the basketball court (of size [0, 94] x 50] fect). Thus /MSE,.. can be interpreted as a model's typical amount of error in feet on the court at a typical timepoint in the forecasting window (but of course forecasting error tends to be lower timepoints closer T than farther from T*).
F Marching Band: Experiment Details, Settings, and Results
F.1 Data generating process
We introduce MarchingBand, synthetic dalaset consisting of individual marching band players (fentities") moving across a 2D field in a coordinated routine to visually form sequence of letters. Each observation is a position x²(1) R2, with the unit square centered at (0.5, 0.5) representing the field. Each entity's position over time the unit square follows the current letter's prescribed movement pattern perturbed by small-scale iid zero-mean Gaussian noise. Each state is stable for 200 timesteps before transitioning in order to the next state. When reaching a field boundary, typically the player is reflected back in bounds. However, with some small chance, an entity will continue out-of out-of-bounds bounds (OOB, æy & (0, 1)). When enough players become OOB (up to a user-controlled threshold), this bottom-up signal triggers the current system state immediately to a special "come together and reset" state, denoted 'C'. During the next 50 timesteps, all players move to the center, then return to repeat the most recent letter before continuing on to remaining letters. Note that the true data- generrtnng process does not come from HSRDM generative model.
For our experiments, we set J 64, the OOB threshold to be 11, and the letter sequence to spell 'LAUGH'. We use N = 10 independent examples for fitting, where each sequence is of a different length according to the number of "C" states triggered. We observe corresponding time lengths for each sequence: (1000, 1050, 1050, 1000, 1000, 1000, 1050, 1050, 1000, 1100). The total time T across all sequences is 10300, where end-times for each example are provided the model {1000, 2050, 3100, 4100, 5100, 6100, 7150, 8200, 9200, 10300).
F.2 Models and Hyperparameter Tuning
HSRDM (with and without recurrence). For all methods, we fairly provide knowledge of the true number of system states, 6 A,U, H. C). Models with system states (HSRDM and its recurrent ablation) set. 5 6 and have K = 4 entity states. We set system-level recurrence g to count the number
38

---

Published in Transactions on Machine Learning Research (08/2025)
of entities out of bounds and the entity-level recurrence f to the identity function. The emission model is set to a Gaussian Autoregressive. All system and entity states are initially set. to uniform distributions, and we use sticky Dirichlet prior on system-level transition parameters (a = 1 and K = 10). For 'smart' initialization, the bottom-half and top- of the model are only trained for a single iteration, such that little to no initialization used for this experiment. The K-Means algorithm is set to a random state with a seed 120. The model is then trained for 10 CAVT iterations with 50 iterations per M-Step. We don't tune any specific hyperparameters. The model is run across independent initializations (with recurrence: seeds 120-123, 126 and without recurrence: seeds 120-124). Seeds (such 124 for the HSRDM with recurrence) that do not lead to any effective optimization (NAN values) are discarded and an alternative seed is attempted. This procedure is the same for all competitor models, which also tend to produce ineffective training due to randomly poor ionss The best system-state classification accuracy across all trials is recorded for the models.
AAR-HMM. A collection of J HAR-HMM models can be fit as a special case of a HSRDM model where the number of system states is taken to be L = 1. The number of entity states is taken to be K = with a uniform prior over possible states. For 'smart' initialization, the bottom-hall half and top-half of the model are only trained for a single itcration, such that little to no initialization is used for this experiment. The K-Means algorithm is set to a random state with seed 120. The model is then trained for 10 CAVI iterations with 50 iterations per M-Step. We don't tune any specific hyperparameters of the AM. The model is trained across 5 independent initializations (seeds 120-122 and 124-126). The resulting most likely states for each entity at each timepoint are clustered with a K-Means algorithm (random state seed = 120). and are compared with the ground truth system-level accuracy for each initialization. The best classification accuracy is recorded.
DSARF. We train the Deep Switching Autoregressive Factorization model (Farnoosh et al., 2021) with several different parameters. We train with several different choices of lags (1) and spatial factors (K), where the number of discrete states (S=0) = 6) is fixed to match our data data-generation. We use the "Concat")) strategy, and conduct a grid search across combinations of and K, specifically for 1 = -11, [],....), Vi E fl, 10, 200} and K E (I...., 10, 15, 25, 30). These numbers for search were selected based on prior work in Ref. (Farnoosh et al., 2021), such that our values match the hyperparameters used for datasets with a similar size. Also, the selections are based on the knowledge that our dataset changes discrete states around every 200 time- steps. The hyperparameters that produce the optimal system classification accuracy are selected for further experimentation, which we observe to be: I= 1= = (1,..., 200} and K = 25.
Once hyperparameters are selected, each model is trained with 200 epochs (as indicated to be sufficient in the hyperparameter search) and a learning rate of 0.01, across 10 independent random initializations (seeds 120-129), The resulting most likely states for each entity at each timepoint are clustered with K-Means algorithm (random state sced 120), and are compared with the ground truth system-level accuracy for each initialization. The best classification accuracy is recorded.
F.3 Experiments with up to 200 entities
To demonstrate that the HSRDM can train and learn the system-level hidden states in the MarchingBand dataset with a larger number of entities J > 64, we run experiment with J = 200 entities. To simplify our experiment, we left out the cluster state C that gets triggered when certain number of entities go out of bounds. Instead, we only have 5 = 5 = (L, A, U, G. H}. All other variables (number of dimensions (D = 2), priors, initialization procedure, and training procedure) are kept the same as in the J = 64 case. Across 5 initializations, we obtain a best system state classification accuracy of 100%, a median accuracy of 100%, and a average accuracy of 90%. Two out of the five initializations only obtained an accuracy of about 70%, which brought down the total average. Varying results based on unlucky initializations for this type of model is common.
39

---

Published in Transactions on Machine Learning Research (08/2025)
F.4 Experiments on Larger Feature Dimensions
Experiments in the main paper focus on interpretable modeling of data with D 2 fealure dimensions recorded each timestep. To demonstrate that the HSRDM can train and learn the system-leve hidden states in the BardhingBand datasct with larger number of dimensions D > 2, we run experiments with = {10, 30) dimensions. While 2 is the true (z, position of the marching band player, D = { f10, 303 are the true positions with added gaussian noise (0, 0.0004). To simplify this experiment, we left out the cluster state C that gets triggered when a certain number of entities go out of bounds. Instead, we only have S = 5 A,, G, H). All other variables (number of entitics (J 64), priors, initialization procedure, and training procedure) are kept the same as in the J = 64 case.
Across 5 initializations, the average and standard deviation system level classification accuracies for D {2, 10,10,30 10) are shown in Fig. F.1. We observe that the best classification accuracies remain high for all D; however, the average decays as increases (D = 2 88%, = 10: : 80%, D = 30 75%). Since all sizes were rum with the same number of initialization itcrations, number of training examples, and the same number of CAVI iterations, we hypothesize that D 30 might achieve better performance with more careful tuning of these dyper-parametee In Fig. F.1, also show the times for single trial of each D size in minutes for an entire training iteration.
100
70
Runtime in Minutes
60
% Accuracy
90
50
85
40
30
75
20
70
10
10
15
20
25
30
5
10
15
20
25
30
Number of Dimensions
Number of Dimensions
(a) Accuracy vs. Dimension Size
(b) Runtime vs. Dimension Size
Figure F.1: MarchingBand system level accuracy and runtimes for {2, 10, 30).
40

---

Published in Transactions on Machine Learning Research (08/2025)
G Soldiers: Experiment Details, Settings, and Results
For the visual socurity experiment, based on a quick exploratory analysis, we set K = 4 and L. = 3. For the sticky Dirichlet prior on system-level transitions, as given in Equation B.2, we set a = 1.0 and K = 50.0, so that the prior would put most of its probability mass on self-transition probabilities between .90 and .99.
41


Reviews:
[{"review_id": "GJv6a26K92", "reviewer": "Reviewer_iSuN", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths\n\nI think this an absolutely fantastic paper and was written beautifully! The motivation for the approach was expertly told and the method was described concisely. I also found the breadth of experiments to be amazing and highlights the potential of this model.\n\n# Weaknesses\n\nThere are no major weaknesses to this work in my opinion. I think the only thing that I would add is from the writing, it is not clear whether the recurrence parameters, $\\phi, \\psi$ are optimized in the varational M-step.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "In this paper, the authors introduce an extension to recurrent auto-regressive hidden Markov Models  that allows for segmentations of systems involving multiple interacting agents. Specifically, the authors introduce a hierarchical latent state that models the group dynamics of the agent along with agent level latent states. This hierarchy, combined with a recurrent connection from the observations, allows for a simple, yet expressive generative model.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I have no requested changes!", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "xlp57h3LYx", "reviewer": "Reviewer_p3rW", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1. Clear and novel integration of top-down and bottom-up feedback mechanisms, effectively capturing complex group dynamics.\n2. Efficient variational inference algorithm with runtime linear in entity count, enabling scalability.\n3. Strong empirical validation with both synthetic and real datasets.\n4. Practical interpretability of learned states and dynamics, suitable for applications requiring domain knowledge integration .\n\nWeaknesses:\n1. Potential scalability concerns for datasets with higher dimensionality of observations.\n2. Experiments largely focus on moderately sized problems (in terms of entities and complexity); larger-scale scenarios are less explored.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces a novel modeling framework, called Hierarchical Switching Recurrent Dynamical Models (HSRDM), designed to infer coordinated group dynamics from multiple interacting time series. Key contributions include:\n1. This paper proposes a hierarchical latent switching-state structure, which includes a system-level Markov chain that represents global state dynamics influencing individual entities, and entity-level chains governing individual trajectories.\n2. This work also proposes system-level states influence entity dynamics (top-down), while observations from entities influence system-level transitions (bottom-up), improving interpretability and forecasting accuracy.\n3. Experiments showing the model’s effectiveness on synthetic datasets (FigureEight, MarchingBand) and real-world scenarios (Basketball player trajectories, soldier behavior), demonstrating competitive forecasting performance and interpretable insights.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The following changes are recommended, primarily to strengthen the work:\n1. Please include experiments with higher-dimensional observational data or a significantly larger number of entities to demonstrate scalability and robustness. (Strengthens work)\n2. Clearly highlight conditions under which the proposed HSRDM would significantly outperform or underperform relative to neural network baselines such as transformers or RNN-based models.\n(Strengthens work)\n3. Please add a detailed ablation study explicitly evaluating the individual impacts of top-down vs. bottom-up recurrence across all experiments. (Strengthens work)\n4. The description of NRI (Kipf et al., 2018) in the related work is not accurate. It always works with the assumption of a fully connected structure. The connectivity of the graph does not have any influence on the runtime. Actually, it works better with dense graphs.\n4. As some idea look similar to the work in the field of relational inference, it would be better to include them in the related work. E.g.: Node-level variational autoencoder for efficiency: A. Wang and J. Pang. Structural inference with dynamics encoding and partial correlation coefficients. ICLR 2024.", "broader_impact_concerns": "The paper currently does not explicitly address potential ethical concerns or misuse scenarios. A broader impact statement would strengthen the submission, explicitly considering issues such as privacy concerns when modeling individual-level data or unintended biases when deploying group dynamic models for sensitive tasks (e.g., military or surveillance applications).", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "f58wdZPDlX", "reviewer": "Reviewer_Pthk", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- The paper is generally very well-written, and the presented problem seems very interesting.\n- The model is parameter efficient, and approximate inference is justified given the prohibitive cost of exact inference.\n- The experiments show direct applications to real-world scenarios, and advantages of the proposed approach in terms of interpretability.\n\n**Weaknesses**\n- The experiments only consider very low-dimensional setups $D=2$ in all settings.\n- The model might suffer from underfitting in certain complex scenarios, due to following aspects:\n  - Assuming only discrete latent variables.\n  - Allowing only first-order dependencies: $p(x_{t}^j | x_{t-1}^j, z_{t-1}^j)$. Could we replace the proposed distribution with a RNN?\n  - Only allowing entity-level interactions, i.e. $x_{t}^j$ only sees $j$, although it can see other variables through $s_t$ via $z_t^j$, but with limitations as $s$ and $z$ are discrete. This looks very restrictive, specially in complex data such as NBA.\n- **Some important baselines might be missing.** Although the aim of the method is to provide a parameter efficient model that provides data interpretability thanks to the latent space simplicity, comparisons with more powerful baselines (e.g. with continuous latent variables) should be present. See **Requested Changes**.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposes Hierarchical Switching Recurrent Dynamical Model (HSRDM), a state-space model with only discrete latent variables, which is parameter efficient and is aimed to improve interpretability of the learned dynamics. The model assumes independent autoregressive hidden markov models, where an additional discrete latent variable induces changes in the global dynamics. Experiments on both synthetic and real-world data demonstrate the improved interpretability of global and individual system dynamics, and show decent forecasting results despite model simplicity.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Text**\n\n- The \"optional exogenous covariates\" mentioned in Eqs. 2.2ab $u_{t-1}^j$ and $v_{t-1}$ might be a bit confusing. Would it be possible to add one sentence connecting these to the example of including the ball's position in NBA?\n- In section 2, the following statement is imprecise:\n> If we remove the top-level system states ... our HSRMD reduces to Linderman's recurrent autoregressive HMM.\n  \n  I believe the statement is not correct, as Linderman does not consider individual element. You should specify this per-entity modelling in your case to avoid confusions. For example, you can talk about some per-entity rAR-HMM instead.\n\n**Experiments**\n\nFor acceptance, I consider the following point should be addressed. From my understanding, it might be expected to see worse results, but still, it is critical to understand whether HSRDM suffers in terms of forecasting performance despite improved interpretability and parameter efficiently.\n\n- Only 3 baselines are present: Agentformer and rAR-HMM in NBA; and rAR-HMM and DSARF in MarchingBrand. It would be very important to include additional baselines to understand **whether the improved interpretability comes at the cost of forecasting performance.** Below are some baselines that could be used. All of them have available code.\n  - NRI from Kipf et al. (2018) [1]: It would be interesting to check how does your work compare to interective systems, considering this work is mentioned in your related work section.\n  - At least 1 baseline that considers Switching Dynamical Systems should be present. Considering you already have rAR-HMM, it might be interesting to consider rSLDS from Linderman et al. (2017) [2], or one of SNLDS (Dong et al. 2020) [3], RED-SDS (Ansari et al. 2021) [4]. \n  - There is this interesting recent paper, Graph Swicthing Dynamical Systems (GRASS) (Liu et al 2023 [5]), which considers SDSs with individual elements, similar to your case. Here a comparison might be very interesting to highlight the advantages of your method in terms of interpretability.\n\nBelow are some miscellaneous points:\n\n- Would it be possible to try the NBA experiment replacing Eq. (2.3) with a Recurrent Neural Network instead?\n-  All the experiments consider 2D trajectories. Would it be possible to devise a synthetic experiment with increasing $D$ to understand the problems outlined in **Section 6: Limitations**?\n\n**References**\n\n[1] Kipf, Thomas, et al. \"Neural relational inference for interacting systems.\" International conference on machine learning. Pmlr, 2018.\n\n[2] Linderman, Scott W., et al. \"Recurrent switching linear dynamical systems.\" arXiv preprint arXiv:1610.08466 (2016).\n\n[3] Dong, Zhe, et al. \"Collapsed amortized variational inference for switching nonlinear dynamical systems.\" International Conference on Machine Learning. PMLR, 2020.\n\n[4] Ansari, Abdul Fatir, et al. \"Deep explicit duration switching models for time series.\" Advances in Neural Information Processing Systems 34 (2021): 29949-29961.\n\n[5] Liu, Yongtuo, et al. \"Graph switching dynamical systems.\" International Conference on Machine Learning. PMLR, 2023.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)