TREE-STRUCTURED RECURRENT SWITCHING LINEAR DYNAMICAL SYSTEMS FOR MULTI-SCALE MODELING
Year: 2019

Abstract:
Many real-world systems studied are governed by complex, nonlinear dynamics. By modeling these dynamics, we can gain insight into how these systems work, make predictions about how they will behave, and develop strategies for controlling them. While there are many methods for modeling nonlinear dynamical systems, existing techniques face a trade off between offering interpretable descriptions and making accurate predictions. Here, we develop a class of models that aims to achieve both simultaneously, smoothly interpolating between simple descriptions and more complex, yet also more accurate models1. Our probabilistic model achieves this multi-scale property through a hierarchy of locally linear dynamics that jointly approximate global nonlinear dynamics. We call it the tree-structured recurrent switching linear dynamical system. To fit this model, we present a fully-Bayesian sampling procedure using Pólya-Gamma data augmentation to allow for fast and conjugate Gibbs sampling. Through a variety of synthetic and real examples, we show how these models outperform existing methods in both interpretability and predictive capability.

heading: 1 INTRODUCTION
text: Complex systems can often be described at multiple levels of abstraction. A computer program can be characterized by the list of functions it calls, the sequence of statements it executes, or the assembly instructions it sends to the microprocessor. As we zoom in, we gain an increasingly nuanced view of the system and its dynamics. The same is true of many natural systems. For example, brain activity can be described in terms of high-level psychological states or via detailed ion channel activations; different tasks demand different levels of granularity. One of our principal aims as scientists is to identify appropriate levels of abstraction for complex natural phenomena and to discover the dynamics that govern how these systems behave at each level of resolution.
Modern machine learning offers a powerful toolkit to aid in modeling the dynamics of complex systems. Bayesian state space models and inference algorithms enable posterior inference of the latent states of a system and the parameters that govern their dynamics (Särkkä, 2013; Barber et al., 2011; Doucet et al., 2001). In recent years, this toolkit has been expanded to incorporate increasingly flexible components like Gaussian processes (Frigola et al., 2014) and neural networks (Chung et al., 2015; Johnson et al., 2016; Gao et al., 2016; Krishnan et al., 2017) into probabilistic time series models. In
1This work was supported by National Science Foundation (NSF IIS-1734910, CCF-1617986, HRD-1612689) and National Institute of Health (NIH R01EB026946). SWL was supported by the Simons Collaboration on the Global Brain (SCGB-418011).
neuroscience, sequential autoencoders offer highly accurate models of brain activity (Pandarinath et al., 2018). However, while these methods offer state of the art predictive models, their dynamics are specified at only the most granular resolution, leaving the practitioner to tease out higher level structure post hoc.
Here we propose a probabilistic generative model that provides a multi-scale view of the dynamics through a hierarchical architecture. We call it the tree-structured recurrent switching linear dynamical system, or TrSLDS. The model builds on the recurrent SLDS (Linderman et al., 2017) to approximate latent nonlinear dynamics through a hierarchy of locally linear dynamics. Once fit, the TrSLDS can be queried at different levels of the hierarchy to obtain dynamical descriptions at multiple levels of resolution. As we proceed down the tree, we obtain higher fidelity, yet increasingly complex, descriptions. Thus, depth offers a simple knob for trading off interpretability and flexibility. The key contributions are two-fold2: first, we introduce a new form of tree-structured stick breaking for multinomial models that strictly generalizes the sequential stick breaking of the original rSLDS, while still permitting Pólya-gamma data augmentation (Polson et al., 2013) for efficient posterior inference; second, we develop a hierarchical prior that links dynamics parameters across levels of the tree, thereby providing descriptions that vary smoothly with depth. The paper is organized as follows. Section 2 provides background material on switching linear dynamical systems and their recurrent variants. Section 3 presents our tree-structured model and Section 4 derives an efficient fullyBayesian inference algorithm for the latent states and dynamics parameters. Finally, in Section 5 we show how our model yields multi-scale dynamics descriptions for synthetic data from two standard nonlinear dynamical systems—the Lorenz attractor and the FitzHugh-Nagumo model of nonlinear oscillation—as well as for a real dataset of neural responses to visual stimuli in a macaque monkey.

heading: 2 BACKGROUND
text: Let xt ∈ Rdx and yt ∈ Rdy denote the latent state and the observation of the system at time t respectively. The system can be described using a state-space model:
xt = f(xt−1, wt; Θ), wt ∼ Fw (state dynamics) (1) yt = g(xt, vt; Ψ), vt ∼ Fv (observation) (2)
where Θ denotes the dynamics parameters, Ψ denotes the emission (observation) parameters, and wt and vt are the state and observation noises respectively. For simplicity, we restrict ourselves to systems of the form:
xt = f(xt−1; Θ) + wt, wt ∼ N (0, Q), (3)
If the state space model is completely specified then recursive Bayesian inference can be applied to obtain an estimate of the latent states using the posterior p (x0:T |y1:T ) (Doucet et al., 2001). However in many applications, the parametric form of the state space model is unknown. While there exist methods that perform smoothing to obtain an estimate of x0:T (Barber, 2006; Fox et al., 2009; Djuric & Bugallo, 2006), we are often interested in not only obtaining an estimate of the continuous latent states but also in learning the dynamics f(·; Θ) that govern the dynamics of the system. In the simplest case, we can take a parametric approach to solving this joint state-parameter estimation problem. When f(·; Θ) and g(·; Ψ) are assumed to be linear functions, the posterior distribution over latent states is available in closed-form and the parameters can be learned via expectationmaximization. On the other hand, we have nonparametric methods that use Gaussian processes and neural networks to learn highly nonlinear dynamics and observations where the joint estimation is untractable and approximations are necessarily imployed (Zhao & Park, 2016; 2018; Frigola et al., 2014; Sussillo et al., 2016). Switching linear dynamical systems (SLDS) (Ackerson & Fu, 1970; Chang & Athans, 1978; Hamilton, 1990; Ghahramani & Hinton, 1996; Murphy, 1998) balance between these two extremes, approximating the dynamics by stochastically transitioning between a small number of linear regimes.
2Source code is available at https://github.com/catniplab/tree_structured_rslds

heading: 2.1 SWITCHING LINEAR DYNAMICAL SYSTEMS
text: SLDS approximate nonlinear dynamics by switching between a discrete set of linear regimes. An additional discrete latent state zt ∈ {1, . . . ,K} determines the linear dynamics at time t,
xt = xt−1 +Aztxt−1 + bzt + wt, wt ∼ N (0, Qzt) (4)
where Ak, Qk ∈ Rdx×dx and bk ∈ Rdx for k = 1, . . . ,K. Typically, zt is endowed with Markovian dynamics, Pr(zt|zt−1 = k) = πk. The conditionally linear dynamics allow for fast and efficient learning of the model and can utilize the learning tools developed for linear systems (Haykin, 2001). While SLDS can estimate the continuous latent states x0:T , the assumption of Markovian dynamics for the discrete latent states severely limits their generative capacity.

heading: 2.2 RECURRENT SWITCHING LINEAR DYNAMICAL SYSTEMS
text: Recurrent switching linear dynamical systems (rSLDS) (Linderman et al., 2017), also known as augmented SLDS (Barber, 2006), are an extension of SLDS where the transition density of the discrete latent state depends on the previous location in the continuous latent space
zt|xt−1, {R, r} ∼ πSB (νt) , (5) νt = Rxt−1 + r, (6)
where R ∈ RK−1×dx and r ∈ RK−1 represents hyperplanes. πSB : RK−1 → [0, 1]K maps from the reals to the probability simplex via stick-breaking:
πSB(ν) = ( π (1) SB(ν), · · · , π (K) SB (ν) ) , π (k) SB = σ(νk) ∏ j<k σ (−νj) , (7)
for k = 1, . . . ,K−1 and π(K)SB = ∏K−1 k=1 σ (−νk) where νk is the kth component of of ν and σ(ν) = (1 + e−ν)−1 is the logistic function (Fig. 1). By including this recurrence in the transition density of zt, the rSLDS partitions the latent space into K sections, where each section follows its own linear dynamics. It is through this combination of locally linear dynamical systems that the rSLDS approximates eq. (3); the partitioning of the space allows for a more interpretable visualization of the underlying dynamics.
Recurrent SLDS can be learned efficiently and in a fully Bayesian manner, and experiments empirically show that they are adept in modeling the underlying generative process in many cases. However, the stick breaking process used to partition the space poses problems for inference due to its dependence on the permutation of the discrete states {1, · · · ,K} (Linderman et al., 2017).

heading: 3 TREE-STRUCUTRED RECURRENT SWITCHING LINEAR DYNAMICAL SYSTEMS
text: Building upon the rSLDS, we propose the tree-structured recurrent switching linear dynamical system (TrSLDS). Rather than sequentially partitioning the latent space using stick breaking, we use a treestructured stick breaking procedure (Adams et al., 2010) to partition the space.
Let T denote a tree structure with a finite set of nodes { , 1, · · · , N}. Each node n has a parent node denoted by par(n) with the exception of the root node, , which has no parent. For simplicity, we initially restrict our scope to balanced binary trees where every internal node n is the parent of two children, left(n) and right(n). Let child(n) = {left(n), right(n)} denote the set of children for internal node n. Let Z ⊆ T denote the set of leaf nodes, which have no children. Let depth(n) denote the depth of a node n in the tree, with depth( ) = 0.
At time instant t, the discrete latent state zt is chosen by starting at the root node and traversing down the tree until one of theK leaf nodes are reached. The traversal is done through a sequence of left/right choices by the internal nodes. Unlike in standard regression trees where the choices are deterministic (Lakshminarayanan, 2016), we model the choices as random variables. The traversal through the tree can be described as a stick breaking process. We start at the root node with a unit-length stick π = 1, which we divide between its two children. The left child receives a fraction πleft( ) = σ(ν ) and the right child receives the remainder πright( ) = 1− σ(ν ) such that ν ∈ R specifies the left/right balance. This process is repeated recursively, subdividing πn into two pieces at each internal node until we reach the leaves of the tree (Fig. 1). The stick assigned to each node is thus,
πn =
{ σ(νpar(n)) I[n=left(par(n))] (1− σ(νpar(n)))I[n=right(par(n))] πpar(n) n 6= , 1 n = .
(8)
We incorporate this into the TrSLDS by allowing νn to be a function of the continuous latent state
νn(xt−1, Rn, rn) = R T nxt−1 + rn, (9)
where the parametersRn and rn specify a linear hyperplane in the continuous latent state space. As the continuous latent state xt−1 evolves, the left/right choices become more or less probable. This in turn changes the probability distribution πk(xt−1,Γ, T ) over the K leaf nodes, where Γ = {Rn, rn}n∈T . In the TrSLDS, these leaf nodes correspond to the discrete latent states of the model, such that for each leaf node k,
p (zt = k | xt−1,Γ, T ) = πk(xt−1,Γ, T ). (10) In general, the tree-structured stick-breaking is not restricted to balanced binary trees. We can allow more than two children through an ordered sequential stick-breaking at each level. In this sense, tree-structured stick-breaking is a strict generalization of stick-breaking. We also note that similar to rSLDS, the model can be made more flexible by introducing a dependence on the previous discrete latent in eq. (9) but for the rest of the paper, we stick to eq. (8).

heading: 3.1 A HIERARCHICAL DYNAMICS PRIOR THAT RESPECTS THE TREE STRUCTURE
text: Similar to standard rSLDS, the dynamics are conditionally linear given a leaf node zt. A priori, it is natural to expect that locally linear dynamics of nearby regions in the latent space are similar. Thus, in the context of tree-structured stick breaking, we impose that partitions that share a common parent should have similar dynamics. We explicitly model this by enforcing a hierarchical prior on the dynamics that respects the tree structure.
Let {An, bn} be the dynamics parameters associated with node n. Although the locally linear dynamics of a discrete state are specified by the leaf nodes, we introduce dynamics at the internal nodes as well. These internal dynamics serve as a link between the leaf node dynamics via a hierarchical prior,
vec([An, bn])| vec([Apar(n), bpar(n)]) ∼ N (vec([Apar(n), bpar(n)]),Σn), (11) where vec(·) is the vectorization operator. The prior on the root node is
vec ([A , b ]) ∼ N (0,Σ ) . (12) We impose the following constraint on the covariance matrix of the prior
Σn = λ depth(n)Σ , (13)
where λ ∈ (0, 1) is a hyper parameter that dictates how "close" a parent and child are to one another. The prior over the parameters can be written as, where the affine term and the vec(·) operator are dropped for compactness,
p({An}n∈T ) = p(A ) ∏
i∈child( )
p(Ai|A ) ∏
j∈child(i) p(Aj |Ai) . . . ∏ z∈Z p(Az|Apar(z)). (14)
It is through this hierarchical tree-structured prior that TrSLDS obtains a multi-scale view of the system. Parents are given the task of learning a higher level description of the dynamics over a larger region while children are tasked with learning the nuances of the dynamics. The use of hierarchical priors also allows for neighboring sections of latent space to share common underlying dynamics inherited from their parent. TrSLDS can be queried at different levels, where levels deeper in the tree provide more resolution.
TrSLDS shares some features with regression trees (Lakshminarayanan, 2016), even though regression trees are primarily used for standard, static regression problems. The biggest differences are that our tree-structured model has stochastic choices and the internal nodes contribute to smoothing across partitions through the corresponding hierarchical prior.
There are other hierarchical extensions of SLDS that have been proposed in the literature. In Stanculescu et al. (2014), they propose adding a layer to factorized SLDS where the top-level discrete latent variables determine the conditional distribution of zt, with no dependence on xt−1. While the tree-structured stick-breaking used in TrSLDS is also a hierarchy of discrete latent variables, the model proposed in Stanculescu et al. (2014) has no hierarchy of dynamics, preventing it from obtaining a multi-scale view of the dynamics. In Zoeter & Heskes (2003), the authors construct a tree of SLDSs where an SLDS with K possible discrete states is first fit. An SLDS with M discrete states is then fit to each of the K clusters of points. This process continues iteratively, building a hierarchical collection of SLDSs that allow for a multi-scale, low-dimensional representation of the observed data. While similar in spirit to TrSLDS, there are key differences between the two models. First, it is through the tree-structured prior that TrSLDS obtains a multi-scale view of the dynamics, thus we only need to fit one instantiation of TrSLDS; in contrast, they fit a separate SLDS for each node in the tree, which is computationally expensive. There is also no explicit probabilistic connection between the dynamics of a parent and child in Zoeter & Heskes (2003). We also note that TrSLDS aims to learn a multiscale view of the dynamics while Zoeter & Heskes (2003) focuses on smoothing, that is, they aim to learn a multi-scale view of the latent states corresponding to data but not suitable for forecasting.
In the next section we show an alternate view of TrSLDS which we will refer to as the residual model in which internal nodes do contribute to the dynamics. Nevertheless, this residual model will turn out to be equivalent to the TrSLDS.

heading: 3.2 RESIDUAL MODEL
text: Let {Ãn, b̃n} be the linear dynamics of node n and let path(n) = ( , . . . , n) be the sequence of nodes visited to arrive at node n. In contrast to TrSLDS, the dynamics for a leaf node are now determined by all the nodes in the tree:
p(xt|xt−1, Θ̃, zt) = N (xt|xt−1 + Āztxt−1 + b̄zt , Q̃zt), (15) Āzt = ∑
j∈path(zt)
Ãj , b̄zt = ∑
j∈path(zt)
b̃j , (16)
We model the dynamics to be independent a priori, where once again the vec(·) operator and the affine term aren’t shown for compactness,
p({Ãn}n∈T ) = ∏ n∈T p(Ãn), p(Ãn) = N (0, Σ̃n), (17)
where Σ̃n = λ̃depth(n)Σ̃ and λ̃ ∈ (0, 1). The residual model offers a different perspective of TrSLDS. The covariance matrix can be seen as representing how much of the dynamics a node is tasked with learning. The root node is given the broadest prior because it is present in eq. (16) for all leaf nodes; thus it is given the task of learning the global dynamics. The children then have to learn to explain the residuals of the root node. Nodes deeper in the tree become more associated with certain regions of the space, so they are tasked with learning more localized dynamics which is represented by the prior being more sharply centered on 0. The model ultimately learns a multi-scale view of the dynamics where the root node captures a coarse estimate of the system while lower nodes learn a much finer grained picture. We show that TrSLDS and residual model yield the same joint distribution (See A for the proof). Theorem 1. TrSLDS and the residual model are equivalent if the following conditions are true: A = Ã , An = ∑ j∈path(n) Ãj , Qz = Q̃z ∀z ∈ leaves(T ), Σ = Σ̃ and λ = λ̃

heading: 4 BAYESIAN INFERENCE
text: The linear dynamic matrices Θ, the hyperplanes Γ = {Rn, rn}n∈T \Z , the emission parameters Ψ, the continuous latent states x0:T and the discrete latent states z1:T must be inferred from the data. Under the Bayesian framework, this implies computing the posterior,
p (x0:T , z0:T ,Θ,Ψ,Γ|y1:T ) = p (x0:T , z1:T ,Θ,Ψ,Γ, y1:T )
p (y1:T ) . (18)
We perform fully Bayesian inference via Gibbs sampling (Brooks et al., 2011) to obtain samples from the posterior distribution described in eq. (18). To allow for fast and closed form conditional posteriors, we augment the model with Pólya-gamma auxiliary variables Polson et al. (2013).

heading: 4.1 PÓLYA-GAMMA AUGMENTATION
text: Consider a logistic regression from regressor xn ∈ Rdx to categorical distribution zn ∈ {0, 1}; the likelihood is
p(z1:N ) = N∏ n=1
( ex T nβ )zn
1 + ex T nβ
. (19)
If a Gaussian prior is placed on β then the model is non-conjugate and the posterior can’t be obtained in closed form. To circumvent this problem Polson et al. (2013) introduced a Pólya-Gamma (PG) augmentation scheme. This augmentation scheme is based on the following integral identity(
eψ )a
(1 + eψ) b = 2−beκψ ∫ ∞ 0 e− 1 2ωψ 2 p(ω)dω (20)
where κ = a− b/2 and ω ∼ PG(b, 0). Setting ψ = xTβ, it is evident that the integrand is a kernel for a Gaussian. Augmenting the model with PG axillary r.v.s {ωn}Nn=1, eq. (19) can be expressed as
p(z1:N ) = N∏ n=1
( ex T nβ )zn 1 + ex T nβ ∝ N∏ n=1 eκnψn ∫ ∞ 0 e− 1 2ωnψ 2 np(ωn)dωn = N∏ n=1 Eωn [e− 1 2 (ωnψ 2 n−2κnψn)].
(21) Conditioning on ωn, the posterior of β is
p(β|ω1:N , z1:N , x1:N ) ∝ p(β) N∏ n=1 e− 1 2 (ωnψ 2 n−2κnψn) (22)
where ψn = xTnβ and κn = zn − 12 . It can be shown that the conditional posterior of ωn is also PG where ωn|β, xn, zn ∼ PG(1, ψn) (Polson et al., 2013).

heading: 4.2 CONDITIONAL POSTERIORS
text: The structure of the model allows for closed form conditional posterior distributions that are easy to sample from. For clarity, the conditional posterior distributions for the TrSLDS are given below:
1. The linear dynamic parameters (Ak, bk) and state variance Qk of a leaf node k are conjugate with a Matrix Normal Inverse Wishart (MNIW) prior
p((Ak, bk), Qk|x0:T , z1:T ) ∝ p((Ak, bk), Qk) T∏ t=1 N (xt|xt−1+Aztxt−1+bzt , Qzt)1[zt=k].
2. The linear dynamic parameters of an internal node n are conditionally Gaussian given a Gaussian prior on (An, bn)
p((An, bn)|Θ−n) ∝ p((An, bn)|(Apar(n), bpar(n))) ∏
j∈child(n)
p((Aj , bj)|(An, bn)).
3. If we assume the observation model is linear and with additive white Gaussian noise then the emission parameters Ψ = {(C, d), S} are also conjugate with a MNIW prior
p((C, d), S|x1:T , y1:T ) ∝ p((C, d), S) T∏ t=1 N (yt|Cxt + d, S).
We can also handle Bernoulli observations through the use of Pólya-gamma augmentation. In the interest of space, the details are explained in Section B.1 in the Appendix.
4. The choice parameters are logistic regressions which follow from the conditional posterior
p (Γ|x0:T , z1:T ) ∝ p (Γ) T∏ t=1 p (zt|xt−1,Γ) = p (Γ) T∏ t=1 ∏ n∈path(zt)\Z (eνn,t) 1(left(n)∈path(zt)) 1 + eνn,t ,
where νn,t = RTnxt−1 + rn. The likelihood is of the same form as the left hand side of eq. (20), thus it is amenable to the PG augmentation. Let ωn,t be the auxiliary Pólya-gamma random variable introduced at time t for an internal node n. We can express the posterior over the hyperplane of an internal node n as: p((Rn, rn)|x0:T , z1:T , ωn,1:T ) ∝ p((Rn, rn)) T∏ t=1 N (νn,t|κn,t/ωn,t, 1/ωn,t)1(n∈path(zt)), (23) where κn,t = 121[j = left(n)]− 1 21[j = right(n)], j ∈ child(n). Augmenting the model with Pólya-gamma random variables allows for the posterior to be conditionally Gaussian under a Gaussian prior.
5. Conditioned on the discrete latent states, the continuous latent states are Gaussian. However, the presence of the tree-structured recurrence potentials ψ(xt−1, zt) introduced by eq. (10) destroys the Gaussinity of the conditional. When the model is augmented with PG random variables ωn,t, the augmented recurrence potential, ψ(xt−1, zt, ωn,t), becomes effectively Gaussian, allowing for the use of message passing for efficient sampling. Linderman et al. (2017) shows how to perform message-passing using the Pólya-gamma augmented recurrence potentials ψ(xt, zt, wn,t). In the interest of space, the details are explained in Section B.2 in the Appendix.
6. The discrete latent variables z1:T are conditionally independent given x1:T thus
p (zt = k|x1:T ,Θ,Γ) = p (xt|xt−1, θk) p (zt = k|xt−1,Γ)∑
l∈leaves(T ) p (xt|xt−1, θl) p (zt = l|xt−1,Γ) , k ∈ leaves(T ).
7. The conditional posterior of the Pólya-Gamma random variables are also Pólya-Gamma: ωn,t|zt, (Rn, rn), xt−1 ∼ PG(1, νn,t) .
Due to the complexity of the model, good initialization is critical for the Gibbs sampler to converge to a mode in a reasonable number of iterations. Details of the initialization procedure are contained in Section C in the Appendix.

heading: 5 EXPERIMENTS
text: We demonstrate the potential of the proposed model by testing it on a number of non-linear dynamical systems. The first, FitzHugh-Nagumo, is a common nonlinear system utilized throughout neuroscience to describe an action potential. We show that the proposed method can offer different angles of the system. We also compare our model with other approaches and show that we can achieve state of the art performance. We then move on to the Lorenz attractor, a chaotic nonlinear dynamical system, and show that the proposed model can once again break down the dynamics and offer an interesting perspective. Finally, we apply the proposed method on the data from Graf et al. (2011).

heading: 5.1 FITZHUGH-NAGUMO
text: The FitzHugh-Nagumo (FHN) model is a 2-dimensional reduction of the Hodgkin-Huxley model which is completely described by the following system of differential equations (Izhikevich, 2007):
v̇ = v − v 3
3 − w + Iext, τ ẇ = v + a− bw. (24)
We set the parameters to a = 0.7, b = 0.8, τ = 12.5, and Iext ∼ N (0.7, 0.04). We trained our model with 100 trajectories where the starting points were sampled uniformly from [−3, 3]2. Each of the trajectories consisted of 430 time points, where the last 30 time points of the trajectories were
used for testing. The observation model is linear and Gaussian where C = (
2 0 0 −2
) , d = [0.5, 0.5]
and S = 0.01I2 where In is an identity matrix of dimension n. We set the number of leaf nodes to be 4 and ran Gibbs for 1,000 samples; the last 50 samples were kept and we choose the sample that produced the highest log likelihood to produce Fig. 2 where the vector fields were produced using the mode of the conditional posteriors of the dynamics.
To quantitatively measure the predictive power of TrSLDS, we compute the k-step predictive mean squared error, MSEk, and its normalized version, R2k, on a test set where MSEk and R 2 k are defined as
MSEk = 1
T − k T−k∑ t=0 ‖yt+k − ŷt+k‖22 , R 2 k = 1− (T − k)MSEk∑T−k t=0 ‖yt+k − ȳ‖ 2 2 , (25)
where ȳ is the average of a trial and ŷt+k is the prediction at time t+ k which is obtained by (i) using the the samples produced by the sampler to obtain an estimate of x̂T given y1:T , (ii) propagate x̂T for k time steps forward to obtain x̂t+k and then (iii) obtain ŷt+k. We compare the model to LDS, SLDS and rSLDS for k = 1, . . . , 30 over the last 30 time steps for all 100 trajectories (Fig. 2I).
5.2 LORENZ ATTRACTOR
Lorenz attractors are chaotic systems whose nonlinear dynamics are defined by,
ẋ1 = σ (x2 − x1) , ẋ2 = x1(ρ− x3)− x2, ẋ3 = x1x2 − βx3.
The parameters were set to σ = 10, ρ = 28 and β = 8/3. The data consisted of 50 trajectories, each of length of 230 where the first 200 time points are used for training and the last 30 are used for testing. The observation model was a projection onto 10 dimensional space with Gaussian noise.We set the number of leaf nodes to be 4 and ran Gibbs for 1,000 samples; the last 50 samples were kept and we choose the sample that produced the highest log-likelihood to produce Fig. 3.
The butterfly shape of the Lorenz attractor lends itself to being roughly approximated by two 2- dimensional ellipsoids; this is exactly what TrSLDS learns in the second level of the tree. As is evident from Fig. 5B, the two ellipsoids don’t capture the nuances of the dynamics. Thus, the model partitions each of the ellipsoids to obtain a finer description. We can see that embedding the system with a hierarchical tree-structured prior allows for the children to build off its parent’s approximations.

heading: 5.3 NEURAL DATA
text: To validate the model and inference procedure, we used the neural spike train data recorded from the primary visual cortex of an anesthetized macaque monkey collected by Graf et al. (2011). The dataset is composed of short trials where the monkey viewed periodic temporal pattern of motions
of 72 orientations, each repeated 50 times. Dimensionality reduction of the dataset showed that for each orientation of the drifting grating stimulus, the neural response oscillates over time, but in a stimulus dependent geometry captured in 3-dimensions (Zhao & Park, 2017). We used 50 trials each from a subset of 4 stimulus orientations grouped in two (140 and 150 degrees vs. 230 and 240 degrees) where each trial contained 140 neurons. Out of the 140 neurons, we selected 63 well-tuned neurons. The spike trains were binarized with a 10 ms window for Bernoulli observation model and we truncated the onset and offset neural responses, resulting in 111 time bins per trial.
We fit TrSLDS with K = 4 leaf nodes and 3-dimensional continuous latent space; the sampler was run for 500 samples where the last sample was used to produce the results shown in Fig. 4. To obtain an initial estimate for x0:T , we smoothed the spike trains using a Gaussian kernel and performed probabilistic PCA on the smoothed spike trains.
From Fig. 4, it is evident that TrSLDS has learned a multi-scale view as expected. It is able to correctly distinguish between the two groups of orientations by assigning them to two different subtrees (green-yellow vs. red-orange). The leaf nodes of each subtree refines the periodic orbit further. From Fig. 4, we can see that TrSLDS also learns two limit cycles that are separated.

heading: 6 CONCLUSION
text: In this paper, we propose tree-structured recurrent switching linear dynamical systems (TrSLDS) which is an extension of rSLDS (Linderman et al., 2017). The system relies on the use of treestructured stick-breaking to partition the space. The tree-structured stick-breaking paradigm naturally lends itself to imposing a hierarchical prior on the dynamics that respects the tree structure. This tree-structured prior allows for a multi-scale view of the system where one can query at different levels of the tree to see different scales of the resolution. We also developed a fully Bayesian sampler, which leverages the Pólya-Gamma augmentation, to learn the parameters of the model and infer latent states. The two synthetic experiments show that TrSLDS can recover a multi-scale view of the system, where the resolution of the system increase as we delve deeper into the tree. The analysis on the real neural data verifies that TrSLDS can find a multi-scale structure.

heading: A PROOF OF THEOREM 1
text: Proof. Let T be a balanced binary tree with K leaf nodes. To show that the models are equal, it suffices to show the equivalence of the likelihood and the prior between models. For compactness, we drop the affine term and the vec(·) operator. The likelihood of TrSLDS is
p(x1:T |z1:T ,Θ) = T∏ t=1 N (xt|xt−1 +Aztxt−1, Qzt), (26)
and the likelihood of the residual model is
p(x1:T |z1:T , Θ̃) = T∏ t=1 N ( xt|xt−1 + Āztxt−1, Q̃zt ) . (27)
where Āzt is defined in eq. (16). Substituting Azt = ∑ j∈path(zt) Ãj into eq. (27) equates the likelihoods. All that is left to do is to show the equality of the priors.
We can express An = ∑ j∈path(n) Ãj recursively
An = Ãn +Apar(n). (28)
Plugging eq. (28) into ln p(An|Apar(n))
ln p(An|Apar(n)) = − 1
2
( An −Apar(n) )T Σ−1n ( An −Apar(n) ) + C (29)
= −1 2
( Ãn +Apar(n) −Apar(n) )T Σ−1n ( Ãn +Apar(n) −Apar(n) ) + C (30)
= −1 2 ÃTnΣ −1 n Ãn + C (31)
= −1 2 ÃTn
( λdepth(n)Σ )−1 Ãn + C (32)
where C is a constant. Because Σ = Σ̃ and λ = λ̃, eq. (32) is equivalent to the kernel of p(Ãn) implying that the priors are equal. Since this is true ∀n ∈ T , the joint distributions of the two models are the same.

heading: B DETAILS ON BAYESIAN INFERENCE
text: B.1 HANDLING BERNOULLI OBSERVATIONS
Suppose the observation of the system at time t follows
p(yt|xt,Ψ) = N∏ n=1 Bern(σ(υn,t)) = N∏ n=1 (eυn,t) yn,t 1 + eυn,t , (33)
υn,t = c T nxt + dn, (34)
where cn ∈ Rdx , dn ∈ R. Equation 33 is of the same form as the left hand side of eq. (20), thus it is amenable to PG augmentation. We introduce PG axillary variables ηn,t. Conditioning on η1:N , eq. (33) becomes
p(yt|xt, η1:N ) = N∏ n=1 e− 1 2 (ηn,tυn,t−2κn,tυn,t) (35)
∝ N∏ n=1 N (cTnxt + dn|κn,t/ηn,t, 1/ηn,t) (36) = N (Cxt +D|H−1t κt, H−1t ) (37)
where Ht = diag([η1,t, . . . , ηN,t]), κt = [κ1,t, . . . , κN,t] and κn,t = yn,t − 12 .
The observation is now effectively Gaussian and can be incorporated into the message passing for x1:T . The emission parameters are also conjugate with the augmented observation potential given a Matrix Normal prior. The conditional posterior on the axillary PG variables ηn,t also follows a PG distribution i.e. ηn,t|(cn, dn), xt ∼ PG(1, υn,t). Note that this augmentation scheme can also work for negative binomial, binomial, and multinomial observations (Polson et al., 2013; Linderman et al., 2015).
B.2 MESSAGE PASSING FOR x1:T
Assuming that the observations, y1:T , are linear and Gaussian, the posterior of the continuous latent states, x0:T , conditioned on all the other variables is proportional to
T∏ t=1 ψ(xt, xt−1, zt)ψ(zt, xt−1)ψ(xt, yt) (38)
where ψ(xt, xt−1, zt) is the potential of the conditionally linear dynamics, ψ(xt, yt) is the potential of the observation and ψ(xt−1, zt) is the recurrence potential. ψ(xt−1, zt) is a product of all the internal nodes traversed at time t
ψ(xt−1, zt) = ∏
n∈path(zt)\Z
ψn(xt−1, zt). (39)
If the potentials in eq. (38) were all linear and Gaussian, then we could efficiently sample from the posetrior of x0:T by passing messsages forward through Kalman Filtering and then sampling backwards; the prescence of the recurrence potentials prevent this because they aren’t Gaussian. By augmenting the model with the PG r.v.’s, the recurrence potential at internal node n becomes
ψn(xt−1, zt, wn,t−1) = N (RTnxt−1 + rn|κn,t−1/ωn,t−1, 1/ωn,t−1) (40)
which is effectively Gaussian , allowing for the use of the Kalman filter for message passing.
C INITIALIZATION
We initialized the Gibbs sampler using the following initialization procedure: (i) probabilistic PCA was performed on the data, y1:T to initialize the emission parameters, {C, d} and the continuous latent states, x1:T . (ii) To initialize the dynamics of the nodes ,Θ, and the hyperplanes, Γ, we propose greedily fitting the proposed model using MSE as the loss function. We first optimize over the root node
arg min A ,b
1
T T∑ t=0 ‖xt+1 − xt −A xt − b ‖22 , (41)
and obtain A∗ , b ∗ (Note that A ∗ , b ∗ can obtained in closed form by computing their corresponding OLS estimates). Fixing A∗ and b ∗ , we then optimize over the second level in the tree
arg min A1,b1,A2,b2,R ,r
1
T T∑ t=0 ‖xt+1 − σ(v )x̂1 − σ(−v )x̂2‖ , (42)
x̂i = xt + (A ∗ +Ai)xt + (b ∗ + bi) , (43)
v = R T xt + r . (44)
This procedure would continue until we reach the leaf nodes of the tree. Θ∗ and Γ∗ are then used to initialize the dynamics and the hyperplanes, respectively. In our simulations, we used stochastic gradient descent with momentum to perform the optimization. (iii) The discrete latent states, z1:T , were initialized by performing hard classification using Γ∗ and the initial estimate of x0:T .

heading: D DEALING WITH ROTATIONAL INVARIANCE
text: A well known problem with these types of model is it’s susceptibility to rotational and scaling transformation, thus we can only learn the dynamics up to an affine transformation Erosheva & Curtis (2017). During Gibbs sampling the parameters will continuously rotate and scale, which can slow
down the mixing of the chains. One possible solution to the issue is if we constrained C to have some special structure which would make the model identifiable; this would require sampling from manifolds which is usually inefficient. Similar to Geweke & Zhou (1996), we use the following procedure to prevent the samples from continuously rotating and scaling:
• Once we obtain a sample from the conditional posterior of the emission parameters {C,D}, we normalize the columns of C.
• RQ decomposition is performed on C to obtain U,O where U ∈ Rdy×dx is an upper triangular matrix and O ∈ Rdx×dx is an orthogonal matrix.
• We set C = U and rotate all the parameters of the model using O.

heading: E SCALABILITY AND COMPUTATIONAL COMPLEXITY OF THE INFERENCE
text: The rSLDS and the TrSLDS share the same linear time complexity for sampling the discrete and continuous states, and both models learn K-1 hyperplanes to weakly partition the space. Specifically, both models incur: an O(TK) cost for sampling the discrete states, which increases to O(TK2) if we allow Markovian dependencies between discrete states; an O(TD3) cost (D is the continuous state dimension) for sampling the continuous states, just like in a linear dynamical system; and anO(KD3) cost for sampling the hyperplanes. The only additional cost of the TrSLDS stems from the hierarchical prior on state dynamics. Unlike the rSLDS, we impose a tree-structured prior on the dynamics to encourage similar dynamics between nearby nodes in the tree. Rather than sampling K dynamics parameters, we need to sample 2K-1. Since they are all related via a tree-structured Gaussian graphical model, the cost of an exact sample isO(KD3) just as in the rSLDS, with the only difference being a constant factor of about 2. Thus, we obtain a multi-scale view of the underlying system with a negligible effect on the computational complexity.
To see how the number of discrete latent states effects the convergence speed of the Gibbs sampler, we fit 3 TrSLDS, with K = 2, 4, 8 respectively, to a Lorenz Attractor described in Sec. but used 250 trajectories to train the model as opposed to 50. To assess convergence, we plotted the logarithm of the joint density as a function of Gibbs samples. The results are shown Fig. 5.

heading: F SYNTHETIC NASCAR R©
text: We ran the TrSLDS on the synthetic NASCAR R© example from (Linderman et al., 2017), where the underlying model is an rSLDS. We trained TrSLDS on 10,000 time points and ran Gibbs for 1,000 samples; the last sample was used to create Fig. 6 where the vector fields where created from the mode of the conditional posterior of the dynamics. Even though the partitions were created using sequential stick-breaking, TrSLDS was able to reconstruct the dynamics of the model.

heading: G TREE SYNTHETIC NASCAR R©
text: To check whether the sampler is mixing adequately, we test TrSLDS on a twist on the synthetic NASCAR R© example where the underlying model is a TrSLDS. We also ran rSLDS on the example to highlight the limitations of seqeuntail stick-breaking. We trained both rSLDS and TrSLDS on 20,000 time points and Gibbs was ran for 1,000 samples; the last sample was used to create Fig. 7. To compare the precditive performance between the models, the R2k was computed for rSLDS and TrSLDS.

Review summaries:
The authors develop a tree structured extension to the recently proposed recurrent switching linear dynamical systems. Like switching linear dynamical systems (sLDS) the proposed models capture non-linear dynamics by switching between a collection of linear regimes. However, unlike SLDS, the transition between the regimes is a function of a latent tree as well as the preceding continuous latent state. Experiments on synthetic data as well as neural spike train data are presented to demonstrate the utility of the model.