Spectral learning of Bernoulli linear dynamical systems models for decision-making
Authors: Iris R Stone
Yotam Sagiv
Il Memming Park
Jonathan W. Pillow
Venue: Accepted by TMLR ()

Abstract: Latent linear dynamical systems with Bernoulli observations provide a powerful modeling framework for identifying the temporal dynamics underlying binary time series data, which arise in a variety of contexts such as binary decision-making and discrete stochastic processes
such as binned neural spike trains. Here we develop a spectral learning method for fast, efficient fitting of probit-Bernoulli latent linear dynamical system (LDS) models. Our approach extends traditional subspace identification methods to the Bernoulli setting via a transformation of the first and second sample moments. This results in a robust, fixed-cost estimator that avoids the hazards of local optima and the long computation time of iterative fitting procedures like the expectation-maximization (EM) algorithm. In regimes where data is limited or assumptions about the statistical structure of the data are not met, we demonstrate that the spectral estimate provides a good initialization for Laplace-EM fitting. Finally, we show that the estimator provides substantial benefits to real world settings by analyzing data from mice performing a sensory decision-making task. 

Published in Transactions on Machine Learning Research (07/2023)
Spectral learning of Bernoulli linear dynamical systems mod- els
Iris Stonef Princeton Neuroscience Institute
istone@@princeton.edy
Princeton University
Yotam Sagiv
ysagiv@pprinceton..dd edu
Princeton Neuroscience Institute Princeton University
II Memming Park Champalimaud Foundation
memminng hatp:m/wwpprrrtttssrrrraeemmmmpp
Jonathan Pillow
pillou@princeton.com edu
Princeton Neuroscience Institute
Princeton University
Reviewed on OpenReview: https: // openreview. net/ form? id= /a
Abstract
Latent lincar dynamical systems with Bernoulli obsorvations provide powerful modeling framework for identifying the temporal dynamics underlying binary time series data, which arise in a variety contexts such as binary decision-making and discrete stochastic processes (e.g., binned neural spike trains). Here we develop a spectral learning method for fast, efficient fitting of Probit-Bernowd approach extends traditional subspace identification methods to the Bernoulli setting via transformation of the first and second sample moments. This results in robust, fixed-cost estimator that avoids the hazards of local optima and the long computation time of iterative fitting procedures like the expetaation-maxiiizattoo maximization (EM) algorithm. In regimes where data is limited or assumptions about the statistical structure of the data are not met, we demonstrate that the spectral estimate provides a good initialization for Laplace-EM fitting, Finally, we show that the estimator provides substantial benefits to real world settings by analyzing data from mice performing sensory drcision-making task.
latent lincar dynamical system (LDS) models. Our
1
Introduction
Latent. linear dynamical system (LDS) models are an important and widely-used tool for characterizing the structure many different types time series data, including natural language sequences (Belanger & Kakade 2015), task bisk-guided exploration (Wagenmaker et al. 2021), and neural population activity (Gao et al. 2015; Nonnenmacher al. 2017 Zoltowski et al. 2020). To fit these models, the typical approach to use inference method such the expectation maximization (FM) algorithm in order to find local maximum
of the likelihood function (Baum et al. 1970 Dempster et al. 1977 Shumway & Stoffer 1982). However, such inference methods are often tine-consuming, computationally expensive, and provide no guarantec of finding the global optimum. To address these issues, an active area of research has emerged using spectral methods to identify the parameters of lincar time invariant (LTT) systems from input-output data (Martens)
2010 Anandkunar et al. 2014 Belanger Kakade, 2015 Hazan et al. 2017 2018). Such methods rely
* Equual contribution

---

Published in Transactions on Machine Learning Research (07/2023)
an approach from control theory known as subspace identification (SSID), which usees moment based estimation to fit. the model parameters without the need for an iterative procedure (Ho & Kalman 1966
van Overschee & de Moor 1994 Viberg 1995 van Overschee & De Moor 1996 Qin 2006). The resulting
parameter estimates are consistent and avoid the problem of local optima that are common with inference
methods such as EM, although they lypically sulfer from lower accuracy. Therefore, a common strategy is to combine both approaches, using the SSID estimates as initializations for EM order to speed convergence and increase the likelihood of finding the global optimum.
Most recent developments in spectral estimation methods for LDS models have focused undriven (no inputs) lincar-Gaussia obsorvation processes. However, there has been work extending SSID for LDS models to handle Poisson observations, which is a useful framework for modeling neural spike-train data at certain timescales (e.g. > 10ms), as well as developing the general case for input-driven models (Buesing et al.) 2012). Yet. there need for additional research explicitly expanding spectral methods to other distribution classes. For example, binary time scries data are common in many fields, including reinforcement learning and decision-making. While there has been some work exploring second-order models for binary variables [Btthge & Berens 2008 Macke et al., 2009 Schein et al. (0166), there have been few attempts to develop a spectral lcarning method for LDS models with Bernoulli obscrvations. Such method would be extremely useful for identifying the dynamics of any binary time series data, including sequences of choice behavior and neural spike-trains with small bin sizes. In addition, recovery accuracy of the input-related parameters has not been well-characterized for non- non-Gaussian emissions models.
Here, we achieve two main innovations. First, we extend the SSID method to the probit-Bernould Bernoulli LDS case, deriving a new estimator: bestLDS (BErnoulli SpecTral Linear Dynamical System). From technical perspective, this involves overcoming difficulties due to redundancies in the moments that are not present in the Poisson or general cases. We show this method yields consistent estimates of the model parameters when fit to input-output data from a wide range of simulated datasets. Furthermore, we demonstrate that using the outputs from bestLDS as initializations for EM significantly accelcrates convergence. Second, we present new analyses, looking al parameter recovery error for the input-related LDS parameters as well as incorporating model-comparison metrics for the model as a whole. To our knowledge, neither of these analyses has been conducted on LDS models with non-Gaussian emissions. Lastly, we demonstrate the benefits of this method when applied to the behavior of mice performing perceptual decision-making task.
2
Background and related work
2.1 Applications of state space models in neuroscience
LDS and related state-space models have a long history in neuroscience. Bernoulli- and Poisson-LDS models have been particularly popular for inferring latent processes underlying neural spike trains (Gao et al. 2015 Nonnenmacher et al., 2017: Zoltowski et al. 2020 Valente et al. 2022). LDS models are also useful for linking neural dynamics to animal behavior (Linderman et al.) 2019), while discrete state-space models like hidden Markov models (HMMs) have become incrcasingly common tools for characterizing both neural (Escola et al.j
2011h and behavioral data (Wiltschko et al. 2015 Calhoun et al.
2019 Wiltschko al., 2020 . More recently,
there has been growing interest in using state-space models to describe behavior: (1) using continuous rather than discrete variables (Johnson et al. 2016
Costacurta et al.
2022);
(2) specifically in two-alternative
forced-choice (2AFC) decision- decision-making making contexts (Roy et al.] [021]] | and (3) while understanding the effects of
inputs such as sensory stimuli (Calhoun et al. 2019 Bolkan et al. 2022 Ashwood et al. 2022). Bernoulli-L.DS
models lie at the intersection of all three objectives, and thus their utility (along with methods that make their inference more efficient) will likely be central to neuroscience research in the future.
2.2 Extensions to LDS models
In addition to standard approaches, there has also been substantial work exlending LDS models to more flexible dynamical systems in order to better capture important features in neuroscience data. For example, Gaussian-process factor analysis (GPFA) was developed to provide insight into neural signals by extracting smooth, low-dimensional trajectories from recorded activity (Yu et al., 2009). Building on GPFA, Orthogonal stochastic linear mixing models (OSLMMs) relax the assumption that correlations across neurons are time-
2

---

Published in Transactions on Machine Learning Research (07/2023)
invariant and provide innovations in the regression framework that makes inference more tractable for large datasets (Meng & Bouchard 2022]. Other approaches employ LDS models constituents of multi-component frameworks in order to capture more complex temporal activity patterns in neural data. One class of examples inchudes latent factor analysis via dynamical systems (LFADS) and its extensions, which use recurrent neural networks (RNNs) to recover neural population dynamics (Pandarinath el al. 2018 Prince el al. 2021: Zhu
ct al.] 2022 Another approach employs variants of switching linear dynamical systems (sLDS), which learn flexible, non-linear models of neural and behavioral data by treating activity as sequences of repeated dynamical modes (Linderman et al., 2019: Nassar et al. 2019 Zoltowski et al.) 2020 Mudrik et al. 2022.
Wang et al.
2022: Weinreb et al. 2023). These advances make clear that LDS models serve valuable and
fundamental purpose in neuroscience and will continue to act as an essential building block for future work.
2.3 Inference methods
Fitting an LDS (as well as any of the extensions mentioned above) requires inferring the parameters governing the latent dynamics and emissions of the system. The most common approaches to inference are the
Expectation Maximization (EM) algorithm (Baum et al. 1970 Dempster et al. 1977 Shumway & Stoffer 1982) and Markov Chain Monte Carlo (MCMC) methods (Metropolis et al., 1953 Hastings 1970 Geman
& Geman 1984 Gelfand & Smith 1990). MCMC is a simulation method that produces samples that are
approximately distributed according to the posterior; these samples are then used to evaluate integrals once convergence is reached. One of the benefits of MCMC is that it tends to work even for complicated distributions. However, it is often difficult to assess accuracy and evaluate convergence. MCMC is also slow and tends to require a large number of samples. EM, on the other hand, is an iterative optimization technique that returns a local maximum of the likelihood. EM is advantageous for a number of reasons, including that the "M-step' equations (in which parameter values are updated) often exist in closed form, the value of the likelihood is guaranteed increase after each iteration of the algorithm, and in many cases EM requires fewer samples than MCMC. Several flexible variants of EM also exist. For example, when the conditional expectation of the log-likelihood is intractable, it is possible instead to use numerical approaches or to compute an analytical approximation using technique such as Laplace's method (Steele, 1996) or variational inference (Blei et al. 2017). In this case, optimization occurs over an 'Evidence Lower Bound' (ELBO) of the likelihood of the observed data. Despite these benefits, one major drawback of EM is that it is only guaranteed to find a local maximum. Thus, it often takes multiple initializations (or else a smart choice of initial parameters) to effectively find the global maximum. For this reason, methods such as bestLDS that identify good initializations for EM stand to greatly improve its computational efficiency. This is in contrast. to MCMC, wherein it is typical to discard the first several samples (due to being poor representations of the posterior distribution) and therefore the initialization scheme is less important.
2.4 Initialization approaches for EM
Given the impact that the choice of initialization has on EM performance, it no surprise that there is an extensive body of work on "smart" EM initialization methods for variety of state-space models, including
both LDS models (Buesing et al.] 2012 Hazan et al.
2017 2018) and HMMs (Siddiqi et al. 2010 Hsu et al.
2012
Anandkumar et al. 2014 Liu & Lemeire 2017
Mattila et al. 2017). Often, these techniques rely on a
combination of momentt-assee estimators (Martens 2010 Buesing et al. 2012 Hsu et al. 2012 Anandkumar et al. 2014; Mattila et al. 2017) and subspace identification methods (Ho & Kalman 1966 van Overschee &
de Moor. 1994 Viberg 1995 van Overschee & De Moor 1996 Andersson & Rydén 2009 Qin 2006). Despite
this work, significant gap has persisted in that these methods (I) often don't account for input-driven data and (2) do not work for binomial distributions. BestLDD addresses both of these shortcomings. Although there have been several developments for augmented approaches to EM in logistic models using Polya-Gamma latent variables (Polson et al. 2013 Schein et al.) 2016), these methods address other inference challenges in binomial-distributed data rather than improvements the initialization scheme. However, it will be an interesting case for future study whether bestLDS can be combined with these methods to further improve inference for logistic models.
3

---

Published in Transactions on Machine Learning Research (07 (07/202) /023)
(O Q)
dynamics noise
inputs
ut-1
Ut
U
B
B
latent state
Lt-1
Lt
T ++1
D
C
C
C
intermediate latents
Zt
Z
2 +++1
bernoulli observations
"O" or "1" choice
y --1
Y₁
y1+1
Figure 1: Generative model schematic.
3 Theory
3.1
Model description
Let denote the g-dimensional observation at time t. Concretely, we assume that Ye is emitted as Bernoulli noise from a latent dynamical system with p-dimensional linear-Gaussian dynamics, potentially driven at each time-step by an m-dimensional input uy. That is:
N(Ho, Qu)
- N(Ax,11 + Bu,, Q)
(1)
2t Cx. + Du,
V/ 21 ~ (ernoulli(f(z.))
where Po, Qo parameterize the initial distribution for the latent variable, A and B refer to the autoregressive and input-driven components of the mean of the latent dynamics, and Q is the covariance of the latent Gaussian noise. The quantity is latent convenience variable describing the state of the LDS in the g-dimensional output space as a linear transformation of the latent state by the loadings matrix C and the input-output interactions matrix D. Finally, the outputs are sampled from a Bernoulli distribution with mean f(z;) for an arbitrary function f (generally taken be either the logistic or probit function). (Figure shows a model schematic.)
3.2 Subspace identification
In general, subspace identification algorithms operate by first extracting estimate of the latent state and then using linear regression to identify the system parameters. Concretely, given the latent state sequence (I1, Ix) and the inputs (uı,..,, one may identify A and B by least-squares regression, using (I1, and (u,,.., -x to predict (I2,..,, Furthermore, given (x1, TN). (un,..,, uN), and the emissions (z, Ey), one can again use least-squares to compute C' and D.
In order achieve this, subspace identification algorithms typically operate on quantities called block Hankel matrices, consisting of tine-lagged subsequences of input ouput data. For example, in the
4

---

Published in Transactions on Machine Learning Research (07 (072))23)
case, the input and output block- Hankel matrices are defined as:
20
21
21
2/2
ZN
Zook
4p 2,)
11
12
22k
22h
40
u1
uN-
U1
1/2
UN
Uojk-1 UuI2K11/ Up
Up =
11p2-1
=
U2k 1
U2k
N
Note that each and above is typically (column) vector, these are block matrices. The top half of 20j28-1 and Uo (i.e., containing rows O through k - 1) represent the "past" relative to the kth row; for simplicity we denote these Zp and Up. respectively. The bottom half of each matrix, by contrast, represents the "future", which we denote by Z5 and Uy, respectively. Note that if N is the total number of data points, N' N 2k + is the number of columns of the matrix. Finally, the time lag parameter & may be called the Hankel size and is defined by the user with only the requirement that k 2 p (Katayama 2005). An in-depth description of this is beyond the scope of this article, but it can be shown that only the top p singular values of the Hankel matrix are greater than zero. Correspondingly, a typical heuristic is build the Hankel matrix with large initial k and then look at its singular value spectrum to pick smaller k more reflective of the underlying dynamics.
In order to see how the block-Hankel matrix relates to the system parameters, we repeatedly expand the linear- Gaussian components of Eq. to yield:
Zp I& (ro
-N-1) UpUp
D
CA
CB
D 0
CA2
Ta
V
CAR-2B
CB D
CAR-1
where F3 and V are referred to as the extended observability and controllability matrices, respectively. An analogous relationship holds for n and Uy.
Different subspace identification algorithms are characterized by their approach to decomposing the Hankel matrix in order to extract the system parameters. In case, we will use the N4SID algorithm (van Overschee & de Moor 1994). Briefly, this involves computing the RQ decomposition of the data matrix:
Up
U
RQ
(2)
2p
7,
It can be shown that processing various sub-blocks of the matrix R permits extraction of the latent state sequence or the bbservabilityy/oontolllbbbb matrices, either of which can then be used to recover each of the
system parameters (Ho & Kalman 1966 van Overschee & de Moor 1994 van Overschee & De Moor 1996
Qin 2006
1 Here overload from Eq. refer the output, can seen as the emission of linear-Gaussan LDS.
5

---

Published in Transactions on Machine Learning Research (07 (07/2023)
3.3 Moment conversion
In the Bernoulli case do not have accoss to {2t), the "intermediate cmissions" that arisc from the lincarly transformed latents and therefore cannot directly construct Zp Z described above. However, are able to infer their moments from the moments of Yp and Y5 (defined analogously), which has been shown in the general case to be sufficient proxy for use in SSID methods (Buesing et al. 2012).
In particular, let 2p denote (z....., =k- and 2 denote (zk: 22k-1)², and define Up uy analogously. Then the covariance of these concatenated naatix-valuue random variables has form that allows us to recover R in Eq. 2 via Cholesky decomposition:
Uhp
Σ = cov
uf
= RR
(3)
Zp
Here, we are exploiting the fact that while the observations are not Gaussian, the intermediate emissions and their time-lagged representations 2p and 2, Gaussian. This requires us to place assumption of normality assumption is violated (see Section 4.2 For now, we begin by describing how to convert the moments of in order to infer Σ.
the inputs, but we will show empirically that parameter recovery is good even when this
The core goal of the moment conversion process is to obtain the mean and covariance of (Eps zy) from the mean and covariance of (Wp: My). In particular, since Zp, 2fr Ups and uy are all jointly normal, is our goal find their joint mean and covariance given the moments of Up Mf. For convenience, we overload the notation Up us)
Zp and y 2
Mp Then the quantities we interested in ys
2
10 11
Duu y
Eus
E
(4)
1=
y=
We begin by relating the moments of to the mean μ² and the covariance " of 2. The required integrals in the AggisticBBernouli case are analytically intractable, therefore we instead consider the probit-Bernoull case. In this selting we have the implication y/: = + z4 2 0, and can therefore conclude:
Ely P(z : 2 0) 1 200 Eff)
(5)
E E[yay, = P(z 0 2 0) = 1 4 (0 p3 Eif 2 2i)
where refers to the ith entry of o, and & and 42 denote the univariate and bivariate cumulative normal distributions, respectively. Unfortunately, since the diagonal second moments are degenerate (that is, E[yay]] = E(y.),, this constitutes an underdetermined system of equations characterizing μ³² and Σ=z, To resolve this identifiability issue, we set the diagonal covariances E without loss of generality and then solve for po and Eij using a numerical root-finder.
In regimes with extremely high-autocorrelationn it is possible that the system in Eq. 5 will not yield solution. In such cases, it is possible to convert these constraints into minimization problem and search for the covariance that most closely matches the target. In practice, however, we expect that such regimes (where the vast majority of the data is comprised of either only Os or only 1s) are rare in real-world applications and did not encounter the need to so our analyses.
Next, we turn our attention the inputs. In our case, μ² and you are immediately available the empirical mean and covariance of the inputs. As such, it only remains to compute the cross-covariance Eas, To that.

---

Published in Transactions on Machine Learning Research (07 (07/202) /203)
end, evaluate the mross-moment
E Elyeu,
up(u, z)dzdu
(So up(u)z)du p(z)dz
(6)
(1th + yus (=- (0)))((p))) p())dz
- Dua u**) 100000 E+t)) + E"a Lo zy(z)d=
The integral in the second term the mean of truncated normal distribution and can be rendered Korotkov & Korotkov 20200:
2p(2)d== V=ß(erf(() - 1) + exp( -8²) 2a² /
with a 0.5 and B = 0..5 From this point, E is converted to R as described in Eq. 3 and given as an input to N4SID. summarize, the steps bestLDS are:
Compute the moments of
2. Convert those moments to moments of 2 as defined in Equation
(a) Compute μ and E by solving the system in Equation 5
(b) Use μ and = to compute nu by solving the system in Equation 6
Cholesky decompose
RR.
Take R as the input to the standard N4SID algorithm to recover the system matrices
4 Results
4.1 bestLDS infers correct parameters on simulated datasets
We begin by illustrating the effectiveness of the moment conversion for a very high- digh-diimnssooal system (observation dimension = 30, latent dimension = 15, Hankel parameter k = 20, and input dimension m = 3). We sampled from a probit-Bernolli LDS model and display the sample covariance the observations, cov(yt, II] the transformed covariance of the latents obtained by moment conversion, cov[z,, 2+1]] as well as the latent covariance matrix (Figure h)). Observe that the output covariance does not match the latent covariance, while the converted moments do.
Next, we demonstrate the relative insensitivity of our algorithm's performance to the Hankel size k by examining the singular value spectrum of the Hankel matrix, which one can use to determine the proper choice of the latent dimension p for independent analysis or prior to using a fitting method like EM (sce section 4.3 For simulated data N = 256, 000 in which = 5, we repeatedly fit bestLDS while varying the Hankel size (A: [3. 5, 7, 10)). The results show that the general structure of of the singular value spectrum is independent of (Figure 2p). Indeed for all four simulations only the first five singular values are non-zero, consistent with the true structure of the data.
To illustrate the effectiveness of bestLDDS, we next show that the estimator returns consistent, accurate estimates of the generative parameters on two simulated datasets, one dov-dimensional and the other high- dimensional. In both datasets, ground-truth data drawn by simulation from a probit-Bernouhe LDS with
For readability, we omit the subscripts and the right equals sign, but each is really 2; and each is really ",

---

Published in Transactions on Machine Learning Research (07 2022)
a
b
observed
transformed
singular values
1.0
7.5
k=3
0.5
singular value
k=5
5.0
k=7
0.0
k=10
2.5
-0.5
0.0
1.0
4
6
8
index
dataset A
dynamics matrix A
d
input-output interactions matrix D
gain matrix
0.6
1.8
avg. diff in EVs
avg. IAI error
0.02
IAI error
0.5
1.4
0.4
0.01
1.0
0.3
0.00
0.6
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
lo (N)
log1o
log
dataset B
f
dynamics matrix A
g
loadings matrix C h
input-output interactions matrix D
gain matrix
diff in EVs
0.04
subspace angle
1.0
avg. IAl error
0.06
avg. IAI error
5.0
0.04
0.02
0.5
2.5
0.02
0.0
0.00
0.0
0.00
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
lN)
l0g (N)
lo
logio(N)
true
best-LDS
2.5
-2.5
2.5
-2.5
2.5
0.0
-2.5
0
200
400
600
800
1000
time steps
Figure 2: Spectral estimator recovers simulation parameters in a variety of settings. (a) Covariance matrix Coviyt+ yt] showing the sample covariance the binary observations y (left), the covariance of obtained by converting the moments y (middle), and the ground truth covariance of (right) for a simulated dataset (g = 30, p = 15, k 20, m 3). Note that the transformed covariance closely matches the ground truth. (b) The singular value spectrum of the Hankel matrix computed at various settings of the Hankel size k. For all settings, only the top p= 5 singular values greater than zero. (c-e) Recovery error for bestLDS inferred parameters vs. training data size for low-dimensional dataset with = 3. For (c) we use the average absolute difference in the eigenvalues as the recovery error, and for (d) and (e) we use the average absolute elementwise difference. (f-i) Recovery error for bestLDS inferred parameters vs. training data size for a high digh-diimnsiooal dataset with = 10. For (g) the recovery error the angle of the subspaces spannod by C and Ć. (i) Latent state trajectorics simulated from true and inforred parameters for one of the simulations in dataset B. Both dalasets exhibited strong temporal correlations and small interaction terms between the latents/ latents/outputs outpuuts and the inputs. See Appendix for more details.

---

Published in Transactions on Machine Learning Research (07/2023)
Table 1: Comparison of the mean elementwise error in the gain matrix for inferred and true parameters across three different datasets. Data was simulated using a trial structure, with five total trials. Inference was performed on training sets composed of four trials aggregated together, and reported error values are the mcans across those training sets. valuee indicate the standard error of the mean. The number of data points listed in the table corresponds to the full size of the dataset (all five trials). See Appendix for details on simulation parameters.
Model Comparisons
Dataset A
Dataset B
Dataset C
Method
50k
256k
50k
256k
50k
256k
bestLDS
3.79 +.22
3.71 11.25
0.30 +0.03
0.19
#0.01
0.20 #0.01
0.17 #0.01
PLDS
N/A
N/A
2.37
#0.90
2.30
#0.28
1.17 #0.11
0.63 #0.03
gaussLDS
4.65 +1.25
3.90
11.16
2.45 10.02
2.25 #0.01
3.02 =0.25
2.24 10.12
No = O and Qu taken to be the stationary marginal covariance of the latents. A, B, and D are sampled independently from a standard normal distribution and then have their eigenvalues thresholded (spectral radius < 1 is needed for to cnsure stability). We take C to be a random orthonormal matrix such that, given A, the stationary marginal covariance of has unit diagonal (results of simulations where we do not make this assumption on the marginal covariance are available in Supplementary Figure A.1
Note that LDS models are not uniquely identifiable and SSID algorithms in general recover parameters only up to a similarity transformation Concretely, let x' = Tx for T an invertible matrix. Then we may recover A' = T-, = T B and C' = CT D is recovered canonically and thus we can look at the mean absolute clementwise error to establish the accuracy of our estimator for this matrix. For the other parameters, however, we must establish alternative error metrics. For A, we measure the error of its estimate À by examining the average absolute difference in their eigenvalues. For C, we look at the angle between the subspaces spanned by the true C' and the estimate Ć. For B, there is not simple error metric. However, we can examine the gain matrix G = C(I - A) + D, which is preserved canonically under the similarity transformation, laking the mean absolute difference across the inferred and true gain matrices as an overall indicator of accuracy in the recovery process. For both datasets, we computed these error metrics and took their average across 30 simulations for N E [1000, 256000]. The time taken to run bestLDS on these datasets is available in Supplementary Figure A.2
The accuracy of the system matrices recovered via bestLDS are displayed in Figure D-i.. In the low- dimensional dataset (dataset A), the errors for each relevant error metric decrease with increasing data size (note subspace angle is not possible to compute for C' when q < p). While the errors for A and G are non-zero even at .N- = 256000, systems in which < p (i.e. high-dimensional latent structure governs low-dimensional outputs), are a particular challenge for LDS models. Therefore it is notable that we nonetheless see decrcasing errors and for cerlain metrics quite small errors in this regime. As a rastingt example, we simulated a high-dimensional dataset with higher-dimension extremely well in this regime. For the system matrices A, C, and D well as the gain matrix G (Figure -iii, the error metrics decay to zero, indicating that estimator is consistent. Furthermore, the errors are quite low across the whole range of dataset-sizes taken example, the mean error in A at N = 1000 is less than 0.05. For completeness, we compared the accuracy of our estimator against running N4SID directly on 2 (i.c., pretending we have access to the Gaussian subset of the LDS and therefore forgoing the need to perform moment conversion) and found that bestLDS did not perform significantly worse on most variables (Supplementary Figure A.3). Thus our estimator is quite effective, even in a sparso-data. regime. To further validate our parameter recovery for dataset B, we simulated noiselessly using both the true and inferred parameters and found that the latent state traces overlapped almost exactly (Figure 23)
latent dynamics (dataset B). Our estimator performs
Lastly, we performed model comparisons on three datasets of different sizes to illustrate the accuracy of our recovery process in different seltings (Table D. Datasel B has the same simulation parameters as in Figure 2 Dataset C has lower autocorrelations, larger instantaneous correlations, and similarly small cant-laaent/outppir interactions. Datasct A * is the average over five datasets with the same generative parameters as dataset A; we did this to make our model comparison results more robust given the peculiarities of the q < p regime. Focusing on the gain matrix, as it is uniquely identifiable and incorporates the major
9

---

Published in Transactions on Machine Learning Research (07 (022))
random
gaussian
bestLDS
a
dataset (256k)
b
dataset D (50k)
C
dataset E
bits per sample
1.41
1.44
0.56
0.72
0
-0.56
-1.41
-1.69
-0.72
-2.82
-1.44
-2.22
10
25
40
55
70
10
25
40
55
10
25
40
55
EM steps
EM steps
EM steps
d
dataset F
dataset G
f
real data
bits per sample
0.36
-0.11
-0.80
0
-0.22
-1.00
-0.34
-0.36
-0.45
-1.20
-0.72
-0.56
10
25
40
55
70
5
10
15
20
5
10
15
20
EM steps
EM steps
EM steps
Figure 3: bestLDS is an effective initialization for EM. For panels a-d, the ELBO (converted to bits/sample) is plotted on each step of EM with random initializations, a bestLDS initialization, and a Gaussian initialization. Panels e-f are the same except that FM returns the log-evidence at each step. (a) Datasot is characterized by slow rotational dynamics and high-magnitude loadings matrix. It exhibits long stretches of Os and Is. Inputs were drawn from a zero-mean normal distribution with covariance Σ,, 10-41. (b) Same as (a) but with N = 50, 000. (c) Dataset E is characterized by extremely fast rotational dynamics. It exhibits switching behavior, alternating emissions of Os and 1s. Inputs were drawn from zero-mean normal distribution with covariance Σn = 0.11. (d) Dataset F is similar to dataset B but has inputs sampled from a Student-t distribution with three degrees of freedom. (e) Dataset G has one-dimensional observations
(q = 1) with slow rotational dynamics. See Appendix for more details on all simulation parameters. (f) Data taken from mice performing two-alternative forcod choice task in virtual reality. See Scction 4.3 for details.
parameters of the system, we computed the mean absolute elementwise error between the true generative parameters and three sets of inferred parameters: the bS inferred parameters, the inferred parameters using the Poisson LDS (pLDS) estimator as detailed in (Bucsing et al.) 2012). and the inferred parameters using a gain error amongst the benchmarks tested, and many instances those errors are extremely low. This indicates that there are conditions in which the bestLDS parameters can suffice on their own in capturing the characteristics of the system even without subsequent optimization. This also demonstrates the need for estimators that account for the specific distribution of the underlying data, as it is not sufficient to simply substitute different estimator and achieve accurate results.
(gaussLDS) estimator In the cases we have examined, bestLDS achieves the smallest
4.2 bestLDS provides good initializations for EM
In regimes where additional accuracy is desired, bestLDS may serve smart initialization method for EM. Indeed, estimates returned by spectral estimators have previously been shown to serve as good EM initialization points (Martens 2010 Buesing et al., 2012). This is particularly important given that EM fitting procedures typically require the user to initialize randomly (or at least, without complete knowledge of what good initialization would be) and therefore a single fit. is rarely sufficient determine that EM has
3 ccomputdd by treating the binary data real-valued and simply fitting Gaussian LDS to the data using N4SID (Lee. skipping the moment. conversion step)
10

---

Published in Transactions on Machine Learning Research (07/2023)
Table 2: Time to convergence for Laplace-EM for different initialization methods. For random initializations (N = 20), both total and mean (total / mean) times are reported. For bestLDS initializations, we report both the time to convergence and the time to acquire the initializations using the bestLDS estimator (estimator time + convergence time). All computations were performed an internal CPU cluster.
Convergence Times for Laplace- EM EM (min.
Initialization
Dataset D (256k)
Dataset D (50k)
Dataset E (256k)
random
902 / 45.10
232.40 / 11.62
571.75 / 28.09
Gaussian
38.67
11.60
39.67
best.DD
0.15 + 34.67
0.09 + 8.40
0.14 + 19.83
Initialization
Dataset F (100k)
Datasct G (256k)
Real Data (54883)
random
514.33 / 25.72
886.67 / /4.33
25.11 1.67
Gaussian
12.67
49.40
5.61
bestLDS
0.10 + 19.67
0.05 + 8.87
0.02 + 0.33
found the global optimum. Instead, users often fit EM to their data many times and take the parameters from those that produced the highest log-likelihoon or ELBO (Evidence Lower Bound). This can be extremely coneucins,,nng, especially when datasets are large. Thus, a method that can generate initializations that (1) converge quickly and (2) consistently converge to the global optimum can substantially reduce compute-time.
Accordingly, we examined how useful bestLDS is in such role. In this section assess the performance of this estimator as an initialization for Laplace-EM" an EM method using the Laplace approximation to efficiently compute the log posterior (Smith & Brown 2003; Kulkarni & Paninski 2007 Paninski et al. 2010 Yuan & Niranjan 2010). We then compare the results to random and "ruea--a initializations when fit to three different simulated datasets. For each dataset and initialization method, we have plotted the ELBO as function of the number of EM steps taken until the system reaches convergencc--thee this is a useful metric to examine since our primary interest training time and not performance. For the sake of clarity, we will briefly describe each dataset considered for this purpose and have provided full accounting of all the gencrative parameters in the Appendix.
First, we examined the case in which the dimensionality of the putputs was greater than that of the latents (i.e., q >> p). In particular, we considered datasets exhibiting properties that one might expect from typical Bernoulli datasots (c.g., mouse pocision-making behavior). For example, datasot D is characterized by long stretches of Os or 1s, which we might expect in behavioral datasets that exhibit relatively high autocorrelation (e.g., habitual behaviors, perseveration, etc.); this was implemented by imposing slow rotational dynamics. Here, the bestLDS initialization substantially outperforms Gaussian and random initializations, eliciting higher ELBOs from the first step of EM and converging to the optimal solution in fewer total steps (Figure 3b), thus saving notable computation time (Table 2 top, first column). This result also holds for the same dataset fit to fewer data points (Figure 3p). In this case, while the Gaussian initialization also achieves the global optimum,, it still takes more time to converge than the bestLDS initialization and in fact offers negligible same-savings over the average random initialization (Table 2 top, second column). Dataset E reflects a different extreme of Bernoulli-distribnted behavior in that it produces observations that rapidly switch between Os and 1s, which we implement by imposing fast rotational dynamics. We might expect this type of data from alternation tasks (i.e., the subject must make the opposite choice from the previous trial). The LestLDS initialization significantly aids performance on this dataset as well (Figure 3b) and reduces convergence time over Gaussian initializations by approximately 50% (Table 2 top, third column). Finally, dataset is almost identical to dataset B (see Figure 2 and Table D except. that inputs are drawn from a Student-t distribution with three degrees of freedom, thus testing the flexibility of our assumptions on the input sampling distribution. Here we find that both non-random initializations perform well, indicating the general utility of these spectral methods on input-output data even when the inputs are non- non-Gaussian (Figure and Table 2 bottom, first column).
"Implemented using the Python ssm package, licensed under the MIT license 5 uusnn the same process as described in the model comparisons section ffor the q 2 p regime, we say that EM has converged when the avg. |A] in the gain matrix successive steps is within tolerance (tol 0.15); for p < g we found that comparing the log-evidence on successive steps to be more reliable (tol 10, in bits/sample).
11

---

Published in Transactions on Machine Learning Research (07/2023)
Next, we turned our attention to the q < p case, specifically focusing on when q due to its practical relevance for neuroscience decision-making experiments, in which subjects are more likely to undergo tasks one-by-one and elicit a single choice output at time. Dataset G is a simulated dataset in this regime, characterized by slow rotational dynamics. As carlier, the bestLDS estimate attains a higher log-evidence than comparative methods from the first step of EM (Figure 3P) and converges much faster than the Gaussian or random initializations (roughly a factor of 5.5 faster; see Table P bottom, second column). Finally, we performed these comparisons on binary data from mice performing a decision-making task (described in further detail in the next section). Here, the bestLDS initialization converges in only three steps of EM (Figure 3h. All other initializations take significantly longer to converge (on the order of 17x slower for Gaussian LDS and 76x slower for the random initialization total time; see Table 2 bottom, third column) and also have initial log-evidence values that are noticeably lower than the bestLDS estimate.
4.3 Application to mouse decision-making data
To asses the performance of our estimator in real- world settings, we turned to a dataset from mice performing a sensory decision-making tas (Bolkan et al.) 2022). In this experiment, mice (n = 13) navigated a T-shaped maze in virtual reality. Mice were trained to detect visual cues that appeared to each side of their field of view while running down the main stem of the maze. At the end of the main stem, they had to turn toward the side of the maze with the most cues in order to receive a water reward (Figure Fa). The mice repeated this task for many trials in row, with an average daily session consisting of approximately 200 trials (N = 54, 883 total trials). On a random subset of 15% of trials within each session, the animals' striatal dopamine 1)2 receptor medium spiny neurons (MSNs) were inactivated optogenetically using a 532mm laser. This transient inhibition provides a mechanism for understanding how the mice perform the task when activity in their striatum, a brain region long associated with decision- decision-making and motor output (Cox & Witten 2019; Tang et al., 2022), is suppressed. In addition to its broad association with decision-making, the striatum is the chief larget of dopamine in the brain and has previously been demonstrated to encode reward prediction error signals (Schultz et al.) 1997). D2 MSNs in particular comprise cell-type in the striatum that has specifically been linked with suppression of movement and constitute the first locus in what is known as the "indirect pathway" of the striatum (see (Cox & Witten
2019 for a review).
When modeling the data from this task, we take the following as inputs u, to bestLDS, given their high likelihood as relevant predictors in this task: 1) the difference between the number of right and left cues on each trial, 2) the animal's previous choice rewarded (coded in the model as #1 if the mouse made a correct right/left choice and 0 otherwise), and 3) the presence of inhibition (coded +1 if delivered to the tight/left hemisphere of the brain and 0 otherwise). We also consider including n-th order previous choices additional inputs (coded as +1 for the /left choice the animal made trials in the past, regardless of reward status). For our preliminary analysis, we vary the total number of previous choices incorporated in the model for n E (0, 1,2,3,4), yielding m E (3, 4,6,6,7, 7)
To determine our choice of p, we first examine the singular value spectrum of the transformed data (Fig- ure 4p). For this dataset, increasing p yields only marginal decreases in the singular value after p 4. We correspondingly take p- 4 as an upper-bound and conduct further analysis within the range p [2, 4), while also varying the number of inputs as previously described. In particular, we fit be to the data and looked at how well the inferred parameters predict subsequent choice (Figure --d) for each combination of and m. For all settings of m, the inferred parameters with p = 3 outperform the other p settings and predict choice well (> 70%). Of these, the best- cest-performing setting is with m = 4 (1st-order previous choico). Therefore, for all remaining analyses we fit bestLDS with m = and
3
It's worth noting that this approach for specifying a priori the dimensionality of the system showcases an additional benefit of bestLDS. Without our method, the typical approach would be to fit many variants of the model using standard inference procedure like EM, choosing the state number by some comparative method (e.g., cross-validation). However, this tends to be quite time-consuming and computationally costly. Here, the singular value spectrum combined with bestLDS fitting immediately gives an informative summary
Data obtained with permission directly
8 Note that because g ==1, 1, this places in the q < p regime, which can sometimes suffer from performance issues despite performing well here; see Figure g
12

---

Published in Transactions on Machine Learning Research (07 (07/2023)
a
d
true x bestLDS
singular value
0.3
m=
0.705
Ne0 000000000000
100000000 ool obk Jeo0o
% correct
choice
9p
0.2
m=
0.700
laser
2=
90
0.1
1=3
pp
9
0.0
0.695
p=4
100
1000
3000000
0
2
6
8 10
4
5
6
7
10
20
30
40
50
index
number of inputs
trials
perseverative
gaussLDS
glm
recovered latents
bestLDS
bestLDS+EM
2.5
0.0
fraction correct
0.75
log- log-evidence
-6000
latent
bioun 8.5
-2.5
0.65
-7000
0.55
-8000
www
0.45
-9000
100
200
300
400
500
model
model
trials
recovered parameters
latent state dynamics
input-laeet interactions
loadings
input-output interactions
1.00
0.4
0.4
0.2
0.50
0.2
0.1
0.2
0.00
0.0
0.0
0.0
-0.50
-0.2
-0.1
-0.2
-0.2
-1.00
-0.4
-0.4
impulse response plots
A cues
prev. choice
prev. rew. choice
inhibition
0.05
0.4
0.10
0.02
-0.05
0.2
0.05
0.01
-0.15
0.0
0.00
0.00
-0.25
0
20
40
0
20
40
0
20
40
20 trials
trials
trials
trials
Figure 4: Analysis of mouse binary dccision-making data with bestLDS. (a) Task schematic. (b) Singular value spectrum of the Hankel matrix with increasing number of previous choices used as inpuls. (c) Prodiction accuracy of bestLDS for various settings and m. (d) Example stretches of predicted and true choices. (e) 5-fold vali--aaldd performance comparisons for variety of methods on this dataset. Left: fraction of choices correctly predicted. Right: test set log-evidence. Black dots indicate each of the individual folds. (f) An example stretch of recovered latents. (g) The fitted parameters estimated by bestLDS + EM. (h) Impulse responses for the fitted system. In each subplot, data has been simulated noiselessly with unit input in the indicated dimension.
of the parameters of the dataset, providing in mere minutes what would otherwise take many hours and hundreds of different EM initializations (that is, dozens of initializations per parameter setting and value). Our spectral estimator can even be applied for this purpose when the ultimate goal is to fit to a different state-space model, such as discrcto-state HMMs.
To further assess bestLDS performance on the real data, we use five-fold cross-validation compare its performance to three other models: simple perseverative model, which predicts a high probability (70%) of making the same choice the previous trial: Gaussian LD and Bernoulli generalized lincar
We convert the real- vell-valudd observations produced by the Gaussian LDS binary observations by taking y 0f < and jt 2, =
13

---

Published in Transactions on Machine Learning Research (07/2023)
model (GLM; Figure 4p). Specifically, we operationalize "performance" in two ways: fraction of choices in the test. set correctly predicted and the log-evidence of the data evaluated under the fitted parameters. Under both metrics, bestLDS significantly outperforms the perserverative and Gaussian LDS models. Furthermore, four out of five test sets achieve equal or greater performance than the GLM on the prediction accuracy metric, though this is not borne out in the log-evidences. This already notable, given that bestLDS only returns an estimate of the system parameters for a multi-state LDS, whereas the GLM is a regression model that is capable of finding fairly accurate weights for a dataset of this size. However, when coupled with EM, which converges in just three steps and takes less than a minute to run, find that the resultant inferred parameters perform better than all other methods regardless of metric.
Not only does bestLDS perform well against comparative models, but it also offers unique scientific insights into the data that don't require subsequently using the inferred parameters as EM initializations. For example, inspection of the recovered latents x shows that one of the latent dimensions dominates the tendency of the overall response; the inferred traces on this latent dimension exhibit significantly higher magnitude than the other two (Figure [), and the corresponding weights in the dynamics matrix A are also of high magnit ude (Figure 4g. left column). As such, "activity" in this dimension can be described as strongly driving both the decisions of the mouse and also the latent state of the LDS underlying it. We can examine how this is affected by the inputs by examining imput-latent interactions matrix B; here we see that the largest magnitude elements of the matrix are those mapping the A cues and previous choice inputs onto the high-magnitude dimension. Overall, this suggests that the first latent dimension is mostly driven by autocorrelation (since its self-weight in A is high and the activity of the other dimensions low), with strong influence due to those inputs, and activity in this dimension largely dominates the contribution from the latent state to the outputs.
Furthermore, examining the input-output interactions matrix reveals the general effect that each input has on choice; more cues on the right side of the maze, previous rightward choice, and previous rightward rewarded choice are all more likely to elicit a rightward choice on the current trial (since they have positive weights, and right choices are coded as 1), whereas hegtshhemipperre inhibition is more likely to clicit a leftward choice (since il has a negative weight). Note that the results obtained from using the bestLDS-inffrred parameters as an initialization for EM produce the same qualitative interpretation of the system, indicating that for this dataset it is possible to accurately characterize the relationship between the inputs, latent states, and animals' choices solely from bestLDS without using an additional inference procedure (see Supplementary Figure A.4) It's also worth noting that it would not be straightforward to reach the same conclusions by assuming linear-Gaussian observations and using standard N4SID methods. In our investigations, such an approach clicits uninterpretable system parameters, including large eigenvalues in A and pathological impulse-response traces that diverge.
A convenient way to analyze the inferred parameters to examine their impulse responses that is, simulating noiseless Bernoulli LDS with the inferred parameters starting with unit input in single dimension at time step 0 (Figure 4n). The resultant responses can be seen as a general summary of how each input affects the output state of the system-and and have been shown to be fully sufficient to characterize LTI systems. For example, a +1 input in the A cues dimension (more right cues than left) makes the system more likely to emit a right-turn, as expected given the task structure. The second trace shows that the mice perseverate (i.e., are, on average, more likely to repeat their previous action than to do the opposite one), as a previous right choice more likely causes a right turn on the current trial. The same effect holds when the previous choice was rewarded, but with smaller magnitude, suggesting that perseveration is not strongly contingent on reward. The negative response in the final trace indicates that inhibiting neural activity in one hemisphere encourages the mice to turn contralaterally (i.c., inhibition of the striatum on the right hemisphere will influence the mice to turn left). Also note the relatively long decay times for choice-related inputs, which speak to the perseverative aspect of the dynamics. These conclusions largely match previously observed patterns in this data using input-output HMMs, and the performance of best +DM in predicting choice in this dataset is similar to that found in the discrete-state case (Bolkan et al., 2022). These findings raise interesting questions about the relationship between discrete- and continuous-state models for decision-making behavior that would be a promising avenue for future study.
14

---

Published in Transactions on Machine Learning Research (07 2023))
5 Discussion
We have presented spectral estimator for driven Bernoulli latent linear dynamical systems. On a variety of simulated datasets, bestLDS returns consistent, high-quality estimates of the generative parameters efficiently on its own and also as an initializer for the rnpectationmmaxiizaatte robust the large data, high-dimensional regime, which is an especially desirable use case of the model given the potential time-savings in fitting. One notable limitation of the method is that accuracy of the estimates may suffer in low-dimensional regimes where the number of observations is smaller than the number of latent dimensions, however show in Section 4.3 that predictive performance is still strong and the returned estimates yield novel scientifie insights, in addition to the significant tine-savings over traditional inference methods (i.e. EM). These results can then inform subsequent analyses and model fitting procedures that require longer compute time (c.g., an asscssment of the optimal number of latent dimensions and the importance different inputs on the latent dynamics). Identifying the syslem matrices the continuous state case may also provide relevant information for discrete state-space models. This would be another interesting area of future work, given the popularity of such models of behavior in neuroscience (Bolkan et al, 2022 Ashwood et al. )022).
algorithm. Estimates are particularly
In the extremely-high dimensional regime, the primary limitation is the need to numerically infer the moments of z. Since the number of parameters increases as ((kg)²), time spent on moment conversion grows quickly with the dimensionality of the data. However, the compute time required for moment conversion is unlikely negate the time-savings garnered from faster EM convergence.
In sum, bestLDS proves an efficient estimator for driven Bernoulli-LDS data. Parameter recovery even without EM is efficient and accurate, potentially averting the need to spend time rumning iterations of expensive search algorithms. Coupled with EM, the method greatly speeds convergence of the final estimate, even in cases where normality of the inputs is violated. Binary time-scries data arise in many contexts, including reinforcement learning and decision-makinng, weather outcomes, finance, and neuroscience. We thus expect that bestLDS will be broadly useful to practitioners looking to quickly acquire a description of the latent dynamics underlying these systems.
Data availability
The data that support the findings of this study are publicly available on figshare at http:://ggggaaaa..mmmm//rddd///s
associated (aata23355670))
Code availability
Code for general use applications of bestLDS analyses developed in this study, including all ap- plications to simulated and real data presented in the manuscript, are available on GitHub http:://gihubb.co///iistoon////
Acknowledgments
LR.S. was supported by the National Institute of Mental Health of the N.L.H. under award number F31MH131304. T.M.P was supported by N.I.H. BRAIN initiative grant (0R01DA056400400)) J.W.P. was supported by grants from the Simons Collaboration on the Global Brain (SCGB AWD543027), the N.I.H. BRAIN initiative (PROIDA056404-0)), and a U19 N.N..N.NN..SS. BRAIN Initiative Award (5U19NS104648),
References
Anima Anandkumar, Rong Ge, Daniel Hsu, Sham M. Kakade, and Matus Telgarsky. Tensor decompositions for learning latent variable models. Journal of Machine Learning Research, 15:2773 2832, 2014.
Sofia Andersson and Tobias Rydén. Subspace estimation and prediction methods for hidden Markov models. The Annals of Statistics, 37(6B):41314411 4152,, 2009.
15

---

Published in Transactions on Machine Learning Research (07/2023)
Zoe C. Ashwood, Nicholas A. Roy, Iris R. Stone, Anne E. Urai, Anne K. Churchland, Alexandre Pouget, and Jonathan W. Pillow. Mice alternate between discrete strategies during perceptual decision-making, Nalure Neuroscience, 25(2):201 212, 2022.
Leonard E. Baum, Ted Petrie, George Soules, and Norman Weiss. A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. The Annals of Mathematical Stalislics, (1):164 171, 1970.
David Belanger and Sham Kakade. A linear dynamical system model for text. In Proceedings of the 32nd International Conference on Machine Learning, 2015.
Matthias Bethge and Philipp Berens. Near-maximum entropy models for binary neural reprosentations of natural images. In Advances in Neural Informalion Processing Systems, volume 20, pp. 97- 97-104, 2008
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: review for statisticians.
Journal of the American Statistical Association, 112((18):859 877, 2017.
Scott S. Bolkan, Iris R. Stone, Lucas Pinto, Zoe C. Ashwood, Jorge M. Iravedra Garcia, Alison L. Herman, Priyanka Singh, Akhil Bandi, Julia Cox, Christopher A. Zimmerman, Jounhong Ryan Cho, Ben Engelhard, Jonathan W. Pillow, and Ilana B. Witten. Opponent control of behavior by dorsomedial striatal pathways depends on task demands and internal state. Nature Neuroscience, 25(3):345-355 2022.
Lars Buesing, Maneesh Sahani, and Jakob H Macke. Spectral learning of linear dynamics from generalised- lincar observations with application to neural population data. In Advonces in Neural Information Processing Systems, volume 25, 2012.
Adam Calhoun, Jonathan W. Pillow, and Mala Murthy. Unsupervised identification of the internal states that shape natural behavior. Nature Nearoscience, 22(12):2040-20 2049, 2019.
Julia C. Costacurta, Lea Duncker, Blue Sheffer, Winthrop Gillis, Caleb Weinreb, Jeffrey Evan Markowitz, Sandeep R. Datta, Alex H. Williams, and Scott Linderman. Distinguishing discrete and continuous behavioral variability using warped autoregressive HMMs. In Adconces in Neural Information Processing Systems, volume 36, 2022.
Julia Cox and Ilana B. Witten. Striatal circuits for reward learning and decision decision-making, Nature Reviews Neuroscience, 20(8):482-495 2019.
A. P. Dempster, N. M. Laird, and D. B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B (Methoolologica)) 39(1):1 22, 1977.
Scan Escola, Alfredo Fontanini, Don Katz, and Liam Paninski. Hidden Markov models for the stimulus- response relationships of multistate neural systems. Neural Compulation, 23(5):1071 1132, 2011.
Yuanjun Gao, Lars Busing, Krishna Shenoy, and John P Cunningham. High-dimensional neural spike train analysis with generalized count linear dynamical systems. In Advances Neural Information Processing Systems, volume 28, 2015.
Alan E. Gelfand and Adrian F. M. Smith. Sampling-base based approaches to calculating marginal densities. Journal of the American Statistical Association, 35(410):398-444 409, 1990.
Stuart Geman and Donald Geman. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Transactions on Pattern Analysis and Machine Intelligence, PAMM-6666:72117441 1984.
W. K. Hastings. Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57 (1):97 109, 1970.
Elad Hazan, Karan Singh, and Cyril Zhang. Learning linear dynamical systems via spectral filtering. In A Aduances in Neural Information Processing Systems, volume 30, 2017.
16

---

Published in Transactions on Machine Learning Research (07/2023)
Elad Hazan, Holden Lee, Karan Singh, Cyril Zhang, and Yi Zhang. Spectral filtering for general linear dynamical systems. In Advances in Neural Information Processing Systems, volume 32, 2018.
B. L. Ho and R. E. Kalman. Editorial: Effective construction of linear state-variable models from input/ / functions. Antomalisieraagsteehnni, 14(1-12):545-548, 1966.
Daniel Hsu, Sham M. Kakade, and Tong Zhang. A spectral algorithm for learning hidden Markov models. Journal of Computer and System Sciences, 78(5):1460 1480, 2012.
Matthew J Johnson, David K Duvenaud, Alex Wiltschko, Ryan P Adams, and Sandeep R Datta. Composing graphical models with neural networks for structured representations and fast inference. In Advances in Neural Information Processing Systems, volume 29, pp. 2946 2954, 2016.
Tohru Katayama, Subspace Methods for System Identification. Communications and Control Engineering. Springer, Berlin; London, 2005. ISBN 978-1-85533-981-44
Nikolai E. Korotkov and Alexander N. Korotkov. Integrals Related to the Error Function. Chapman and Mall/CRC, New York, 2020. ISBN 978-0-367-509223-..
Jayant E. Kulkarni and Liam Paninski. Common-input models for multiple neural spike-train data. Network: Computation in Neural Systems, 18(4):375 407, 2007.
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. elegans. bioRwie, 2019.
Tingting Liu and Jan Lemeire. Efficient and effective learning of HMMs based on identification of hidden states. Mathematical Problems in Engineering, 2017:e7318940, 2017.
Jakob H. Macke, Philipp Berens, Alexander S. Ecker, Andreas S. Tolias, and Matthias Bethge. Generating spike trains with specified correlation coefficients. Neural Computation, 21(2):397- 423, 2009.
James Martens. Learning the linear dynamical system with ASOS. In Proceedings of the 27th International Conference on Machine Learning, 2010.
Robert Mattila, Cristian R. Rojas, Vikram Krishnamurthy, and Bo Wahlberg. Identification of hidden Markov models using spectral learning with likelihood maximization. In IEFE 56th Annual Conference on Decision and Control, pp. 5859 5864. Institute of Electrical and Electronics Engineers (IEEE), 2017.
Rui Meng and Kristofer Bouchard. Bayesian inference in high-dimensional time-series with the orthogonal stochastic linear mixing model. arXiv, 2022.
Nicholas Metropolis, Arianna W. Rosenbluth, Marshall N. Rosenbluth, Augusta H. Teller, and Edward Teller. Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21:1087-1092, 1953.
Noga Mudrik, Yenho Chen, Eva Yezerets, Christopher J. Rozell, and Adam S. Charles. Decomposed linear dynamical systems (dLDS) for learning the latent components of neural dynamics. arXiv, 2022.
Josue Nassar, Scott W. Linderman, Monica Bugallo, and II Memming Park. Tree-structured recurrent switching lincar dynamical systems for multi-scale modeling. In 7th International Conference on Learning Representations, 2019.
Marcel Nonnenmacher, Srinivas Turaga, and Jakob Macke. Extracting low-dimensional dynamics from multiple large-scale neural population recordings by learning to predict correlations. In Advances in Neural Information Processing Systems, volume 21, 2017.
Chethan Pandarinath, Daniel J. O'Shea, Jasmine Collins, Rafal Jozefowicz, Sergey D. Stavisky, Jonathan C. Kao, Eric M. Trautmann, Matthew T. Kaufman, Stephen I. Ryu, Leigh R. Hochberg, Jaimic M. Henderson, Krishna V. Shenoy, L. F. Abbott, and David Sussillo. Inferring single-trial neural population dynamics using sequential ial auto-encoders. Nature Methodss 15(10):805 815, 2018.
17

---

Published in Transactions on Machine Learning Research (07/2023)
Liam Paninski, Yashar Ahmadian, Daniel Gil Ferreira, Shinsuke Koyama, Kamiar Rahnama Rad, Michael Vidne, Joshua Vogelstein, and Wei Wu. A new look at state-space models for neural data. Journal of Computational Neuroscience, 29(1-2)1107 126, 2010.
Nicholas G. Polson, James G. Scott, and Jesse Windle. Bayesian inference for logistic models using Pólya-Gamma Gamma latent variables. Journal of the American Statistical Association, 108(504):1339-13)) 2013.
Luke Y. Prince, Shahab Bakhtiari, Colleen J. Gillon, and Blake A. Richards. Parallel inference of hierarchical latent dynamics in two-photon calcium imaging of neuronal populations. bioRzie, 2021.
S. Joe Qin. An overview of subspace identification. Computers & Chemical Engineering, 30(10):1502-1613, 2006.
Nicholas A. Roy, Ji Hyun Bak, International Brain Laboratory, Athena Akrami, Carlos D. Brody, and Jonathan W. Pillow. Extracting the dynamics of behavior in sensory decision 109(4):597-611, 2021.
decision-making making experiments. Neuron,
Aaron Schein, Hanna Wallach, and Mingyuan Zhou. Poisson-Gamma dynamical systems. In Advances in Neural Information Processing Systems, volume 29, 2016.
W. Schultz, P. Dayan, and P. R. Montague. A neural substrate of prediction and reward. Science, 275(5306): 1593-1599, 1997.
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.
Sajid M Siddiqi, Byron Boots, and Geoffrey J Gordon. Reduced-rank hidden markov models. Proceedings of Machine Learning Research, 9:741-744 2010.
Anne C. Smith and Emery N. Brown. Estimating a state-space model from point process observations. Newral Computation, 15(5):965 991, 2003.
Brian M. Stocle. A modified EM algorithm for estimation in generalized mixed models. Biometrics, 52(4): 2295-1310 1996.
Hua Tang, Vincent D. Costa, Ramon Bartolo, and Bruno B. Averbeok. Differential coding of goals and actions in ventral and dorsal corticostriatal circuits during goal-directed behavior. Cell Reports, 38(1), 2022.
Adrian Valente, Srdjan Ostojic, and Jonathan W. Pillow. Probing the relationship between latent lincar dynamical systems and low-rank recurrent neural network models. Neural Computation, 34:1871-1892, 2022.
Peter van Overschee and Bart de Moor. N4SID: Subspace Algorithms for the Identification of Combined Searministic-Stocassttt Systems. Automatica, 30(1):75-93, 1994.
Peter van Overschee and Bart De Moor. Subspace Identification for Lineor Systems. Springer US, Boston, MA, 1996. ISBN 978- 11881--613-0061-00
Mats Viberg. Subspace-basee methods for the identification of linear time-invariant systems. Automatica, 31 (12):1835-1851, 1995.
Andrew Wagenmaker, Julian Katz-Samuels, and Kevin Jamicson. Experimental design for regret minimization in linear bandits. In Proceedings of The 24th International Conference on Artificial Intelligence and Statislics, pp. 3088 3096. PMLR, 2021.
Emily T. Wang, Marina Vannucci, Zulfi Hancef, Robert Moss, Vikram R. Rao, and Sharon Chiang. A Bayesian switching linear dynamical system for estimating seizure chronotypes. Proceedings of the National Academy of Sciences, 11944)::e2000822119 2022.
18

---

Published in Transactions on Machine Learning Research (07/2023)
Caleb Weinreb, Mohammed Abdal Monium Osman, Libby Zhang, Sherry Lin, Jonah Pearl, Sidharth Annapragada, Eli Conlin, Winthrop F. Gillis, Maya Jay, Ye Shackai, Alexander Mathis, Mackenzie Weygandt Mathis, Talmo Pereira, Scott W. Linderman, and Sandeep Robert Datta. Keypoint-MoSeq: parsing behavior by linking point tracking to pose dynamics. bioRriu, 2023.
Alexander B. Wiltschko, Tatsuya Tsukahara, Ayman Zeine, Rockwell Anyoha, Winthrop F. Gillis, Jeffrey E. Markowitz, Ralph E. Peterson, Jesse Katon, Malthew J. Johnson, and Sandeep Robert Datta. Revealing the structure of pharmacobehavioral space through motion sequencing. Nature Neuroscience, pp. 1-11, 2020.
Alexander B. Wiltschko, Matthew J. Johnson, Giuliano Turilli, Ralph E. Peterson, Jesse M. Katon, Stan L.. Pashkovski, Victoria E. Abraira, Ryan P. Adams, and Sandeep Robert Datta. Mapping sub-second structure in mouse behavior. Newron, 186))11121 2015.
Byron M. Yu, John P. Cunningham, Gopal Sant Santhanam, Stephen I. Ryu, Krishna V. Shenoy, and Maneesh Sahani. Gaussian-prose process factor analysis for low-dimensional single-trial analysis of neural population activity. Journal of Neurophysiology, 102(1):614 635, 2009.
Ke Yuan and Mahesan Niranjan. Estimating a state-space model from point process observations: note on convergence. Neural Computation, 22(8):1993 2001, 2010.
Feng Zhu, Harrison A. Grier, Raghav Tandon, Changjia Cai, Anjali Agarwal, Andrea Giovannucci, Matthew T. Kaufman, and Chethan Pandarinath. deep learning framework for inference of single-trial neural population dynamics from calcium imaging with subframe temporal resolution. Nature Neuroscience, (12):1724 1734, 2022.
David M Zoltowski, Jonathan W Pillow, and Scott W Linderman. A general recurrent state space framework for modeling neural dynamics during dccision-making. In Proceedings of the 37th International Conference on Machine Learning, 2020.
19

---

Published in Transactions on Machine Learning Research (07/202) 2023)
A Appendix
A.1 Simulated dataset parameters
Table A.1: Generative parameters for simulated datasets.
Simulated dataset parameters
Dataset
A
B
D
Q
q
p
m
A
eigenvalues
in
orthonormal
orthornormal
orthonormal
.11
1
[99.99]
scaled by
scaled by .1 orthornormal scaled by .1 orthonormal scaled by 10 orthonormal
scaled by .1 orthonormal .17
eigenvalues
orthonormal
[99,99]
scaled by orthonormal
scaled by orthonormal
C
eigenvalues in [55, .9]
.11
8
scaled by .1 orthonormal scaled by .1 orthonormal scaled by .1
scaled by .1 orthonormal
D
rotation matrix with = T/ /48
.00017
5
2
3
scaled by le4 orthonormal
scaled by .1
E
rotation matrix with = 7/2
orthonormal
.17
5
3
scaled by le4 orthornormal scaled by .1
scaled by .1 orthonormal
cigenvalucs in [99.99]
orthonormal scaled by scaled by .01
.11
scaled by .1
G
rotation matrix with 0 = / 400
scaled by .25
scaled by .2
.001I
1
A.2
Parameter recovery with non-unitized diagonal covariance
unitized
nonunitized
a
dynamics matrix A
b
loadings matrix C
input-output interactions
d
gain matrix G
matrix
0.08
0.04
8
avg. diff in EVs
subspace angle
1.0-
avg. IAl error
0.06
avg. IAI error
6
4
0.02
0.04
0.5
2
0.02
0
0.00
0.0
0.00
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
3.2 log1o(N)
4.2
5.2
log1o(N)
(og
log
Figure A.1: Stationary marginal covariance of & having non-unit diagonal covariance affects inference of C but mostly not other paramcters. (a-d) Average recovery error for bestLDS inferred parameters under unit-diagonal ("unitized") vs non-unil-diagonal ("nonunitized") regimes.
20

---

Published in Transactions on Machine Learning Research (07 (07/202) /023)
A.3
Estimator run time function of dataset size
a
dataset A
b
dataset
15.0
avg run time (s.
avg. run time (s.)
6-
12.5
10.0
2
7.5
3.2
4.2
5.2
3.2
5.2
(o
(o
Figure A.2: bestLDS run time as a function of dataset size.
A.4
Effect of moment conversion
moment conversion
without moment conversion
a
dynamics matrix A
b
loadings matrix C
C
input-output interactions
d
gain matrix G
matrix
0.08
0.04
avg. diff in EVs
subspace angle
1.0
avg. IAI error
0.06
avg. IAI error
4
0.02
0.04
0.5
0.02
0
0.00
0.0
0.00
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
3.2
4.2
5.2
log1o(N)
log1o(N)
log1o(M)
Figure A.3: bestLDS underperforms compared model that has direct access 2. (a-d) Average recovery error for bestLDS inferred parameters ("with moment conversion") vs. running N4SID directly on ("without moment conversion").
21

---

Published in Transactions on Machine Learning Research (07 (07/2023)
A.5
Comparison of recovered system parameters for bestLDS and bestLDS+EM
latent state dynamics
input-latent interactions
loadings
input-output interactions
bestLDS recovered parameters
1.00
0.4
1.0
0.6
0.50
0.2
0.3
0.5
0.00
0.0
0.0
0.0
-0.50
-0.2
-0.3
-0.5
-0.6
-1.00
-0,4
-1.0
bestLDS+EM recovered parameters
1.00
0.4
1.0
0.6
0.50
0.2
0.5
0.00
0.0
0.0
0.0
-0.50
-0.2
-0.3
-0.5
-0.6
-1.00
00.
-1.0
b
bestLDS
bestLDS+EM
A cues
prev. choice
prev. rew. choice
inhibition
0.4
0.10
0.02
0.05
-0.05
0.2
0.05
0.01
00.55
0.0
0.00
0.00
-0.25
0
20
40
0
20
40
20
40
20
40
trials
trials
trials
trials
Figure A.4: bestLDS recovered parameters closely match the inferred parameters using EM with bestLDS initializations. (a) The recovered parameters using bestLDS (top) and bestLDS+EM (bottom). (b) Impulse responses for each of the fitted systems. In each subplot, data has been simulated noiselessly with a unit input in the indicated dimension.
22


Reviews:
[{"review_id": "W8mxFWxOcE", "reviewer": "Reviewer_1ZBy", "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": "The article is well written and the literature revision is thorough. \n\nDespite the fact that the article is clearly presented, I had a tough time understanding the novel contribution of the manuscript. It is not entirely clear when the authors describe their setup and when they describe their proposal.\n\nThe authors proposed a model for time series---eq. (1)---together with a technique to analyse that model. Then, this proposal is tested on synthetic and real-world data. In this context, I can't identify a contribution to the machine learning community from the proposal of the authors. The fixed model together with the mechanism to analyse it do not constitute a novelty within ML.\n\nThe experimental validation is also far from convincing. The proposal is tested on synthetic data (Secs 4.1 and 4.2), where best LDS is compared against the ground truth or initialisation heuristics, but not other benchmarks. It is not clear to me if the known parameters used to produce realisations of the models are used in this experiments.\n\nIn the last experiment, though real-world data were considered, there are also no comparison: all claims are based on inspection of the results of best LDS rather than on a quantitative comparison against other methods. As a non-expert on what to expect from the mice experiment, the feeling I have after reading about the experiment is inconclusive: the method was applied and we obtained results, however, I wouldn't know how to assess them.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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": "The authors propose a spectral learning method for latent linear dynamical systems (LDS) with Bernoulli observations. The proposed method operates on transformations of the first and second moments of the time series. The proposal is validated experimentally. \n\n", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "Perhaps this article is out of the scope of TMLR, I would have expected a more general model formulation, a more principled methodology to analyse the model and a thorough experiment tal validation.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ck5pp92rez", "reviewer": "Reviewer_Ex51", "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 well-written. \n* The paper deals with a common setting in sequential data analysis, and some TMLR readers would appreciate this work. \n* Extensive experiments are performed to support the proposed method, including an experiment with a real-world behavioural dataset. \n\n## Weakness\n* The evaluation of the method could be made stronger. I would have preferred to see more synthetic experiments to investigate failure modes of the method. \n  * The synthetic data experiment seems to deal with a well-specified setting only. The method assumes the unit covariance diagonal of the marginal of the latent. What happens if this assumption is violated? \n  * For synthetic data, we can generate $z$. It might be helpful to compare against the exact N4SID to see the effect of the moment conversion. ", "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 work proposes an estimation method for a linear Gaussian state space model, where the continuous observation in the typical formulation is treated as latent, and the observation is instead given as a Bernoulli variable depending on the latent. \nA technical contribution of this work is to extend the N4SID method, a technique for continuous models, to this binary setting. \n\nSpecifically, the N4SID method requires the QR decomposition of the (unavailable) input-latent Hankel matrix. Authors proposed to use the Cholesky decomposition of the covariance of the input-latent Hankel matrix as the R matrix. \nMajor findings of the paper are the following: in a well-specified setting, the method (empirically) seems to yield a consistent estimate; the authors also show that this serves as a good initialisation for the EM algorithm. \n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please address the weakness comment above. ", "broader_impact_concerns": "The work is methodological, and I do not see any potential ethical implication. \n", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Fd1vizI4DO", "reviewer": "Reviewer_WsBP", "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": "Positives:\n1) The problem is an important one in domains such as neuroscience, and the proposed solution seems useful and practical.\n2) The paper is largely well written\n3) The model shows good performance in the experiments, which are relatively thorough\n\nNegatives: \n1) As far as I can tell, this paper is a fairly straightforward extension of Buesing et al from the Poisson case to the Bernoulli case. The authors here use the same principle: the observations Y are a transformation of Gaussian intermediate outputs Z, and learning the parameters of the system involves inferring the mean-covariance structure of Z from Y. The main technical problem here is that for Bernoulli variables, since the mean equals the mean square, the system of equations is underdetermined. The authors solution is just to set the diagonal of the covariance of Z to be one and solve for the remaining covariance variables. I'm not convinced this is a significant methodological innovation. The paper also allows inputs u to the linear-dynamical system, but dealing with this is also a fairly immediate application of Buesing et al.\n\n2) The proposed method only works for probit transformations, where the sysyem of nonlinear equations relating the moments of Y and Z can be written explicitly.\n\n3) Evaluation metrics seem ad hoc and not well motivated. Why is recovery of A measured using eigenvalues, and C using the angle between the subspaces? I don't think the paper clearly enough discusses the effect on estimating C of fixing the diagonal of the covariance of z.\n\nQuestions:\n1) Can we directly apply the estimator for the Poisson observations from Buesing et al, treating the Bernoulli outputs as outputs from a Poisson distribution with small rate parameter? Table 1 might be more stronger if the model is compared with that, right now there are no other baselines to compare performance with.\n\n2) It might help to mention how long it takes to find the bestLDS initialization, if we are going to claim that it help faster convergence. It would also be nice to see total time to convergence, including computing the initialization using bestLDS.\n\n3) I'm not really clear what the take-away for the real-world experiment is: the last paragraph before the discussion is very dense and hard to follow. Here too, it might be worth comparing with the Poisson LDS approach of Buesing et al.\n", "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 a spectral learning method for (possibly controlled) linear dynamical HMMs with Bernoulli observations, building on existing work for Gaussian/Poisson observations.The proposed method is fast, avoids hazards of local minima, and can be used to initialize EM algorithms for better solutions. The authors study their method on a number of synthetic tasks, and also apply it to analyze data from mice performing a decision-making task.\n\n", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1) I think the paper would benefit with a clearer discussion of the technical novelty, in its current form it does not seem like a significant advancement over Buesing et al. The authors might also consider comparing with this approach.\n\n2) I think the paper would benefit from a slightly clearer discussion of how the nonlinear systems of equations are solved. RIght now, it seems like this is a high-dimensional problem, but if I understand correctly, this is just a collection of uncoupled low-dimensional equations. It would also be nice if the paper also dwelt a little more on the situation when the equations don't have a solution, this is mentioned only very perfunctorily.\n\n3) I would like to see some theoretical and empirical discussion of the cost of solving the nonlinear moment conditions.\n\n4) I think \"for decision-making\" in the title is unnecessary and a bit confusing. The method also only applied for probit-transformations, which is not mentioned clearly upfront (in the abstract/introduction).\n", "broader_impact_concerns": "I see no major broader impact concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes spectral estimation for linear Gaussian state space models with binary observables.  It can be used as low-computation (and somehow interpretable) estimator, as well as initialization of the EM algorithm.  Experiments show good performance.

Reviewers originally mentioned the following pros and cons:

Pros:
- Well written
- Relevant problem setting
- Experiments support the proposed method.

Cons:
- Limited novelty (straight-forward extension of Buesing et al.)
- Weak experiments (lack of baselines in some experiments, miss-specified case)
- Evaluation metrics seem ad hoc.

The authors revised the manuscript, addressing most of the concerns.  The novelty is at least clear.  Although it's not innovation, it matches the TMLR criterion.
In the revision, baseline methods are added, and the miss-specified cases are discussed.  The authors' arguments on the evaluation metrics are convincing.

Overall, the paper has been significantly improved after revision, and the revised paper satisfies the TMLR criteria for acceptance.

I have a question, which is the reason why the decision is conditional.  Please make this point clear with minor revision:
In Table 1, bestLDS indeed gives the highest evidence.  My question is why it is even better than "true".  I can't find clear explanation about what "true" means, and I'd assume that it is a kind of oracle, which gives an upper-bounds of performance of practical methods.  If this would be the case, the reason why bestLDS gives higher evidence should be discussed.
If not and "true" is just another baseline method, which bestLDS outperforms, then the method "true" should be clearly explained, and what makes the "true" method weak even though it uses the true parameters should be discussed.
Decision: Accept (poster)