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

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 models 1 . 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.

1. INTRODUCTION:

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 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 fully-Bayesian 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.

2. BACKGROUND:

Let x t ∈ R dx and y t ∈ R dy denote the latent state and the observation of the system at time t respectively. The system can be described using a state-space model:

EQUATION

where Θ denotes the dynamics parameters, Ψ denotes the emission (observation) parameters, and w t and v t are the state and observation noises respectively. For simplicity, we restrict ourselves to systems of the form:

EQUATION

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 (x 0:T |y 1: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 x 0: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.

(sequential) stick breaking tree-structured stick breaking 

2.1. SWITCHING LINEAR DYNAMICAL SYSTEMS:

SLDS approximate nonlinear dynamics by switching between a discrete set of linear regimes. An additional discrete latent state z t ∈ {1, . . . , K} determines the linear dynamics at time t,

EQUATION

where

A k , Q k ∈ R dx×dx and b k ∈ R dx for k = 1, . . . , K. Typically, z t is endowed with Markovian dynamics, Pr(z t |z t-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 x 0:T , the assumption of Markovian dynamics for the discrete latent states severely limits their generative capacity.

2.2. RECURRENT SWITCHING LINEAR DYNAMICAL SYSTEMS:

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

EQUATION

where R ∈ R K-1×dx and r ∈ R K-1 represents hyperplanes. π SB : R K-1 → [0, 1] K maps from the reals to the probability simplex via stick-breaking:

EQUATION

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 z t , 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) .

3. TREE-STRUCUTRED RECURRENT SWITCHING LINEAR DYNAMICAL SYSTEMS:

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 z t is chosen by starting at the root node and traversing down the tree until one of the K 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,

EQUATION

)

We incorporate this into the TrSLDS by allowing ν n to be a function of the continuous latent state

EQUATION

where the parameters R n and r n specify a linear hyperplane in the continuous latent state space. As the continuous latent state x t-1 evolves, the left/right choices become more or less probable. This in turn changes the probability distribution π k (x t-1 , Γ, T ) over the K leaf nodes, where Γ = {R n , r n } 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

(z t = k | x t-1 , Γ, T ) = π k (x t-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).

3.1. A HIERARCHICAL DYNAMICS PRIOR THAT RESPECTS THE TREE STRUCTURE:

Similar to standard rSLDS, the dynamics are conditionally linear given a leaf node z t . 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 {A n , b n } 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,

EQUATION

where vec(•) is the vectorization operator. The prior on the root node is vec

EQUATION

)

We impose the following constraint on the covariance matrix of the prior

EQUATION

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,

EQUATION

)

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 z t , with no dependence on x t-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.

3.2. RESIDUAL MODEL:

Let { Ãn , bn } 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:

EQUATION

We model the dynamics to be independent a priori, where once again the vec(•) operator and the affine term aren't shown for compactness,

EQUATION

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 n = j∈path(n) Ãj , Q z = Qz ∀z ∈ leaves(T ), Σ = Σ and λ = λ

4. BAYESIAN INFERENCE:

The linear dynamic matrices Θ, the hyperplanes Γ = {R n , r n } n∈T \Z , the emission parameters Ψ, the continuous latent states x 0:T and the discrete latent states z 1:T must be inferred from the data. Under the Bayesian framework, this implies computing the posterior, p (x 0:T , z 0:T , Θ, Ψ, Γ|y 1:T ) = p (x 0:T , z 1:T , Θ, Ψ, Γ, y 1:T ) p (y 1:T ) .

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) .

4.1 PÓLYA-GAMMA AUGMENTATION Consider a logistic regression from regressor x n ∈ R dx to categorical distribution z n ∈ {0, 1}; the likelihood is

EQUATION

)

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

EQUATION

where κ = a -b/2 and ω ∼ PG(b, 0). Setting ψ = x T β, it is evident that the integrand is a kernel for a Gaussian. Augmenting the model with PG axillary r.v.s {ω n } N n=1 , eq. ( 19) can be expressed as

p(z 1:N ) = N n=1 e x T n β zn 1 + e x T n β ∝ N n=1 e κnψn ∞ 0 e -1 2 ωnψ 2 n p(ω 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

EQUATION

where ψ n = x T n β and κ n = z n -1 2 . It can be shown that the conditional posterior of ω n is also PG where

ω n |β, x n , z n ∼ PG(1, ψ n ) (Polson et al., 2013).

4.2. CONDITIONAL POSTERIORS:

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 (A k , b k ) and state variance Q k of a leaf node k are conjugate with a Matrix Normal Inverse Wishart (MNIW) prior

p((A k , b k ), Q k |x 0:T , z 1:T ) ∝ p((A k , b k ), Q k ) T t=1 N (x t |x t-1 +A zt x t-1 +b zt , Q zt ) 1[zt=k] .

2. The linear dynamic parameters of an internal node n are conditionally Gaussian given a Gaussian prior on We can also handle Bernoulli observations through the use of Pólya-gamma augmentation.

(A n , b n ) p((A n , b n )|Θ -n ) ∝ p((A n , b n )|(A par(n) , b par(n) )) j∈child(n) p((A j , b j )|(A n , b n )).

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 (Γ|x 0:T , z 1:T ) ∝ p (Γ) T t=1 p (z t |x t-1 , Γ) = p (Γ) T t=1 n∈path(zt)\Z (e νn,t ) 1(left(n)∈path(zt))

1 + e νn,t , where ν n,t = R T n x t-1 + r n . 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((R n , r n )|x 0:T , z 1:T , ω n,1:T ) ∝ p((R n , r n )) T t=1 N (ν n,t |κ n,t /ω n,t , 1/ω n,t ) 1(n∈path(zt)) , (23) where κ n,t = 1 2 1[j = left(n)] -1 2 1[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 ψ(x t-1 , z t ) 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, ψ(x t-1 , z t , ω 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 ψ(x t , z t , w n,t ). In the interest of space, the details are explained in Section B.2 in the Appendix.

6. The discrete latent variables z 1:T are conditionally independent given x 1:T thus

p (z t = k|x 1:T , Θ, Γ) = p (x t |x t-1 , θ k ) p (z t = k|x t-1 , Γ) l∈leaves(T ) p (x t |x t-1 , θ l ) p (z t = l|x t-1 , Γ)

, k ∈ leaves(T ).

7. The conditional posterior of the Pólya-Gamma random variables are also Pólya-Gamma:

ω n,t |z t , (R n , r n ), x t-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.

5. EXPERIMENTS:

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) . 

5.1. FITZHUGH-NAGUMO:

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) :

EQUATION

)

We set the parameters to a = 0.7, b = 0.8, τ = 12.5, and I ext ∼ 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.01I 2 where I n 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, MSE k , and its normalized version, R 2 k , on a test set where MSE k and R 2 k are defined as

EQUATION

)

where ȳ is the average of a trial and ŷ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 xT given y 1:T , (ii) propagate xT for k time steps forward to obtain xt+k and then (iii) obtain ŷ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 ). Going one level lower in the tree, we see that in order to capture the nuances of the dynamics, each of the ellipsoids must be split in half. A trajectory from the Lorenz attractor is shown for comparison.

5.2. LORENZ ATTRACTOR:

(e) Plotting the dynamics, it is evident that the leaf nodes improve on it's parent's approximation. (f)

The R 2 k demonstrates the predictive power of TrSLDS.

Lorenz attractors are chaotic systems whose nonlinear dynamics are defined by,

ẋ1 = σ (x 2 -x 1 ) , ẋ2 = x 1 (ρ -x 3 ) -x 2 , ẋ3 = x 1 x 2 -βx 3 .

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 2dimensional 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.

5.3. NEURAL DATA:

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 x 0: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. The model is able to separate the limit cycles for each orientation group (green-yellow vs. red-orange) and refine them further with the leaf nodes. (right bottom) Two model generated predictive trajectories showing two stable limit cycles that resemble the two periodic orbits.

6. CONCLUSION:

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.

The observation is now effectively Gaussian and can be incorporated into the message passing for x 1: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 |(c n , d n ), x t ∼ 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 x 1:T

Assuming that the observations, y 1:T , are linear and Gaussian, the posterior of the continuous latent states, x 0:T , conditioned on all the other variables is proportional to

EQUATION

where ψ(x t , x t-1 , z t ) is the potential of the conditionally linear dynamics, ψ(x t , y t ) is the potential of the observation and ψ(x t-1 , z t ) is the recurrence potential. ψ(x t-1 , z t ) is a product of all the internal nodes traversed at time t

EQUATION

)

If the potentials in eq. ( 38) were all linear and Gaussian, then we could efficiently sample from the posetrior of x 0: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

EQUATION

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, y 1:T to initialize the emission parameters, {C, d} and the continuous latent states, x 1: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

EQUATION

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

EQUATION

EQUATION

EQUATION

)

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, z 1:T , were initialized by performing hard classification using Γ * and the initial estimate of x 0:T .

D DEALING WITH ROTATIONAL INVARIANCE:

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 ∈ R dy×dx is an upper triangular matrix and O ∈ R dx×dx is an orthogonal matrix.

• We set C = U and rotate all the parameters of the model using O. Figure 5 : The logarithm of the joint density was computed for all the samples generated from the 3 TrSLDS and smoothed using a trailing moving average filter. The sampler seems to converge to a mode rather quickly for all the three instantiations of the TrSLDS.

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(T K) cost for sampling the discrete states, which increases to O(T K 2 ) if we allow Markovian dependencies between discrete states; an O(T D 3 ) cost (D is the continuous state dimension) for sampling the continuous states, just like in a linear dynamical system; and anO(KD 3 ) 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 is O(KD 3 ) 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 .

F SYNTHETIC NASCAR R:

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. 

G TREE SYNTHETIC NASCAR R:

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 R 2 k was computed for rSLDS and TrSLDS. 

:

This 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).

Source code is available at https://github.com/catniplab/tree_structured_rslds

Figures and Tables:

Figure 1: State probability allocation through stick-breaking in standard rSLDS and the TrSLDS.

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|x 1:T , y 1:T ) ∝ p((C, d), S) T t=1 N (y t |Cx t + d, S).

Figure 2: TrSLDS applied to model the FitzHugh-Nagumo nonlinear oscillator. (a) The model was trained on 100 trajectories with random starting points. (b) The model can infer the latent trajectories. (c) The true vector field of FHN is shown where color of the arrow represents log-speed. The two nullclines are plotted in yellow and green. (d-f)The vector fields display the multi-scale view learned from the model where color of the arrows dictate log-speed The background color showcases the hierarchical partitioning learned by the model where the darker the color is, the higher the probability of ending up in that discrete state. As we go deeper in the tree, the resolution increases which is evident from the vector fields. (g) A deterministic trajectory from the leaf nodes (colored by most likely leaf node) with affine transformation onto a trajectory FHN (gray). (h) Plotting w and v over time, we see that the second level captures some of the oscillations but ultimately converges to a fixed point. The model learned by the leaf nodes captures the limit cycle accurately. (i) Performances compared for multi-step prediction. We see that TrSLDS outperforms rSLDS.

Figure 3: (a) The 50 trajectories used to train the model are plotted where the red "x" displays the starting point of the trajectory. (b) The inferred latent states are shown, colored by their discrete latent state. (c) We see that the second layer approximates the Lorenz attractor with 2 ellipsoids. A trajectory from the Lorenz attractor starting at the same initial point is shown for comparison. (d)Going one level lower in the tree, we see that in order to capture the nuances of the dynamics, each of the ellipsoids must be split in half. A trajectory from the Lorenz attractor is shown for comparison. (e) Plotting the dynamics, it is evident that the leaf nodes improve on it's parent's approximation. (f) The R 2 k demonstrates the predictive power of TrSLDS.

Figure 4: Modeling primary visual cortex spike trains. (top) Example spike raster plots in response to a drifting grating of orientations 150 and 240 degrees. Our data consisted of 200 such trials. (bottom) The average inferred latent trajectories over time for orientations 140 and 150 degrees colored by the most likely discrete latent state. (right top) Same plotted in space.The model is able to separate the limit cycles for each orientation group (green-yellow vs. red-orange) and refine them further with the leaf nodes. (right bottom) Two model generated predictive trajectories showing two stable limit cycles that resemble the two periodic orbits.

Figure 6: TrSLDS applied to the synthetic NASCAR R example. (left top) The true continuous latent states colored by their discrete state assignment. (left bottom) TrSLDS is able to infer the continuous and discrete latent states. (middle top) The dynamics were constructed such that the trajectories produced creates an oval track. (middle bottom) Although the partitioning is done through sequential stick-breaking, TrSLDS is still able to recover the dynamics. (right) We can see that TrSLDS can indeed recover a multi-scale view. The root nodes captures the rotation. The second level seperates the track into two rotations of different speeds.

Figure 7: TrSLDS and rSLDS applied to the tree version of to the synthetic NASCAR R . (left top)The true continuous latent states colored by their discrete state assignment. (left bottom) TrSLDS is able to infer the continuous and discrete latent states. (middle) TrSLDS was able to learn the underlying dynamics and partitions, indiicating that the sampler is mixing well. (right top)Due to the sequential nature of rSLDS, the model can't adequately learn the dynamics of the model. (right bottom) We can see from the k-step R 2 that TrSLDS outperforms rSLDS.

Peer review (iclr2019):
 This paper introduces a probabilistic model to model nonlinear dynamic systems with multiple granularities. The nonlinearity is achieved by using multiple local linear approximations. The method is an extension to rSLDS (recurrent switching linear dynamical systems), which in turn is an extension to SLDS. Pros: 1. Introducing the tree structure is a neat way of extending the existing rSLDS model to multiscale scenarios. 2. The paper is written clearly. The background is well illustrated and the idea rises naturally from there. The paper is also solid in the part describing the model. Con: 1. from the rSLDS paper (https://arxiv.org/pdf/1610.08466.pdf), the authors there was experimenting with some settings similar to those used in this paper. However, I am not able to find some explicit comparison between the TrSLDS and rSLDS in this work. I think it should be needed since TrSLDS itself is derived out from rSLDS, it would be good to show explicitly the advantage of the new model.